Olhando e comparando outros modelos para descrever propriedades de redes de interações

01

Continuando o que estávamos vendo no post anterior, vamos olhar novamente somente a rede do exemplo 4, acima.
Aqui temos bastante interações e desde polinizadores muito generalistas a poucos generalistas.
Apesar de essa rede ter sido gerada por acaso, é muito interessante o fato que as redes biológicas vão ser fruto das estratégias das espécies que a compõe, e ao conseguir descrever essa rede e falar ela é assim, podemos passar a pensar quais forças evolutivas estão agindo como um todo na comunidade. Claro que a pressão evolutiva vai ser exercida nos indivíduos resultando no que é o jeitão comum da espécie e não vemos esse processo, mas aqui estamos olhando para o resultado dele.

Então continuando, nos ajustamos o modelo exponencial para a “escadinha” que se forma ao contarmos quantas espécies tem um determinado número de interações.

Lembrando que esse é a função que a gente fala o degrau da escadinhas e a função da a altura do degrau:

 p(k) = b \cdot exp(-\gamma \cdot k)

E ajustando o modelo no R usando a função nls temos como resultado

Formula: P ~ b * exp(-gamma * k) Parameters: Estimate Std. Error t value Pr(>|t|) gamma 0.14839 0.01197 12.39 2.15e-07 *** b 1.26079 0.06602 19.10 3.37e-09 *** ---

E olhando esses números no gráfico temos

02

Temos nossos 2 parâmetros, que o mais interessante é o gamma (lembre-se que gamma é essa letra grega \gamma), que é a inclinação da escadinha. Mais interessante no sentido que ele esta informando muita coisa sobre a comunidade, ele que esta informando se temos mais especialistas, menos especialistas, mais generalistas, quantos de cada tipo, isso fica impresso naquele gráfico.

Existem ainda duas possibilidades como proposto pelo Jordano (2003).
A primeira é a power law, aqui o power é de potencia, um número elevado ao outro, que vemos na formula da função que é o modelo.

 p(k) = b \cdot k^{-\gamma}

Novamente ajustando no R temos:

Formula: P ~ b * k^(-gamma) Parameters: Estimate Std. Error t value Pr(>|t|) gamma 0.51585 0.09301 5.546 0.000245 *** b 1.18185 0.14298 8.266 8.83e-06 *** ---

E olhando no gráfico temos

03

Uma coisa legal de se pensar em relação a esse power law é que ele não enverga muito bem, ele tem uma capacidade de previsão muito pior que o modelo exponencial aqui nesse exemplo, mas talvez ele seja muito mais adequado para descrever vários generalistas, ou uma predominância de generalismo, já que não “envergar” bem o modelo quer dizer que os degraus (lembre-se do gráfico de barras do post anterior) vão continuar do mesmo tamanho, ou seja o mesmo número de espécies varias vezes com o mesmo número de interações, generalistas.

Agora vamos olhar a terceira opção, Truncated Power Law.
Bem se é power law, é algo a potência do outro, agora com esse “Truncated” que eu não sei o que faz direito. Mas vamos olhar ele melhor.

 p(k) = b \cdot (k^{-\gamma}) \cdot exp(-k/kx)

E temos um novo parâmetro, que deve ser o tal do “Truncated”, e vamos la e ajustamos ele no R e temos…

Formula: P ~ b * (k^(-gamma)) * exp(-k/kx) Parameters: Estimate Std. Error t value Pr(>|t|) gamma -0.29097 0.08258 -3.523 0.00648 ** b 1.22771 0.04615 26.605 7.24e-10 *** kx 4.48568 0.47416 9.460 5.67e-06 *** ---

Agora temos um gamma negativo, mas olhando o gráfico que é mais simples de ver as coisas que esses números, de inicio vemos que…

04

Ele parece que tem um ajuste melhor que todos os outros modelos, a linha passa bem pertinho dos pontos, mas por outro lado ele tem um parâmetro a mais. Parâmetros deixam a gente mexer na forma da linha logo, quanto mais parâmetros sempre a gente vai ter a possibilidade de um melhor ajuste, lembrando que possibilidade num é garantia.

Mas a gente pode usar o AIC para verificar a qualidade dos ajustes, e colocar em número isso que estamos olhando nesses gráficos.
Vamos fazer então um gráfico com os 3 modelos e olhar os respectivos AICs

05

#Exponencial AIC(m1.exe4) [1] -27.11152 #Power Law (m2.exe4) [1] -5.733472 #Truncated Power Law AIC(m3.exe4) [1] -36.24854

Bem o aic é uma medida de ajuste do modelo, o quanto a linha fica bem pertinho dos pontinhos do gráfico mais ou menos. Quanto menor é melhor, e ao olhar-mos aqui vemos o seguinte, olhando no gráfico da para ver facilmente que o Power Law é o pior nesse caso, talvez com mais generalistas ele fosse melhor, mas para esses dados que temos aqui ele é ruim, ou seja tem o Maior AIC, -5.

Agora olhando no gráfico, o modelo Exponencial e o Truncated Power Law tem um bom ajuste, e tem valores baixos de AIC, -27 e -36 respectivamente. No entanto podemos pensar um pouco ainda sobre simplesmente usar a regra de pegar o menor, o modelo truncado.
Bem o modelo truncado tem um parâmetro a mais. No exponencial as coisas são mais simples de imaginar, o b vai acabar empurrando a linha mais para cima ou mais para baixo no gráfico, e o gamma vai ser o quão íngreme é a subida da curva. Agora esse modelo truncando tem um outro parâmetro que deixa tudo mais difícil de pensar, ele parece que torce um pouco a curva, uma coisa a se fazer é ficar mudando valores de kx, b e gamma e tentar entender como tudo funciona se você como eu não conseguir imaginar como a curva vai ficar.

Outra coisa é que qualquer um ponto a mais no gráfico muda esse valor de AIC, qualquer interação que vimos a mais, ou que não conseguimos ver pode trocar quem é o modelo com ajuste melhor. Apesar da diferença dos dois modelos para o power law serem grandes, entre eles talvez esse 10 de diferença de AIC seja espurico. Mas podemos tentar fazer um grau de confiança se realmente esses modelos são diferentes, pois se ele não forem diferentes é mais fácil a gente ficar com o mais simples, que é mais fácil de pensar, sempre lembramos da navalha de Occam nesses momentos. Usando a estatística f vemos o seguinte ao comparar os modelos.

anova(m1.exe4,m3.exe4,test="f") Analysis of Variance Table Model 1: P ~ b * exp(-gamma * k) Model 2: P ~ b * (k^(-gamma)) * exp(-k/kx) Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F) 1 10 0.044500 2 9 0.017591 1 0.026909 13.767 0.004842 ** ---

E parece que os modelos são diferentes, o modelo truncando é um melhor preditor e ficamos com ele. Mas esse teste f para comparar os modelos não é muito perfeito. Principalmente porque comparar modelo não linear é muito complicado. Nos fóruns as pessoas são bastante relutantes em usar essas coisas com modelos não lineares, como os que estamos ajustando aqui, mas a principio não estamos vendo nada que não esparavamos ao avaliar visualmente os gráficos e a analise dar um resultado que você já espera ao olhar o gráfico é uma coisa boa.

Então assim finalmente resolvemos usar o modelo truncado para descrever essa comunidade. Se tivessemos mais comunidades poderíamos comparar os parâmetros delas agora quanto a essa propriedade descrita por esse gamma do modelo e começar a pensar o que esta acontecendo, porque que esses valores fiquem assim, ou se tiver alguma diferença, o que poderia gerar essa diferença em comunidades.

Bem um ultimo detalhe é que se olhar o gráfico do bipartite ele fica diferente do que fizemos aqui.

06

Talvez para uma comunidade muito maior que a desse exemplo, o gráfico ficaria muito grande do jeito que eu fiz, então ele usa a escala logarítmica para plotar o gráfico. Quer dizer que ele calcula o log do valor de x e do valor de y para plotar. Eu achei mais ilustrativo usar os dados brutos para explorar aqui o que esta acontecendo, mas se a gente fizer questão de ver a mesma coisa é so colocar log nos valores de x e y

07

#Semente para os resultados serem reproduzíveis
set.seed(13)
###################################################
#Função para gerar rede de interação baseado no modelo de Erdős–Rényi
#Creditos a Edwin Grappin, vi o post no R-bloggers
#http://probaperception.blogspot.com.br/2012/12/have-you-tried-to-understand-your.html
###################################################
generateER = function(n = 100, p = 0.5){
map = diag(rep(1, n))
link = rbinom(n*(n-1)/2, 1,p)
t = 1
for(j in 2:n){
for(i in 1:(j-1)){
map[i,j] = link[t]
t = t + 1
}
}
dimnames(map)<-list(paste("Pol",1:n),paste("Sp",1:n))
return(map)
}
###########################
#Minha função para transformar matriz em lista para o igraph
##############################
mat2lis<-function(matriz) {
lista<-dimnames(matriz)
dados<-data.frame(Pol=rep(lista[[1]],length(lista[[1]])),
SP=rep(lista[[2]],each=length(lista[[2]])),
int=c(matriz))
return(dados[dados$int==1,1:2])
}
############################
library(igraph)
 
#Gerando exemplos
exe4<-generateER(n=15,p=0.8)
 
#Transformando na forma de dados que o igraph usa
exe4.l<-mat2lis(exe4)
exe4.g <-graph.data.frame(exe4.l, directed=F)
 
#Plots da rede
#01
plot(exe4.g,layout=layout.fruchterman.reingold,
vertex.color=rep(c("blue","green"),each=15),
vertex.label=c(paste("P",1:15),paste("F",1:15)),
vertex.shape=rep(c("square","circle"),each=15),
edge.color="black",edge.width=3,vertex.label.cex=0.5,
vertex.label.color="black"
)
title("Exemplo 4")
 
#Para ver os modelos consulte…
#Jordano, P., Bascompte, J. and Olesen, J. M. (2003)
#Invariant properties in coevolutionary networks of
#plant-animal interactions. Ecology Letters 6, 69–81
 
#Exemplo 4
k<-sum(exe4)
pol.n<-rowSums(exe4)
P.pol<-sapply(sort(unique(pol.n)), function(x) sum(pol.n>=x))
dados.pol4 <- cbind.data.frame(k = sort(unique(pol.n)), P = P.pol/max(P.pol))
m1.exe4<-nls(P~b*exp(-gamma*k),start=list(gamma = 0, b = 1),data=dados.pol4)
m2.exe4<-nls(P~b*k^(-gamma),start=list(gamma = 0, b = 1),data=dados.pol4)
m3.exe4<-nls(P~b*(k^(-gamma))*exp(-k/kx),start=list(gamma = 0, b = 1,kx=5),data=dados.pol4)
 
#Grafico modelo exponencial
plot(dados.pol4[,2]~dados.pol4[,1],pch=19,cex=2,ylab="Distribuição Acumulada",
xlab="K",frame.plot=F,ylim=c(0,1),xlim=c(0,18))
curve(coef(m1.exe4)[2]*exp(-coef(m1.exe4)[1]*x),1,15,add=T,lwd=2,lty=2,
cex=1.5,col="red")
legend("topright",lwd=2,lty=2,col=c("red"),
legend=c("Exponential"))
#grafico do modelo power law
plot(dados.pol4[,2]~dados.pol4[,1],pch=19,cex=2,ylab="Distribuição Acumulada",
xlab="K",frame.plot=F,ylim=c(0,1),xlim=c(0,18))
curve(coef(m2.exe4)[2]*x^-coef(m2.exe4)[1],1,15,add=T,lwd=2,lty=2,
cex=1.5,col="blue")
legend("topright",lwd=2,lty=2,col=c("blue"),
legend=c("Power Law"))
#grafico do modelo truncado
plot(dados.pol4[,2]~dados.pol4[,1],pch=19,cex=2,ylab="Distribuição Acumulada",
xlab="K",frame.plot=F,ylim=c(0,1),xlim=c(0,18))
curve(coef(m3.exe4)[2]*(x^(-coef(m3.exe4)[1]))*exp(-x/coef(m3.exe4)[3]),1,15,add=T,lwd=2,lty=2,
cex=1.5,col="green")
legend("topright",lwd=2,lty=2,col=c("green"),
legend=c("Truncated Power Law"))
 
#os 3 modelos juntos
plot(dados.pol4[,2]~dados.pol4[,1],pch=19,cex=2,ylab="Distribuição Acumulada",
xlab="K",frame.plot=F,ylim=c(0,1),xlim=c(0,18))
curve(coef(m1.exe4)[2]*exp(-coef(m1.exe4)[1]*x),1,15,add=T,lwd=2,lty=2,
cex=1.5,col="red")
curve(coef(m2.exe4)[2]*x^-coef(m2.exe4)[1],1,15,add=T,lwd=2,lty=2,
cex=1.5,col="blue")
curve(coef(m3.exe4)[2]*(x^(-coef(m3.exe4)[1]))*exp(-x/coef(m3.exe4)[3]),1,15,add=T,lwd=2,lty=2,
cex=1.5,col="green")
legend("topright",lwd=2,lty=2,col=c("red","blue","green"),
legend=c("Exponential","Power Law","Truncated Power Law"))
 
#resultados dos modelos e seus AICs
summary(m1.exe4)
summary(m2.exe4)
summary(m3.exe4)
AIC(m1.exe4)
AIC(m2.exe4)
AIC(m3.exe4)
 
#comparação dos modelos, leia nos foruns primeiros
#sobre comparar modelos não lineares antes de confiar
#nesse resultado
anova(m1.exe4,m3.exe4,test="f")
 
#Grafico e ajuste usando o bipartite para o exemplo 4
library(bipartite)
degreedistr(exe4, plot.it=T)
 
#grafico como apresentado no bipartite
plot(dados.pol4[,2]~dados.pol4[,1],pch=19,cex=2,ylab="Distribuição Acumulada",
xlab="K",frame.plot=F,log="xy")
abline(h=1,lty=2)
curve(coef(m1.exe4)[2]*exp(-coef(m1.exe4)[1]*x),1,15,add=T,lwd=2,lty=2,
cex=1.5,col="red",log="xy")
curve(coef(m2.exe4)[2]*x^-coef(m2.exe4)[1],1,15,add=T,lwd=2,lty=2,
cex=1.5,col="blue",log="xy")
curve(coef(m3.exe4)[2]*(x^(-coef(m3.exe4)[1]))*exp(-x/coef(m3.exe4)[3]),1,15,add=T,lwd=2,lty=2,
cex=1.5,col="green",log="xy")
legend(2,0.3,lwd=2,lty=2,col=c("red","blue","green"),
legend=c("Exponential","Power Law","Truncated Power Law"))
 
#########################################

Usando o AIC e Maximum likelihood, um pequeno exemplo.

Ha muitas armadilhas para interpretar resultados estatísticos, ainda mais se tratando de valores p como ja visto em um post que escrevi aqui.
Mas valores p não são a única forma de se obter conclusões, ou pensar em relação a análise de dados.
Como a gente já viu em outro post anterior, podemos pensar no ajuste do modelo aos dados.

Isso tem raízes na teoria da informação eu acredito, assim a máxima verosimilhança e o critério de informação de Akaike podem ser usadas como medidas de de o quão bom é o ajuste do modelo aos dados.

Então podemos pensar ao invés do famoso “Existe diferença significativa”, qual o modelo que melhor explica nossos dados, dentro de uma gama de vários modelos que podemos produzir e então utilizar o mais adequado que veio de todos os modelos que produzimos.

Vou tentar ilustrar com um exemplo de como usar esse tipo de abordagem funciona, pelo menos para esses simples exemplos.

Então vamos supor uma uma regressão, onde temos o modelo “real” que é resposta=a+b*preditor. Simulando isso para um a=5 e b=2 temos:

Certo, vamos então ajustar 2 modelos.
Um modelo vai ser a regressão normal, resposta = a + b * preditor e o segundo modelo vai ser resposta = a , ou seja o preditor não prevê nada (o que sabemos que não é verdade).

Gerado os dados, podemos usar o comando logLik() para calcular um valor pro “quão bom esse modelo é”. Minha terminologia deve ser meio ruim, o primeiro lugar a verificar definições mais corretas é na documentação da própria função usando o comando ?logLik() no R.
Mas continuando. Faremos isso para os 2 modelos que ajustamos:

> logLik(modelo1) 'log Lik.' -46.87914 (df=3) > logLik(modelo2) 'log Lik.' -78.03321 (df=2)

Então temos que o valor do modelo que sabemos que é o correto (já que simulamos ele, mas a natureza nunca nos conta isso…) é de -46, enquanto o modelo mais simples (e errado) foi de -78, vemos então que quanto melhor o ajuste (e teoricamente melhor o modelo) temos um valor maior de logLik.

Agora o AIC é mais ou menos esse valor, só que ele da uma penalidade por parâmetro estimado.
Como assim parâmetro? Assim, o modelo 1 tem a+bx, então precisamos de um parâmetro, o “a”, dois parâmetros, o “b” e a quantidade de desvio do modelo (deviance), variância. Olhem la em cima que o df=3, ou seja o modelo tem esses 3 ingredientes.

Agora o modelo 2 tem menos ingredientes, ele só tem um “a” e desvio, a variação dos dados, olhem como o df=2. Ou seja, se ele explicar bem, mesmo que um pouco pior que o modelo 1, ele deveria ganhar uma “colher de chá”, já que ele tenta explicar os mesmo dados com menos parâmetros, de forma mais simples.

Esse é o aic, que podemos calcular usando o comando extractAIC(), para ver certinho a forma como é calculado, entre outras informações, use o comando ?extractAIC().

Mas utilizando ele no nosso exemplo temos:

> extractAIC(modelo1) [1] 2.00000 12.62196 > extractAIC(modelo2) [1] 1.0000 72.9301

Agora o modelo 1, com 2 parâmetros(a e b, variação todo modelo tem que ter senão não seria probabilístico então não contamos aqui) tem um valor baixo de 12 enquanto o modelo 2 (que sabemos que é o errado) tem 1 parâmetro e o valor de 72 de AIC, então valores menores devem ser bons. Somos inclinados a usar o modelo 1 que tem os 2 parâmetros, intercepto e inclinação.

Mas podemos ainda fazer um teste F, vendo a razão entre os logLik para saber se os modelos são diferentes.

> anova(modelo1,modelo2) Analysis of Variance Table Model 1: resposta ~ preditor Model 2: resposta ~ 1 Res.Df RSS Df Sum of Sq F Pr(>F) 1 28 39.99 2 29 319.11 -1 -279.12 195.44 3.745e-14 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

E vemos o que? Que existe diferença entre os modelos, então ficamos com aquele com mais parâmetros, o modelo 1, o mais complexo, não podemos abandonar ele pelo mais simples, ja que ele explica muita coisa que o modelo 2 mais simples não deu conta de explicar.

Podemos ainda fazer isso usando o comando step(), que vai pegar o modelo inicial que você fornecer, e se vc mandar ele simplificar com o argumento direction =”backward”, ele vai tirando um por um dos parâmetros, e ve se o modelo fica diferente, se continuar explicando bem os dados, ele descarta esse parâmetro, so deixar o trace=1 para ver o processo, como o modelo mais complexo aqui, o modelo 1 tem apenas 2 parâmetros, ele só faz isso uma vez, mas se você tiver muitos parâmetros ele vai fazendo isso pra você.
Então temos:

ps(trace=1 e nesse caso o direction=”backward” são os default da função)

> mod.final1<-step(modelo1) Start: AIC=12.62 resposta ~ preditor Df Sum of Sq RSS AIC 39.99 12.622 - preditor 1 279.12 319.11 72.930

ele tira um parâmetro, o modelo cresce muito o valor de AIC ai ele fala, não da pra tirar esse parâmetro senão o modelo vai ficar ruim, então ele deixa o modelo como está.

> mod.final1 Call: lm(formula = resposta ~ preditor) Coefficients: (Intercept) preditor 4.696 2.069

Até agora nos vimos um caso onde havia uma relação, ou seja havia uma inclinação, o b tinha um valor de 2 nos dados que simulamos, mas e o contrario será que funciona?
Então simulamos um caso onde não existe relação:

Podemos começar como no exemplo anterior, vamos chamar agora o modelo de regressão de modelo 3, e o modelo que é somente uma média de modelo 4, sabemos (porque simulamos) que o modelo 4 é o correto, e o 3 é errado, então olhando os valores de logLik.

> logLik(modelo3) 'log Lik.' -73.34481 (df=3) > logLik(modelo4) 'log Lik.' -73.34932 (df=2)

Agora eles são praticamente os mesmo, os dois modelos parecem explicar bem os dados, a diferença é que o modelo 3 intercepto, inclinação e variação e o modelo 4 usa somente intercepto e variação, ele explica o mesmo sendo mais “económico”.

Utilizando o AIC vemos:

> extractAIC(modelo3) [1] 2.00000 65.55331 > extractAIC(modelo4) [1] 1.00000 63.56234

Aha, agora o modelo 3 com 2 parâmetro tem o valor de AIC maior. Ou seja apesar de eles explicarem bem os dados, ele é penalizado por ter um parâmetro a mais e estar explicando a mesma coisa do modelo 4.

Se fizermos um teste F entre os modelos agora nos vemos que:

> anova(modelo3,modelo4) Analysis of Variance Table Model 1: resposta2 ~ preditor2 Model 2: resposta2 ~ 1 Res.Df RSS Df Sum of Sq F Pr(>F) 1 28 233.45 2 29 233.52 -1 -0.070269 0.0084 0.9275

Eles são iguais, logo a gente escolhe o mais simples, se ambos explicam igualmente os dados. Qual a lógica disso? Uma introdução rápida a ideia da navalha de Occlam, imagine que cada parâmetro vai explodir, mas ele também é a chave para a verdade, se você explodir não vai dar para ver a verdade, mas parâmetros podem valer muito, então a gente pesa, será que eu fico com 2 parâmetros e tenho 2 chances de explodir ou como eu sei que um não serve para muita coisa jogo ele fora, e fico com uma chance de ver a verdade e uma chance somente de explodir. Algo mais ou menos assim.

Podemos fazer esse processo com o comando step() também e ficamos com:

> mod.final2<-step(modelo3) Start: AIC=65.55 resposta2 ~ preditor2 Df Sum of Sq RSS AIC - preditor2 1 0.070269 233.52 63.562 233.45 65.553 Step: AIC=63.56 resposta2 ~ 1 > mod.final2 Call: lm(formula = resposta2 ~ 1) Coefficients: (Intercept) 19.69

Agora ele viu que o modelo fica igual, quase não aumenta o AIC então ele tira a inclinação do modelo e vai para o próximo passo, como já só tem um parâmetro no próximo passo ele para, e nosso modelo final selecionado pelo step() é somente intercepto, igual ao modelo 4.

Bem essa é somente uma ideia de como funciona o AIC, mas é preciso cuidado já que os modelos usados nos exemplos são extremamente simplistas. Com modelos mais complexos as coisas precisam de mais pensamento e reflexão, e nem step() ainda é melhor que um pouco de raciocínio, sempre temos que olhar se o resultado ou respostas que estamos tendo fazem sentido com a reliadade. Mas a partir daqui já da para ter uma ideia do que se trata AIC, logLik e compania, já que esses valores saem junto com modelos que ajustamos na maioria das funções no R e em outros programas estatísticos.
Outra coisa a notar, é que a conclusão sobre o que é importante aqui é a mesma usando valores p ou AIC/logLik, so dar um summary() nos modelos e veremos como no final a nossa conclusão sobre o que é importante é a mesma. Temos muitas possibilidades hoje de chegar a conclusões, mas é importante ter pelo menos uma ideia de todas para facilitar a leitura de artigos e saber as possibilidades existentes para resolver nossos próprios problemas:

set.seed(5)
#exemplo1
#simulação
preditor<-runif(30,5,10)
resposta<-rnorm(30,5+2*preditor)
 
#grafico
plot(resposta~preditor,pch=19,cex=1.6,xlab="Preditor",ylab="Resposta")
 
#ajuste dos modelos
modelo1<-lm(resposta~preditor)
modelo2<-lm(resposta~1)
 
#Qualidade do ajuste
logLik(modelo1)
logLik(modelo2)
 
extractAIC(modelo1)
extractAIC(modelo2)
 
anova(modelo1,modelo2)
 
#seleção automatica de modelo
mod.final1<-step(modelo1)
mod.final1
 
#compare com os valores p
summary(modelo1)
 
#exemplo2
#simulação, note que eu uso a media e desvio da resposta anterior
preditor2<-runif(30,5,10)
resposta2<-rnorm(30,mean(resposta),sd(resposta))
 
#grafico
plot(resposta2~preditor2,pch=19,cex=1.6,xlab="Preditor",ylab="Resposta")
 
#ajuste dos modelos
modelo3<-lm(resposta2~preditor2)
modelo4<-lm(resposta2~1)
 
#avaliando a qualidade do ajuste
logLik(modelo3)
logLik(modelo4)
 
extractAIC(modelo3)
extractAIC(modelo4)
 
anova(modelo3,modelo4)
 
#seleção automatica de modelo
mod.final2<-step(modelo3)
mod.final2
 
#compare com os valores p
summary(modelo3)