Generalized least squares (GLS) – Um exemplo com auto-correlação espacial

Vamos dar uma olhada como funciona esse tal de Generalized least squares (GLS), pra começar, vamos ver um exemplo, que está no livro do Michael J. Crawley, o R book segunda edição.

Bem esse livro é bem legal, mas eu sempre olhei ele mais por capítulos, da pra ler um capítulo separado de forma relativamente tranquila, pelo menos os que eu olhei.
Bem esse está no capítulo de estatística espacial. Como a gente ja estava falando de correlação filogenética, vamos ver o uso do gls sem ser com a correlação filogenética primeiro.

Podemos abrir os dados que ele usa no exemplo direto da internet.

1
spatialdata <- read.table("http://www.bio.ic.ac.uk/research/mjcraw/therbook/data/spatialdata.txt",header=T)

E temos os seguintes dados.

str(spatialdata) > 'data.frame': 224 obs. of 5 variables: $ Block : int 1 1 1 1 1 1 1 1 1 1 ... $ variety : Factor w/ 56 levels "ARAPAHOE","BRULE",..: 12 2 50 7 1 14 15 16 4 52 ... $ yield : num 29.2 31.6 35 30.1 33 ... $ latitude : num 4.3 4.3 4.3 4.3 4.3 4.3 4.3 8.6 8.6 8.6 ... $ longitude: num 19.2 20.4 21.6 22.8 24 25.2 26.4 1.2 2.4 3.6 ...

São testes com 56 variedades de trigo, mas num experimento em escala geográfica, com fazendas distribuídas ao longo de um gradiente de latitude e longitude. O yield é a resposta, a produtividade, variety é um fator, com 56 níveis que é o nome da variedade, block são os blocos de amostras.

Bem a gente tem 4 amostras por variedade de trigo

> table(spatialdata$variety) ARAPAHOE BRULE BUCKSKIN CENTURA CENTURK78 CHEYENNE CODY 4 4 4 4 4 4 4 COLT GAGE HOMESTEAD KS831374 LANCER LANCOTA NE83404 4 4 4 4 4 4 4 NE83406 NE83407 NE83432 NE83498 NE83T12 NE84557 NE85556 4 4 4 4 4 4 4 NE85623 NE86482 NE86501 NE86503 NE86507 NE86509 NE86527 4 4 4 4 4 4 4 NE86582 NE86606 NE86607 NE86T666 NE87403 NE87408 NE87409 4 4 4 4 4 4 4 NE87446 NE87451 NE87457 NE87463 NE87499 NE87512 NE87513 4 4 4 4 4 4 4 NE87522 NE87612 NE87613 NE87615 NE87619 NE87627 NORKAN 4 4 4 4 4 4 4 REDLAND ROUGHRIDER SCOUT66 SIOUXLAND TAM107 TAM200 VONA 4 4 4 4 4 4 4

E podemos começar fazendo uma figura, to tipo boxplot, para ver como é a distribuição dessas amostras.

01

Bem, aparentemente, temos uma diferença entre algumas variedades pelo menos, com uma inspeção visual, a gente consegue ver que tem variedades de baixa produtividade e variedades de alta produtividade.

Mas ao realizar uma análise de variância, para avaliar esses dados, vemos o seguinte.

> model1 <- aov(yield~variety,data=spatialdata) > summary(model1) Df Sum Sq Mean Sq F value Pr(>F) variety 55 2387 43.41 0.73 0.912 Residuals 168 9990 59.47

Não esperamos uma diferença na produtividade dessa variedades, o que é estranho, dado a aparência do gráfico, mas veja que estamos quebrando fortemente uma das premissas da analise de variância aqui, que é a da homogeneidade das variâncias, (Homoscedasticity). Tudo bem quebrar um pouquinho, mas aqui acho que foi demais.

Se a gente olhar a produtividade por bloco do experimento

> tapply(spatialdata$yield,spatialdata$Block,mean) 1 2 3 4 27.57500 28.81091 24.42589 21.42807

A gente vai ver que o bloco 2 é quase 10 unidades mais que o bloco 4, existe uma variação grande ai entre os blocos. Então talvez se tentarmos tirar a variação desses blocos, algo como uma anova pareada, de certo, e tentamos isso.

> Block <- factor(spatialdata$Block) > model2 <- aov(yield~Block+variety+Error(Block),data=spatialdata) > summary(model2) Error: Block Df Sum Sq Mean Sq Block 1 1469 1469 Error: Within Df Sum Sq Mean Sq F value Pr(>F) variety 55 2393 43.51 0.853 0.75 Residuals 167 8516 50.99

E ainda continuamos no mesmo problema, mas veja que o valor F ja inflou um pouquinho, e o valor p diminuiu. Como temos as coordenadas de onde o experimento foi feito, a coordenada de cada amostras, vamos dar uma olhada nelas em relação a produtividade.

02

Veja que a gente tem o que parece ser uma tendência, ao longo do espaço aqui.

Podemos adicionar ela ao modelo, adicionando as variáveis latitude e longitude

> model3 <- aov(yield~Block+variety+latitude+longitude,data=spatialdata) > summary(model3) Df Sum Sq Mean Sq F value Pr(>F) Block 1 1469 1469.1 43.730 5.00e-10 *** variety 55 2393 43.5 1.295 0.109 latitude 1 692 691.8 20.594 1.09e-05 *** longitude 1 2281 2280.8 67.891 5.09e-14 *** Residuals 165 5543 33.6 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

E elas são bastante significativas, o valor F vai la em cima, e agora variedade quase está se mostrando com uma diferença. Mas aqui a gente ainda está pecando em uma coisa, a gente está assumindo que a variação na latitude e longitude é linear, e na figura la em cima da latitude e longitude, aquilo parece tudo menos linear.

Agora vamos fazer uma figura, mostrando o resultado do experimento ao longo do espaço.

03

Aqui, a cor representa a variedade, os eixos x e y são a latitude e longitude, e o tamanho da bolinha é proporcional a produtividade daquela amostras, além disso, do lado de cada amostra eu coloquei o número do bloco. Pra mim essa figura é muito reveladora.
Primeiro veja que no canto inferior direito, algo acontece, a produtividade la é muito baixa, independente da variedade estudada. Mas como o bloco é delimitado de forma constante, numa sequência, a gente da o azar do o que quer que ocorre naquela região, afeta parte do bloco 3 e do bloco 4, e como de cima pra baixo a gente tem uma diminuição, e da esquerda pra direita também, pelo menos no canto, bloco, latitude e longitude são significativos, mas não em relações lineares, mas quando mais próximo ali daquele ponto, a produtividade é menor, e longe dali as coisas melhoram.

E finalmente decidimos partir para o GLS.

Começamos assim:

> model4 <- gls(yield~variety,spatialdata) > anova(model4) Denom. DF: 168 numDF F-value p-value (Intercept) 1 2454.621 <.0001 variety 55 0.730 0.9119

Veja que aqui, na verdade não fizemos nada de diferente do primeiro modelo, veja que todos os valores são iguais. Isso porque se a gente olhar o help do gls (?gls), o tipo de correlação default é o uncorrelated.

Que vai ser algo assim:

   \left( \begin{array}{ccc}  s^2 & 0 & 0 \\  0 & s^2 & 0 \\  0 & 0 & s^2 \end{array} \right)

Lembra-se aqui, quando a gente colocava aquela matriz para gerar valores da distribuição multivariada normal? Então essa é a tal da matriz de variância-covariância. Tudo se encaixa não?

No caso, quando a gente altera aqueles valores de zero para outra coisa, estamos colocando a correlação entre os valores. Essa é uma matriz quadrada, e o tamanho dela é o tamanho das amostras.

Bem, agora chegar a essa matriz é difícil, a está a alguns posts, aqui e aqui, batalhando para entender como é o cara na filogenia, mas a algum tempo, olhamos os dados do Jon Snow, e aqui demos uma olhada em autocorrelação espacial.

Bem, para aquilo que começamos a fazer la, existem funções prontas, uma classes de funções na verdade, que no R são as corStruct, que descrevem essa matriz de covariância, e fornecem os dados para o gls levar em conta a correlação entre os dados, seja la qual for.

Primeiro vamos olhar a função de correlação do variograma. Não vou explicar ela aqui, porque nem eu entendo direito, mas mais ou menos como fizemos nos dados do Snow, vemos a correlação dos dados de acordo com distância, calculada a partir da latitude e longitude.

04

Vemos uma forte relação da distância, basicamente, só olhar naquele gráfico la em cima, tem uma área de produtividade baixa no mapa. Então vamos incorporar isso no modelo, vamos falar, olha todo mundo por aqui vai ter valores meio baixos, la pra cima todo mundo tem valor mais alto, e assim vai.

> model5 <- update(model4,corr=corSpher(c(28,0.2),form=~latitude+longitude,nugget=T)) > anova(model5) Denom. DF: 168 numDF F-value p-value (Intercept) 1 82.46320 <.0001 variety 55 1.87472 0.0012

Fazemos um update do modelo anterior (com a estrutura de covariância independente) para uma estrutura de covariância corSpher (?corSpher para entender melhor) e agora vemos uma diferença entre as variedades.

Mas essa não é a única estrutura de covariância possível, podemos testar outras.

anova(model6) > Denom. DF: 168 numDF F-value p-value (Intercept) 1 30.399419 <.0001 variety 55 1.850939 0.0015

Mas veja que a conclusão final ficou igual. o que é bom.
Podemos então comparar o modelo 5 com o 6 para ver qual explicou melhor os dados

> anova(model5,model6) Model df AIC BIC logLik model5 1 59 1185.863 1370.177 -533.9315 model6 2 59 1183.278 1367.592 -532.6389

Veja que o modelo 6 tem um valor menor de AIC, assim perguntamos, vale a pena levar em conta a auto-correlação espacial? (Comparando com o modelo 4, que não tem nenhuma estrutura de covariância nos dados)

> anova(model4,model6) Model df AIC BIC logLik Test L.Ratio p-value model4 1 57 1354.742 1532.808 -620.3709 model6 2 59 1183.278 1367.592 -532.6389 1 vs 2 175.464 <.0001

E sim, vale a pena, é um modelo melhor, considerar a correlação espacial compensa o aumenta na complexidade do modelo necessária.

Veja que tiramos a dependência espacial dos dados, para corretamente comparar a produtividade das variedades.

05

E se a gente olhar os resíduos, da pra ver que o modelo não apresenta, a principio, nenhum problema (o resíduos tem uma distribuição aparentemente normal).

06

07

E uma curiosidade que vi lendo essa parte, é como é mais fácil que eu pensava, fazer comparações específicas, teste de contrastes nas variedades

> levels(spatialdata$variety) [1] "ARAPAHOE" "BRULE" "BUCKSKIN" "CENTURA" "CENTURK78" [6] "CHEYENNE" "CODY" "COLT" "GAGE" "HOMESTEAD" [11] "KS831374" "LANCER" "LANCOTA" "NE83404" "NE83406" [16] "NE83407" "NE83432" "NE83498" "NE83T12" "NE84557" [21] "NE85556" "NE85623" "NE86482" "NE86501" "NE86503" [26] "NE86507" "NE86509" "NE86527" "NE86582" "NE86606" [31] "NE86607" "NE86T666" "NE87403" "NE87408" "NE87409" [36] "NE87446" "NE87451" "NE87457" "NE87463" "NE87499" [41] "NE87512" "NE87513" "NE87522" "NE87612" "NE87613" [46] "NE87615" "NE87619" "NE87627" "NORKAN" "REDLAND" [51] "ROUGHRIDER" "SCOUT66" "SIOUXLAND" "TAM107" "TAM200" [56] "VONA"

Para comparar a primeira e a terceira variedade, basta a gente mandar um argumento a mais na anova.

> anova(model6,L=c(-1,0,1)) Denom. DF: 168 F-test for linear combination(s) (Intercept) varietyBUCKSKIN -1 1 numDF F-value p-value 1 1 7.565719 0.0066

Bem é isso ai, o script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo, ou mande um e-mail e até mais.

Referência:
Michael J. 2013 Crawley – The R Book, 2nd Edition Wiley 1076 pp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
library(nlme)
 
#Abrindo os dados do site do Crawley
spatialdata <- read.table("http://www.bio.ic.ac.uk/research/mjcraw/therbook/data/spatialdata.txt",header=T)
str(spatialdata)
 
#Olhando os dados
table(spatialdata$variety)
#Figura 1
boxplot(yield~variety,data=spatialdata,frame=F)
 
#Anova das variedades
model1 <- aov(yield~variety,data=spatialdata)
summary(model1)
 
#Existe uma diferença entre blocos, então adicionamos eles na analise
tapply(spatialdata$yield,spatialdata$Block,mean)
 
Block <- factor(spatialdata$Block)
model2 <- aov(yield~Block+variety+Error(Block),data=spatialdata)
summary(model2)
 
#Figura 2 - Efeito do espaço
par(mfrow=c(2,1))
plot(spatialdata$latitude,spatialdata$yield,pch=19,xlab="Latitude",ylab="Produtividade")
smooter <- loess(yield~latitude,data=spatialdata,span=0.75)
lines(predict(smooter), col='red', lwd=2)
 
plot(spatialdata$longitude,spatialdata$yield,pch=19,xlab="Longitude",ylab="Produtividade")
smooter <- loess(yield~longitude,data=spatialdata,span=0.75)
lines(predict(smooter), col='red', lwd=2)
 
#Levando em conta a latitude e longitude
model3 <- aov(yield~Block+variety+latitude+longitude,data=spatialdata)
summary(model3)
 
#Figura 3 - Explorando melhor como é o experimento
paleta<-palette(rainbow(length(levels(spatialdata$variety))))
plot(spatialdata$latitude,spatialdata$longitude,col=paleta[spatialdata$variety],pch=19,
     cex=2*(spatialdata$yield/max(spatialdata$yield)),frame=F,xlab="Latitude",ylab="Longitude",xlim=c(4,50),ylim=c(0,30))
text(spatialdata$latitude+1.5,spatialdata$longitude,spatialdata$Block,cex=0.8)
 
#Modelo inicial gls
model4 <- gls(yield~variety,spatialdata)
anova(model4)
 
#Figura 4 - Variograma
plot(Variogram(model4,form=~latitude+longitude))
dev.off()
 
#Adicionando duas estruturas de correção espacial
model5 <- update(model4,corr=corSpher(c(28,0.2),form=~latitude+longitude,nugget=T))
anova(model5)
 
model6 <- update(model4,corr=corRatio(c(12.5,0.2),form=~latitude+longitude,nugget=T))
anova(model6)
 
#Comparando modelos
anova(model5,model6)
anova(model4,model6)
 
#Obsevando o resultado final
anova(model6)
 
#Figura 5 - Variograma do modelo final
plot(Variogram(model6,resType="n"))
 
#Figura 6 - Residuos
plot(model6,resid( ., type="n")~fitted(.),abline=0)
dev.off()
 
#Figura 7 - Distribuição dos residuos
hist(resid(model6))
 
#Analise de constrastes
levels(spatialdata$variety)
anova(model6,L=c(-1,0,1))

As vezes, menos é melhor.

A gente ja viu vários posts com vários tipos de modelos, de poisson, binomial, mixtos e até expansões como overdispersion e zeros inflados.
Mas antes de mais nada, antes de sair tentando tudo e ver o que da mais certo. A gente tem que pensar, será que os parâmetros que estamos estimando fazem sentido. A gente consegue discutir aquele monte de número? Eles tem algum sentido biológico (para quem é da biologia)?

Vamos olhar um conjunto de dados, contagens de uma população ao longo do tempo, pensa que vemos o seguinte.

01

Não precisamos mentir. A maioria das pessoas (inclua eu nessa) pensaria em tacar uma regressão linear nesse conjunto de dados, um pensamento instantâneo, como se estivéssemos condicionados a fazer isso. Não que isso seja totalmente errado. Vamos ver como esse modelo ficaria.

02

Olha, não ficaria tão ruim, até que o modelo ajustaria. Mas e ai, o que estaríamos fazendo? Vendo como a população dessa espécie muda ao longo do tempo? E a inclinação é a intensidade como a população está mudando e o intercepto é como ela começou?

Eu acho que não é um interpretação tão ruim assim. Mas podemos melhorar mais e mais esse modelo, usando polinômios por exemplo.
Polinômios de segundo grau são uma parábola, de terceiro grau fica com uma curvinha, e cada vez o modelo vai ficando melhor, até que com 6 graus temos o seguinte.

03

Agora ta perfeito. Mas infelizmente, perfeito somente para esse grupo de dados, conforme a gente aumenta o número de parâmetros, a gente sempre melhora o ajuste, diminui o espaço para erro de medida, o problema é que o que os parâmetros fazem? O que significam? Nada mais, a gente so tem um ajuste bom para um modelo ininterpretável.
Olhe o que criamos, ele pode parecer bonito no grafico, mas não significa muito. E pense que fazemos modelos para pensar nos dados que coletamos, mas nos dados também que não coletamos, que outras pessoas vão coletar, e mais que isso, para entender como o mundo funciona, e será que temos muita capacidade de predição com esse modelo?
Se parar para pensar por um momento, a gente lembra que existem modelos de crescimento populacional que vão além de retas e linhas. Que cada letrinha, cada parametro tem um sentido especial que discutimos um bocado por aqui. Então se a gente passa tanto tempo lendo sobre eles, deveríamos gastar eles. até porque dai fica fácil discutir os dados.

04

Olha ae, talvez ele não seja tão bonito quanto ao ajuste, como o modelo polinimial, talvez até pior que o ajuste da regressão linear, mas é muito mais o que pensamos, e temos que pensar mais para frente, no tempo a frente.

05

Veja o que preveríamos para o futuro de acordo com cada modelo. O modelo de regressão linear preveria um crescimento constante, e o pior de tudo, o modelo polinomial preveria uma queda brusca na população. Agora você acha que tem algum evento obvio que faria a população diminuir. No fim das contas o modelo de crescimento populacional é o que faz a previsão mais razoável, que a população vai continuar crescendo. E é o que a gente acha. Os modelos de regressão são muito bons localmente, ali aonde temos dados, mas eles são péssimos fora da região onde estão os dados. Isso é algo que é bem comum, porque a regressão não assume muita coisa sobre o sistema biológico em si, o preço que se paga é isso.
Claro que se a gente não sabe nada, seria uma informação importante saber o que está acontecendo nesse período de tempo, e como está acontecendo, mas enquanto estamos construindo um conhecimento “bottom-up”, não sabemos nada, e estamos começando. Agora se a gente já tem modelos, a gente pode por eles a prova, algo mais “top-down”. A gente parte do que já sabe e vê se funciona na prática.

Agora quando falamos de chi-quadrado aqui, falamos também de como usar ele para seleção de modelos, mas aqui ferrou, porque para a seleção de modelos, precisamos de modelos aninhados. Por exemplo, modelos de Evolução de DNA temos vários modelos, mas um é uma versão mais simplificada de outro, um caso especial. A mesma coisa temos com os modelos de regressão linear, polinomial ou não. Agora o modelo de crescimento polinomial não aninha com os modelos de regressão, então acho que temos que escolher na argumentação, e não vai ter um valor de p dando sopa para justificar sem discussão essa escolha. Mas é a hora que pensar, apesar de no começo da curva de crescimento populacional (a vermelha) ser ruim no começo, lembramos que no começo a população é pequena, e coisas ao acaso acontecem mais comumente, so lembrar de genetic drift. Depois a população deve estabilizar mais e a curva fica melhor. Muito melhor porque claramente a população ta começando a crescer, e vai crescer bastante, e o modelo de crescimento populacional é o único que está dizendo isso.

Bem é isso ai, uma semana sem posts por conta de provas, mas agora voltando a atividades :), os scripts sempre no repositório recologia e é isso, até o próximo post.

tempo<-seq(0,12,by=2)
N<-c(230,1205,1150,2091,3010,6106,7323)
 
#figura 1
plot(N~tempo,xlab="Tempo",ylab="Tamanho populacional",frame=F,pch=19)
 
#modelo de regressão linear
modelo1<-lm(N~tempo)
modelo1
 
#figura 2
plot(N~tempo,xlab="Tempo",ylab="Tamanho populacional",frame=F,ylim=c(0,8000),pch=19)
curve(coef(modelo1)[1]+coef(modelo1)[2]*x,0,12,add=T,lwd=2,lty=1,col="green")
segments(x0=tempo, y0=N, x1 =tempo , y1 = predict(modelo1))
 
#modelo de regressão polinomial
modelo2<-lm(N~poly(tempo, degree=6, raw=TRUE))
modelo2
 
#figura 3
plot(N~tempo,xlab="Tempo",ylab="Tamanho populacional",frame=F,ylim=c(0,8000),pch=19)
curve(coef(modelo2)[1]+coef(modelo2)[2]*x+coef(modelo2)[3]*x^2+coef(modelo2)[4]*x^3+coef(modelo2)[5]*x^4+
      coef(modelo2)[6]*x^5+coef(modelo2)[7]*x^6,0,12,add=T,lwd=2,lty=1,col="blue")
segments(x0=tempo, y0=N, x1 =tempo , y1 = predict(modelo2))
 
#modelo não linear
modelo3<-nls(N~N0+exp(r*tempo),start=list(N0=1,r=1))
modelo3
 
#figura 4
plot(N~tempo,xlab="Tempo",ylab="Tamanho populacional",frame=F,ylim=c(0,8000),pch=19)
curve(coef(modelo3)[1]+exp(coef(modelo3)[2]*x),0,12,add=T,lwd=2,col="red")
segments(x0=tempo, y0=N, x1 =tempo , y1 = predict(modelo3))
 
#figura 5
plot(N~tempo,xlab="Tempo",ylab="Tamanho populacional",frame=F,ylim=c(0,10000),xlim=c(0,20))
curve(coef(modelo1)[1]+coef(modelo1)[2]*x,0,20,add=T,lwd=2,lty=1,col="green")
curve(coef(modelo2)[1]+coef(modelo2)[2]*x+coef(modelo2)[3]*x^2+coef(modelo2)[4]*x^3+coef(modelo2)[5]*x^4+
      coef(modelo2)[6]*x^5+coef(modelo2)[7]*x^6,0,20,add=T,lwd=2,lty=1,col="blue")
curve(coef(modelo3)[1]+exp(coef(modelo3)[2]*x),0,20,add=T,lwd=2,col="red")
points(N~tempo,pch=19,cex=1.2)
legend("topleft",col=c("green","blue","red"),lwd=2,bty="n",title ="Modelo",
       legend=c("Regressão Linear","Polinomial de grau 6","Crescimento populacional"))
 
logLik(modelo1)
logLik(modelo2)
logLik(modelo3)

A distribuição Chi-quadrado

Ola, vamos dar uma olhada aqui na distribuição chi-quadrado e tentar entender melhor porque da para usar ela para selecionar modelos.

A gente já falou da distribuição f aqui e da distribuição t aqui, e agora adicionando a distribuição chi-quadrado ao repertório.

A distribuição chi-quadrado sempre tem um papel importante no teste de hipóteses sobre frequências.

Basicamente, se a gente tem um conjunto de n elementos \{Z_1,Z_2,\cdots,Z_n\} independentes e normalizados, a soma dos quadrados deles vai ter uma distribuição aleatória chi-quadrado com n graus de liberdades.

Então:

\chi^2=\sum\limits_{i=1}^n\{Z_i^2\}

Mas o que quer dizer normalizado? No post anterior, aqui, a gente realizou essa operação, é só diminuir o valor pela média e dividir pelo desvio.

Na distribuição normal a gente faz:

{x-\alpha} \over \sigma

Agora no teste de qui-quadrado a gente faz o famoso:

{Observado\ -\ Esperado} \over {Esperado}

Certo, agora vamos contextualizar isso.
Na graduação, muita gente deve ter participado em algum momento de um trabalho, pelo menos eu participei, de reproduzir moscas Drosophila melanogaster. Famosas nos trabalhinhos de genética.

Chega uma hora que a gente deixa elas cruzarem e depois ia ver se havia realmente uma segregação independente dos cromossomos. No caso do cromossomos ligados ao sexo, temos fêmeas XX e machos XY.

Ai lembrando as leis do nosso colega Gregor Mendel , podemos fazer a quadrado de punnet.

Fêmea X X M X XX XX a c Y XY XY h O

Ou seja esperamos 50% de machos e 50% de fêmeas, uma razão de 1:1.

Então se nasceu 100 moscas, dessas 50 machos e 50 fêmeas, isso é perfeitamente o que esperamos.
Mas se nascer 49 macho e 51 fêmeas, não é exatamente o esperado, mas esta bem próximo, agora se a gente ver 5 machos e 95 fêmeas, isso é bem longe do que esperamos. Mas como transformar esse próximo ou longe em um número?
Bang, chegamos ao chi-quadrado. Vamos fazer uma função de chi-quadrado no R para o nosso caso e fazer alguns teste.

chi2_moscas<-function(observado) {
    if(is.vector(observado) && length(observado)==2) {
        esperado<-rep(sum(observado)/2,2)
        estatistica<-((observado[1]-esperado[1])^2/esperado[1])+((observado[2]-esperado[2])^2/esperado[2])
        return(estatistica)
    } else {
        print("Entrada errada. Os dados devem ser entrados como vetor numericos de tamanho 2!")
    }
}

então essa função recebe como entrada um vetor de dois números, a nossa contagem de machos e fêmeas, e vai calcular a estatística chi-quadrado, dado que nossa hipótese para as mosquinhas realmente seja uma razão sexual de 1:1, como o senhor Mendel propôs, não estamos tirando da cartola esse valor de 50% para cada sexo.

Veja que na função, a gente adicionou um teste para ver que a gente está entrando com os dados corretamente, caso alguém mande os dados errados, é comum a gente ver verificações em funções, fizemos uma aqui.
Usando nossa função a gente vê.

> chi2_moscas(c(200)) [1] "Entrada errada. Os dados devem ser entrados como vetor numericos de tamanho 2!"

Opa, entramos com um dado que só tem uma contagem, isso não faz o menor sentido, então não calculamos nada, a gente só faz conta quando faz sentido.

Agora vamos supor que pegamos nossa garrafinha la, e contamos 221 moscas macho e 475 moscas fêmeas, nossa isso ta bem fora da casinha, ta muito diferente de 1:1, se a gente calcula o chi-quadrado a gente vê…

> chi2_moscas(c(221,475)) [1] 92.6954

Um número grandão. Mas é grande mesmo? Vamos pensar que em outra garrafa, a gente pegou 332 machos e 375 fêmeas, isso é bem mais razoável segundo a hipótese não? Ta mais com cara de 1:1, e quando a gente calcula o valor de chi-quadrado a gente vê…

> chi2_moscas(c(332,375)) [1] 2.615276

Um valor baixinho. Legal, então um valor grandão quer dizer uma discrepância muito grande da nossa hipótese, e um valor baixo quer dizer que estamos diante de algo bem próximo da hipótese.
Isso não é difícil de enxergar com a formula que estamos fazendo.

Estamos calculando na função chi-quadrado o seguinte

{\chi^2}={\sum\limits_{i=1}^2} { {(Observado\ -\ Esperado)^2} \over {Esperado} }

Veja que vai até 2 a somatória, porque temos machos e fêmeas aqui nesse problema. O quadrado é para não ter problemas com números negativos, mas o que acontece é que se o observado é próximo do número esperado, vai dar um número pequenininho ali no numerador, se for exatamente o mesmo número, da zero, e zero dividido por qualquer coisa da zero, se for um número pequeno, vai ser algo pequeno dividido por um número que vai dar um número pequeno, agora conforme essa diferença aumenta, o numerador aumenta, e cada vez o número fica maior.

Continuando com esse número total de 100 moscas, para ficar fácil fazer contas, podemos calcular o chi-quadrado para os casos de 49-51, 48-52 … até 1-99 e ver o que acontece.

01

Olha ai, conforme a discrepância em relação a nossa hipótese aumenta, o valor de chi-quadrado aumenta também, exponencialmente.

Mas num é qualquer discrepância que a gente vai para e assumir que o senhor Mendel tava errado. A gente tem que lembrar como as coisas estão acontecendo.

Cada mosca tem 50% de chance de ser macho o fêmea, independente da contagem atual, pode já ter nascido 99 moscas fêmeas, a próxima, se Mendel está certo, tem 50% de ser macho. Nos, seres humanos, somo péssimos geradores de números aleatórias por isso, se você contar 99 fêmeas, você vai primeiro pensar que algo acontece, as vezes sem levar em conta que esse resultado é perfeitamente possível dentro da teoria de Mendel, ele é raro, mas perfeitamente possível.

Agora a gente pode fazer um simulação no R e ver o quão raro é um caso desse, alias, a gente pode ao invés de guardar quantas moscas de qual sexo nasceram, simular os nascimentos e guardar somente o valor de chi-quadrado para o resultado.

Então é so a gente ir simulando os nascimento usando a função sample

> sample(c("Macho","Femea"),100,replace=T) [1] "Femea" "Macho" "Femea" "Femea" "Macho" "Femea" "Macho" "Macho" "Macho" "Femea" "Femea" [12] "Femea" "Femea" "Femea" "Femea" "Femea" "Femea" "Femea" "Macho" "Macho" "Femea" "Femea" [23] "Femea" "Femea" "Macho" "Macho" "Femea" "Macho" "Macho" "Macho" "Macho" "Macho" "Macho" [34] "Macho" "Macho" "Macho" "Macho" "Macho" "Macho" "Femea" "Macho" "Macho" "Femea" "Macho" [45] "Femea" "Macho" "Macho" "Femea" "Femea" "Macho" "Femea" "Macho" "Macho" "Macho" "Femea" [56] "Macho" "Femea" "Femea" "Femea" "Macho" "Femea" "Femea" "Femea" "Macho" "Macho" "Macho" [67] "Macho" "Femea" "Macho" "Macho" "Macho" "Femea" "Macho" "Femea" "Macho" "Macho" "Femea" [78] "Femea" "Femea" "Macho" "Macho" "Femea" "Macho" "Femea" "Macho" "Macho" "Femea" "Macho" [89] "Macho" "Femea" "Femea" "Macho" "Macho" "Femea" "Macho" "Femea" "Femea" "Femea" "Macho" [100] "Macho"

E ir guardando o resultado num vetor, já que é somente um número.

for(i in 1:10000) {
    experimento<-sample(c("Macho","Femea"),100,replace=T)
    macho<-sum(experimento=="Macho")
    femea<-sum(experimento=="Femea")
    estatistica[i]<-chi2_moscas(c(macho,femea))
}

Agora para facilitar a, a gente olha o resultado num histograma.

02

Humm, então existem poucas chances de valores extremos, mas valores próximos de zero são bem comuns.
Bem quando alguém descreve uma distribuição, o cara faz uma probability density function ou PDF, que nada mais é que uma função que descreve esse resultado. O R ja tem pronta um monte desses PDF, normalmente é d grudado no nome da distribuição, aqui é o dchisq.

Vamos dar uma olhada, se é parecido com o que encontramos.

03

Olha ai que legal, fica bem parecido com o que encontramos, ou seja antes da gente fazer a simulação, ja tinha como saber como seria o resultado com essa distribuição. Eu marquei o quantile de 95% também. O que esse risco azul ta mostrando, que 95% dos casos estão a esquerda desse número, e 5% a direita.

É fantástico pensar que isso funciona, podemos comparar esse numero que calculamos com o quantile com os valores de chi-quadrado que calculamos e ver se ele funciona de verdadinha.

> sum(estatistica<qchisq(0.95,1))/length(estatistica)
[1] 0.9456

E num é que funciona, 0.9456 são menor que ele, isso é praticamente 95%, muito boa predição.

Agora se o nosso valor que encontramos na garrifinha for muito extremo, podemos começar a suspeitar da teoria do Mendel, não que isso implique que ela é errada, mas podemos suspeitar, isso é o tal do p ser significativo, é a gente começar a suspeitar que algo estranho ta acontecendo, esse algo é permitido acontecer, mas incomum, muito incomum.

Vamos pensar em dois casos aqui, vamos pensar que contamos 43 machos e 57 femeas, se a gente calcular o valor de qui-quadrado

> experimento1<-c(43,57) > chi2_moscas(experimento1) [1] 1.96

Da 1.96, isso é bem baixinho, se a gente compara com os resultados da nossa simulação

> sum(chi2_moscas(experimento1)<estatistica)/length(estatistica)
[1] 0.1277

tem pelo menos 13% dos valores tão baixos quanto 1.96, não tem muito porque levantar suspeitas sobre a teoria de Mendel, e podemos comparar o procedimento que estamos fazendo com o teste de chi-quadrado já pronto do R, usando a função chisq.test

> chisq.test(experimento1,p=c(0.5,0.5)) Chi-squared test for given probabilities data: experimento1 X-squared = 1.96, df = 1, p-value = 0.1615

E a gente tem o mesmo valor da estatística calculada, 1.96 (estamos fazendo a mesma continha), e um valor de p que traz a mesma conclusão, 0.16.

Agora vamos pensar que em outra experimento vemos um resultado bem diferente, 21 machos e 79 fêmeas.
Aqui começa a ter cara de um desvio na razão sexual nervoso. A gente pode calcular o valor de chi-quadrado e vemos o seguinte.

> experimento2<-c(21,79) > chi2_moscas(experimento2) [1] 33.64

Nossa, 33.64, se a gente comprar com os valores que obtemos ao acaso, fazendo a simulação onde machos e fêmeas tem 50% de chance cada de nascer.

> sum(chi2_moscas(experimento2)<estatistica)/length(estatistica)
[1] 0

Nenhum, absolutamente nenhuma das 10 mil simulações foi tão extrema assim. Ou seja, algo deve estar acontecendo, ou a teoria do Mendel esta errada aqui, ou alguém esta pregando uma peça no laboratório, não sabemos o que, varias coisas podem influenciar esse resultado, mas sabemos que se for para confiar no azar para ver um resultado assim, esse resultado é muitoooo raro.

Podemos fazer o teste de chi-quadrado com a função do R aqui também.

> chisq.test(experimento2,p=c(0.5,0.5)) Chi-squared test for given probabilities data: experimento2 X-squared = 33.64, df = 1, p-value = 6.631e-09

E olha ai, um valor de p de 6.631e-09, lembrando que isso está em notação cientifica, esse e-09 significa que esse número é…

> format(6.631e-09,scientific=F) [1] "0.000000006631"

Muito legal né, então podemos usar essa distribuição para compara se algo está próximo da hipótese ou longe. Agora outro coisa, além de razão sexual, que tem uma distribuição que segue uma distribuição chi-quadrado são valores de verosimilhança, ou likelihood, e AIC também (AIC é loglikelihood corrigido meu, claro que tem que ser igual o likelihood quanto a distribuição).

E onde entra o chi-quadrado aqui?

Imagine que temos duas variáveis, um preditor e uma resposta, mas que nao tem relação nenhuma, os dois são amostras ao acaso.

> preditor<-rnorm(30) > resposta<-rnorm(30)

A gente pode tentar modelar eles tanto como variáveis que tem uma relação entre si, como um modelo onde a resposta simplesmente vem de uma distribuição normal, ou seja, assumir que ela simplesmente vem de uma distribuição normal é mais simples que assumir que ela tem relação com algo.

> modelo1<-lm(resposta~preditor) > modelo2<-lm(resposta~1)

A gente pode avaliar esses dois modelos e seus parâmetros.

> summary(modelo1) Call: lm(formula = resposta ~ preditor) Residuals: Min 1Q Median 3Q Max -1.7239 -0.6622 -0.1236 0.5255 2.1173 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.06622 0.17136 0.386 0.702 preditor 0.00845 0.19107 0.044 0.965 Residual standard error: 0.9379 on 28 degrees of freedom Multiple R-squared: 6.984e-05, Adjusted R-squared: -0.03564 F-statistic: 0.001956 on 1 and 28 DF, p-value: 0.965 > summary(modelo2) Call: lm(formula = resposta ~ 1) Residuals: Min 1Q Median 3Q Max -1.7339 -0.6600 -0.1162 0.5272 2.1214 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.06592 0.16826 0.392 0.698 Residual standard error: 0.9216 on 29 degrees of freedom

Mas e ai, qual é melhor? Ue o que tiver o melhor likelihood, o menor custo, aquele que acumula menos resíduo. Mas essa diferença pode ser ao acaso, a diferença de likelihood, e podemos testar isso usando a distribuição chi-quadrado.

> anova(modelo1,modelo2,test="Chisq") Analysis of Variance Table Model 1: resposta ~ preditor Model 2: resposta ~ 1 Res.Df RSS Df Sum of Sq Pr(>Chi) 1 28 24.628 2 29 24.630 -1 -0.0017202 0.9647

Olha ai, comparamos esses modelos e vemos que eles são iguais, tem a mesma capacidade de predição. Então por parcimônia, é melhor eu escolher o mais simples uai, porque ficar com algo mais complexo, se a complexidade nao ajuda em nada, essa é a idea da navalha de occan, e vamos seguir ela, mas tendo um teste para sustentar essa escolha, pode diminuir a discussão entre nos e nossos pares, afinal se todo mundo concorda com esses métodos, vamos todos chegar as mesmas conclusões.

Ficamos por aqui, agora com a distribuição chi-quadrado melhor compreendida, espero eu, os scripts sempre no repositório recologia e é isso, até o próximo post.

set.seed(123)
chi2_moscas<-function(observado) {
 
    if(is.vector(observado) && length(observado)==2) {
 
        esperado<-rep(sum(observado)/2,2)
        estatistica<-((observado[1]-esperado[1])^2/esperado[1])+((observado[2]-esperado[2])^2/esperado[2])
        return(estatistica)
 
    } else {
        print("Entrada errada. Os dados devem ser entrados como vetor numericos de tamanho 2!")
    }
 
}
 
 
chi2_moscas(c(200))
chi2_moscas(c(221,475))
chi2_moscas(c(332,375))
 
 
dados<-vector()
 
for(i in 1:50) {
    dados[i]<-chi2_moscas(c(50-i,50+i))
}
 
#figura 1
plot(dados,frame=F,type="p",pch=19,cex=0.5,xaxt="n",ylab="Estatisthica Chi-quadrado",xlab="Razão de Machos para Femeas")
axis(1,at=c(1,13,25,38,49),label=c("1:1","1:13","1:25","1:38","1:49"))
 
 
sample(c("Macho","Femea"),100,replace=T)
 
estatistica<-vector()
 
for(i in 1:10000) {
    experimento<-sample(c("Macho","Femea"),100,replace=T)
    macho<-sum(experimento=="Macho")
    femea<-sum(experimento=="Femea")
    estatistica[i]<-chi2_moscas(c(macho,femea))
}
 
#figura 2
hist(estatistica,prob=T,main="Estatística Chi-Quadrado",xlab="Estatística")
 
 
 
#figura 3
hist(estatistica,prob=T,main="Estatística Chi-Quadrado",xlab="Estatística")
curve(dchisq(x,1),0,15,add=T,lwd=3,lty=2,col="red")
abline(v=qchisq(0.95,1),col="blue",lty=3)
legend("topright",lwd=c(3,1),col=c("red","blue"),lty=c(2,3),legend=c("Distribuição Chi-Quadrado","Quantile de 95%"),
       bty="n")
 
 
sum(estatistica<qchisq(0.95,1))/length(estatistica)
 
 
experimento1<-c(43,57)
chi2_moscas(experimento1)
sum(chi2_moscas(experimento1)<estatistica)/length(estatistica)
chisq.test(experimento1,p=c(0.5,0.5))
 
 
experimento2<-c(21,79)
chi2_moscas(experimento2)
sum(chi2_moscas(experimento2)<estatistica)/length(estatistica)
chisq.test(experimento2,p=c(0.5,0.5))
 
format(6.631e-09,scientific=F)
 
preditor<-rnorm(30)
resposta<-rnorm(30)
 
modelo1<-lm(resposta~preditor)
modelo2<-lm(resposta~1)
 
summary(modelo1)
summary(modelo2)
 
anova(modelo1,modelo2,test="Chisq")

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)