Integral pela regra de Simpson

A gente já viu aqui, como calcular uma integral pela regra do trapézio. Agora vamos tentar ver como funciona a regra de Simpson.

A regra de Simpson consiste dividir o intervalo que queremos calcular a integral, [a,b] em aproximadamente n subintervalos, onde n tem que ser par, então em cada par de conjunto de subintervalos consecutivos, aproximamos o comportamento de f(x) por uma parábola (polinômio de segundo grau) ao contrário de uma reta, como na regra do trapézio.

Seja u < v < w[/latex] qualquer três pontos distantes [latex]h[/latex] entre eles. Para [latex] x \in [u,w] [/latex] nos queremos aproximar [latex]f(x)[/latex] por uma parábola que passa pelos pontos [latex](u,f(u))[/latex], [latex](v,f(v))[/latex] e [latex](w,f(w))[/latex]. Existe exatamente uma parábola [latex]p(x)[/latex] que faz isso e é dada pela formula:  [latex]p(x) =f(u) \frac{(x-v)(x-w)}{(u-v)(u-w)} + f(v) \frac{(x-u)(x-w)}{(v-u)(v-w)} + f(w) \frac{(x-u)(x-v)}{(w-u)(w-v)}[/latex]  Pode parecer assustador, mas essa é a formula da parábola que vai passar por [latex](u, v, w)[/latex].  Se definirmos uma função     3f1c9cc348a3fb7b57d0b5898f87ca79000   E três pontos quaisquer  <pre>  <div style='font-family:Courier New;font-size:10pt;'>  u<-1  v<-2  w<-3  </div>  </pre>  Veja que essa equação define a parábola passando certinho em cima desses pontos, e tome cuidado com os parenteses, eu apanhei um pouquinho porque estava esquecendo alguns e a reta não passava em cima dos pontos heheh.   3f1c9cc348a3fb7b57d0b5898f87ca79001   <a href="http://recologia.com.br/wp-content/uploads/2015/01/044.jpg"><img src="http://recologia.com.br/wp-content/uploads/2015/01/044.jpg" alt="04" width="480" height="480" class="aligncenter size-full wp-image-2812" srcset="http://recologia.com.br/wp-content/uploads/2015/01/044.jpg 480w, http://recologia.com.br/wp-content/uploads/2015/01/044-150x150.jpg 150w, http://recologia.com.br/wp-content/uploads/2015/01/044-300x300.jpg 300w" sizes="(max-width: 480px) 100vw, 480px" /></a>   Então antes era trapézio, retos, e agora temos curvas. Como uma aproximação para área abaixo da curva [latex]y=f(x), nos usamos \int_{u}^w p(x)dx. Uma demonstração grande e lenta (que eu não sei fazer) chega ao seguinte:

\int_{u}^w p(x)dx = \frac{h}{3} (f(u)+4f(v) + f(w))

Agora, assumindo que n é par, nos somamos as aproximações para os subintervalos [x_{2i},x_{2i+2}] para obter a aproximação de Simpson S para a integral \int_{a}^b f(x)dx

Regra de Simpson:

S=\frac{h}{3}(f(x_0)+4f(x_1)+2f(x_2)+4f(x_3)+ \cdots + 4f(x_{n-1})+f(x_n))

Note que f(x_i) para um i par é sempre multiplicado por 4, enquanto o f(x_i) para um i impar (exceto 0 e n) é multiplicado por 2 já que eles aparecem em dois subintervalos.

A regra de Simpson da resultados exatos se f(x) é uma função quadrática já que ela é uma aproximação para cada parte de f(x) por uma parábola, Surpreendentemente, ela também da o valor exato se f(x) é uma função cúbica. Mas de forma geral ela tem melhores resultado que a regra do trapézio.

Novamente, fazer isso funcionar é so usar a regra diretamente no R

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
simpson_n <- function(ftn, a, b, n = 100) {
    #Integração númerica de ftn de a até b
    #Usando a regra de simpson para n subdivisões
    #
    #ftn é uma função de uma unica variavel
    #assumimos que a < b e n é um número par positivo
    n <- max(c(2*(n %/% 2), 4))
    h <- (b-a)/n
    x.vec1 <- seq(a+h, b-h, by = 2*h)
    x.vec2 <- seq(a+2*h, b-2*h, by = 2*h)
    f.vec1 <- sapply(x.vec1, ftn)
    f.vec2 <- sapply(x.vec2, ftn)
    S <- h/3*(ftn(a) + ftn(b) + 4*sum(f.vec1) + 2*sum(f.vec2))
    return(S)
}

Veja então que primeiro a gente faz a divisão do n por 2, %/% é a parte inteira da divisão no R, e multiplica por 2, para saber quantos pontos usaremos, veja que isso está dentro de um vetor com 4, já que se você escolher um n de 1, 2 ou 3 não tem como, tem que ser pelo menos 4, temos que ter 3 pontos pelo menos e tem que ser par, e lembre-se que tudo que é multiplicado por 2 vira par.

Depois calculamos o h, que é o incremento do intervalo de a a b para dar n partições. Lembrando que na formula, só tem peso os pontos no meio, por isso a gente começa do a+h e termina no b-h o vetor, que é onde vai ter peso, o a e b estão separados ali em baixo.

Mas ai calculamos f(x) para esses valores e jogamos na formula da regra de Simpson.

Então para calcular a seguinte integral

\int_{0}^1 4x^3 dx

No R:

1
2
3
ftn <- function(x) {
    return(4 * x^3)
}

01

> simpson_n(ftn, 0, 1, 20) [1] 1

Que podemos confirmar com a função integrade do R que é esse valor mesmo

> integrate(ftn,0,1) 1 with absolute error < 1.1e-14

E na verdade, essa função da para integrar na mão.

\int 4x^3 dx = x^4 + constante

como

\int_{a}^b f(x) dx = F(b)-F(a) onde F é a antiderivada de f temos

\int_{0}^1 4x^3 dx = 1^4-0^4 = 1

Muito legal isso, pena que nem toda função é simples assim, não que isso importe para todo mundo. Carl Friedrich Gauss, o minininho matemágico, que somou de 1 a 100 mais rápido que todo mundo, tem entre os muitos atos prodigiosos de sua vida ter compilado a tabelinha de \phi(z)=\int_{-\inf}^z \frac{1}{\sqrt{2\pi}}e^\frac{-x^2}{2}dx a mão com muitas casas decimais de precisão. Mas com o computador, calcular isso fica fácil, inclusive no R, a função \phi(z) ja vem implementada na função pnorm.

Podemos definir a função do Gauss, ou distribuição Gaussiana, ou distribuição normal, que tem um monte de nomes para nos confundir, bem podemos definir ela no R

1
phi <- function(x) return(exp(-x^2/2)/sqrt(2*pi))

Bem como sua integral.

1
2
3
4
5
6
7
Phi <- function(z) {
    if (z < 0) {
        return(0.5 - simpson_n(phi, z, 0))
    } else {
        return(0.5 + simpson_n(phi, 0, z))
    }
}

E usar

> phi(1.96) [1] 0.05844094 > Phi(1.96) [1] 0.9750021

Veja que essas são as funções de densidade e de probabilidade da distribuição normal

> dnorm(1.96) [1] 0.05844094 > pnorm(1.96) [1] 0.9750021

E o número magico 1.96 é o famoso intervalo de confiança de 95%

> Phi(1.96)-Phi(-1.96) [1] 0.9500042 > pnorm(1.96)-pnorm(-1.96) [1] 0.9500042

02

Resta ainda pensar na acurácia da regra de Simpson. Podemos estimamar \int_{0.01}^1 f(\frac{1}{x})dx = -log(0.01) para uma sequência crescente de valores de n, o número de partições.

03

No gráfico de log(erro) contra o log(n) da pra ver uma inclinação negativa de aproximadamente -4 para valores de n, indicando que o erro decai como n^{-4}. Isso pode de fato ser demostrado como comum, de forma geral, para funções f com uma quarta derivativa contínua.

Bem é isso ai, hora que tomar coragem de começar a postar sobre os modelos em evolução, essa preocupação com cálculo vai começar a ser recompensada, você verá.
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:
Owen Jones, Robert Maillardet & Andrew Robinson 2009 – Introduction to Scientific Programming and Simulation Using R. CRC Press 499pp

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
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
#Post
#http://recologia.com.br/2015/01/integral-pela-regra-de-simpson
 
ftn <- function(x) {
    return(4 * x^3)
}
 
p <- function (x,u,v,w) {
    saida<- ftn(u) * ( ((x-v)*(x-w))/((u-v)*(u-w)) ) +
            ftn(v) * ( ((x-u)*(x-w))/((v-u)*(v-w)) ) +
            ftn(w) * ( ((x-u)*(x-v))/((w-u)*(w-v)) )
    return(saida)
}
 
u<-1
v<-2
w<-3
 
#figura 1
pontos<-c(u,v,w)
plot(pontos,ftn(pontos),xlim=c(-2,6),ylim=c(0,120),frame=F,pch=19)
 
vetor<-seq(0,4,by=0.1)
points(vetor,p(x=vetor,u=1,v=2,w=3),type="l",lty=3)
 
 
#Itegral pela Regra de Simpson
 
simpson_n <- function(ftn, a, b, n = 100) {
    #Integração númerica de ftn de a até b
    #Usando a regra de simpson para n subdivisões
    #
    #ftn é uma função de uma unica variavel
    #assumimos que a < b e n é um número par positivo
    n <- max(c(2*(n %/% 2), 4))
    h <- (b-a)/n
    x.vec1 <- seq(a+h, b-h, by = 2*h)
    x.vec2 <- seq(a+2*h, b-2*h, by = 2*h)
    f.vec1 <- sapply(x.vec1, ftn)
    f.vec2 <- sapply(x.vec2, ftn)
    S <- h/3*(ftn(a) + ftn(b) + 4*sum(f.vec1) + 2*sum(f.vec2))
    return(S)
}
 
ftn <- function(x) {
    return(4 * x^3)
}
 
#Figura 2
vetor<-seq(-1,2,by=0.01)
plot(vetor,ftn(vetor),type="l")
vetor<-seq(0,1,by=0.01)
polygon(x=c((vetor),rev(vetor)),y=c(ftn(vetor),rep(0,length(vetor))),col="red")
abline(h=0,lty=3)
 
 
simpson_n(ftn, 0, 1, 20)
integrate(ftn,0,1)
 
#Distribuição normal
 
#função densidade
phi <- function(x) return(exp(-x^2/2)/sqrt(2*pi))
 
#Função cumulativa (probabilidade)
Phi <- function(z) {
    if (z < 0) {
        return(0.5 - simpson_n(phi, z, 0))
    } else {
        return(0.5 + simpson_n(phi, 0, z))
    }
}
 
phi(1.96)
Phi(1.96)
 
dnorm(1.96)
pnorm(1.96)
 
Phi(1.96)-Phi(-1.96)
pnorm(1.96)-pnorm(-1.96)
 
z <- seq(-5, 5, by = 0.1)
phi.z <- sapply(z, phi)
Phi.z <- sapply(z, Phi)
jpeg("02.jpg")
plot(z, Phi.z, type = "l", ylab = "", main = "phi(z) and Phi(z)")
lines(z, phi.z)
dev.off()
 
 
ftn <- function(x) {
    return(1/x)
}
 
S <- function(n) {
    simpson_n(ftn, 0.01, 1, n)
}
 
n.vec <- seq(10, 1000, by = 10)
S.vec <- sapply(n.vec, S)
 
#Figura 4
par(mfrow = c(1, 2), pty="s", mar=c(4,4,2,1), las=1)
plot(n.vec, S.vec + log(0.01), type = "l",
xlab = "n", ylab = "erro")
plot(log(n.vec), log(S.vec + log(0.01)), type = "l",
xlab = "log(n)", ylab = "log(erro)")

Séries de cursos de R no Coursera e no EdX

Ola.
Esse post é só para compartilhar duas séries de cursos que acho bem interessante, para aprender mais sobre análise de dados e R.

Primeiro, no Coursera, a John Hopkins University está oferecendo uma série de cursos, na forma de uma especialização em Data Science.

São nove cursos no total, que tem duração de uma semana cada um. Os tópico de cada curso são bem específicos, mas eles se complementam muito bem, principalmente porque eles não são apenas sobre análise de dados, mas englobam todas as partes do processo de análise de dados.

Por exemplo, a gente tem o curso “Getting data and Cleaning data“. Esse curso achei muito legal, porque ele mostra varias formas de ler dados no R, das mais variadas fontes. Na maioria dos cursos que a gente faz, que eu fiz na faculdade pelo menos, a gente acaba ficando limitado a ler dados de arquivos do tipo delimitados por tabulações, e arquivos delimitados por vírgulas, e raramente a gente ainda entra no porque existe um read.csv e um read.csv2 no R. E essas coisas a gente vê nesse curso, ainda como abrir dados diretos do excel, que é talvez o lugar mais comum aonde os dados estão guardados na biologia em geral, mas também como ler dados de banco de dados sql, arquivos xml e páginas html, além de como ajeitar os dados, com o plyr por exemplo, e ler arquivos muito grandes, onde a read.table da vida pode falhar, já que ele tem que ler os dados e dar um parse neles, então ele ocupa bastante memoria.

Outro curso muito legal dessa série é o “Exploratory Data Analysis“, que ensina a fazer gráficos bonitos, além de esclarecer bastante coisas sobre como salvar e exportar os gráficos de várias formas, e usar o gráficos do tipo lattice, que quando os dados começam a ter muitas variáveis, começa a ficar mais difícil de representar graficamente, mas o lattice é uma mão na roda. Além de ver o ggplo2, que faz gráficos muito bonitos também, tem uma série de ferramentas legais, para fazer várias coisas bacanas, e que no final a gente acaba precisando para pelo menos conseguir ler código de outras pessoas, por mais que você prefira ou não usar o sistema base para gráficos. E gráficos são essenciais, quando alguém quer me mostrar alguma coisa, ou pede ajuda, eu sempre tento fazer alguma figura primeiro, para começar a entender o problema, nós, como seres visuais, ganhamos muito fazendo gráficos, e toda habilidade para fazer eles melhores é valida.

Não vou falar de todos os cursos aqui, mas a forma de avaliação também é bem interessante, porque além dos quizzes, perguntas de múltipla escolha para avaliação, todos praticamente envolvem alguma atividade prática, que você tem que fazer algum tipo de análise e escrever um relatório, que normalmente é usando R Markdown (um exemplo meu aqui, o relatório final do curso Practical Machine Learning, que tem que ser entregue num repositório do github), ou programar alguma coisa, como gráficos dinâmicos usando Shinny (minha primeira figura do shinny nesse repositório).
Isso é muito legal porque você recebe um feedback de pessoas, vê o trabalho de outros alunos para comparar, que vária de coisas ruins e erradas até trabalhos excepcionais, que faz você se perguntar porque a pessoa está fazendo o curso se está naquele nível tão avançado de conhecimento.

Os cursos recomeçam todo mês praticamente, então da pra fazer de um em um, ou mais de um por vez, e se não conseguir acompanhar, é só começar denovo, tendo a possibilidade de pagar para receber um certificado verificado, o que não deve ser muito útil no Brasil, mas eu pagaria somente pela qualidade do material (mas minha situação financeira não anda boa hehe).

Antes desses cursos, eu fiz um monte de cursos do pessoal da John Hopkins, e eles tem melhorado cada vez mais na qualidade dos cursos, acho que de tanto oferecer disciplinas no coursera, eles estão cada vez fazendo cursos mais legais, e num formato legal para atender o público em geral.

Bem, desde o ano passado eu comecei a tentar alguns cursos no EdX também, que é uma plataforma de ensino massivo online como o Coursera, mas que começou com o MIT e Harvard, talvez para competir com a popularidade do Coursera, que foi fundado por dois professores de Standford, Andrew Ng and Daphne Koller.

Bem o EdX tem muitos cursos que envolvem matemática, por causa do MIT provavelmente, que é bem legal para quem curte, mas ele abriu uma série de curso sobre R também, que são do professor Rafael Irizarry, uma coisa legal é que esse cara era professor la no mesmo departamento do povo que da os cursos da John Hopkins no Coursera, mas ele mudou para Havard, ai foi parar no EdX.

Mas ele vai dar uma série de cursos no EdX, todos usando R para análise de dados, começando com o Statistics and R for the Life Sciences, e na descrição desse curso, você pode ver os oito cursos que fazem parte dessa série. São mais coisas de análise para bioinformática, mas pode ser um boa oportunidade de aprender mais sobre o biocondutor, que é um projeto bem grande junto com o R.
Eu estou no primeiro curso, mas está parecendo bem interessante, é bem simples até, mas vale a pena seguir, até para seguir a série e ver a linha de raciocínio dele, e praticar.

Bem já que estamos falando de EdX, outros dois cursos que fiz la, que achei muito legais são o de Introdução ao Linux, muito legal para aprender as coisas que todo mundo acha que você sabe e você não tem ideia, pelo menos nos livros, as pessoas assumente sempre um conhecimento de linux que eu pelo menos nunca aprendi formalmente. Mas esse curso é bem básico, simples, mas ao mesmo tempo legal, ensinando pequenas coisas como pipes e sistema de arquivos do linux de uma forma que da para entender bem, além do que o curso é oferecido pela Linux Foundation. Outro curso é o Quantitative Biology Workshop que é uma salada mista de tópicos, mas bem legal para ver a matemática aplicada a biologia, eu já fiz até um tópico sobre um modelinho que vi neste curso aqui, mas ele, como um monte de curso do EdX, é feito usando Matlab, mas com um pouco de paciência, da para encarar o Matlab, e ir traduzindo o código para o R para conseguir fazer outras coisas, ou mesmo só porque é legal.

Bem as vezes são raros os cursos sobre análise de dados ou R que podemos participar, seja por tempo para estar presente em aulas, seja por falta de grana para viajar para onde esses cursos são ofertados, seja por disponibilidade de professores ou cursos mesmo nas instituições, então esses cursos online do tipo MOOC são uma grande oportunidade, tenho feito cursos assim a algum tempo e eu acho bem legal, o acesso que é proporcionado por eles, por exemplo, você ter aulas de caras muito bons na área, com uma ótima didática, já que as universidades não desviar um professor ruim para montar um curso desses, até porque eles são também uma propaganda da instituição, e o fato de que você pode falhar nos cursos sem grandes problemas, enquanto se você falha num curso na pós-graduação por exemplo, você pode perder até bolsa de estudos, dependendo das normas dos cursos, num MOOC ninguém nem precisa ficar sabendo. Mas na especialização em data science do coursera por exemplo, eu desiste de vários cursos várias vezes, antes de terminar, talvez porque as vezes eu era ganancioso e queria fazer muitos cursos de uma vez, e não dava conta de acompanhar tudo, ou porque perdia prazos de entregar os exercícios, mesmo o relatório mais simples, não da para fazer em uma hora, eu não consegui! E sempre da para aprender muito nas discussões nos fóruns de mensagens desses cursos.

Mas é isso ai, só queria comentar as possibilidades de cursos legais sobre R que tem por ai, pelo menos que eu conheço. Além do que eu acho cursos assim mais dinâmicos de ler livros de cara, a informação vem mais mastigadinha. Assim quando a gente vai ler um livro sobre o assunto, ou um artigo científico, tudo começa a ficar mais fácil de digerir, com a boa base que esses cursos vão dando, mas é claro que não podemos ficar somente neles. Mas eu sempre tento acompanhar pelo menos um curso, é sempre legal.

Até o próximo post, e qualquer sugestão de curso, ou oportunidade legal de aprendizado, deixa um comentário ou mande um e-mail, que eu quero saber 🙂

Crescimento populacional com atraso para o efeito de dependência da densidade

Voltando ao crescimento populacional, tópico que tem alguns posts aqui já, so dar uma olhada nessa busca aqui do blog.

Mas quando a gente olha o crescimento populacional em modelos como o do crescimento logístico.

  N_{t+1} = r\cdot N_{t} \cdot (1- \alpha * N_{t} )

Lembrando que aqui, o alpha (\alpha) é o efeito de um indíviduo sobre a população, que é a mesma coisa que um sobre a capacidade suporte, ( k= {1\over \alpha} ).

Bem podemos escrever essa formula como normalmente é vista nos livros, com capacidade suporte ao invés de efeito da população sobre ela mesma.

  {dN\over dt} = r\cdot N\cdot (1-{N\over K})

Uma coisa legal que eu ouvi, é que a capacidade suporte não precisa ser pensada somente no efeito da própria população nela mesma, como falta de espaço de uma população de alga para crescer numa placa de petri, mas podemos pensar nela também como um proxy de outras condições que mudam no ambiente, como a “deterioração do ambiente”, que tem seu valor próprio, mas pense que esse valor tem uma relação direta com o tamanho populacional, ele ferra a população, logo ao invés de medir ele, podemos usar o próprio tamanho populacional, que representa esse valor, e simplifica em muito a modelagem, diminuindo uma variável por exemplo.

Mas ok, o ponto que queremos chegar é que aqui

 (1-{N\over K})

é a hora que acontece a diminuição do crescimento de acordo com o tamanho da população, mas isso é baseado no tamanho populacional aqui e agora.
Mas e se, pensando numa população de ratinhos na mata, quando o ratinho tem sua ninhadinha, logo que ele teve os filhos, não temos ainda o efeito da deterioração do ambiente por exemplo, esses filhotes tem que crescer, e viver e se reproduzir para deteriorar o ambiente, mas na vida deles eles não sentem essa deterioração, quem sente isso são seus descendentes. Quem vai lidar com os problemas causados são os filhos, netos e bisnetos. Quando agora os ratinhos estão la desfrutando da fartura e comendo todas as sementes, quem vai se ferrar são as próximas gerações, que se todas as sementes são comidas, não serão germinadas novas arvores e arbustos que deveriam fornecer comida para outras gerações.
Agora como será que isso impacta a dinâmica populacional?

Bem vamos fazer uma simulação aqui para tentar ver o efeito de dependência da densidade atrasado (Delayed density dependence).

Então a gente precisa dos parâmetros populacionais que usamos no modelo logístico:

1
2
3
4
5
N0 <- 1 #População inicial
r  <- 0.1 #Taxa de crescimento Populacional
K  <- 100 #Capacidade Suporte
t  <- 200 #Número de gerações
tau<- 20  #Atraso de resposta para o efeito dependente de densidade

E mais um parâmetro especial, que é o tau, ou o lag, o quanto vai demorar para a população sentir o efeito do tamanho populacional de agora.

E vamos implementar uma simulação discreta da seguinte forma:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
for(i in 1:length(generation)) {
        if(i == 1){
            population[i] <- N0
            dNdt[i] <- r*population[i]*(1-population[i]/K)
        } else {
            if(i <= tau) {
                population[i] <- population[i-1]+dNdt[i-1]
                dNdt[i] = r*population[i]*(1-population[i]/K)
             } else {
                 population[i] = population[i-1]+dNdt[i-1]
                 dNdt[i] = r*population[i]*(1-population[i-tau]/K)
             }
        }
    }

Bem temos dois if nessa simulação discreta (para noção da diferença entre discreto e continuo, olhe aqui e aqui).

O primeiro, aqui

1
2
3
4
        if(i == 1){
            population[i] <- N0
            dNdt[i] <- r*population[i]*(1-population[i]/K)
        }

A gente so está tratando da geração 0.

1
2
3
4
5
6
        } else {
            if(i <= tau) {
                population[i] <- population[i-1]+dNdt[i-1]
                dNdt[i] = r*population[i]*(1-population[i]/K)
             } 
          }

Depois aqui a gente ta tratando do inicio da dinâmica, lembre-se que falar que temos um lag de 10 gerações, ou 10 tempos, significa que o tamanho populacional de agora afeta a taxa de crescimento de daqui 10 gerações, so que não da para isso acontecer na primeira geração, nem na segunda, nem na decima, porque não tem 10 gerações para trás ainda, então temos que dar um jeito nisso, que aqui é usar o crescimento logístico sem lag nesse inicio.

1
2
3
4
             } else {
                 population[i] = population[i-1]+dNdt[i-1]
                 dNdt[i] = r*population[i]*(1-population[i-tau]/K)
             }

E por ultimo, a gente trata do caso quando começamos a poder implementar o lag, no caso do lag 10, a partir do momento que ja temos 10 gerações para trás.

Veja que a diferença está somente aqui na formula

 (1-{N\over K})

O efeito da diminuição da taxa de crescimento. A gente usa a população do tamanho lag para traz

 (1-{N_{t-lag}\over K})

ou seja

  N_{t+1} = r\cdot N_{t} \cdot (1- \alpha * N_{t-lag} )

A taxa de crescimento, vezes o tamanho populacional agora, mas o efeito na diminuição da taxa de crescimento é da população t-lag gerações atrás.

Bem eu vou colocar tudo isso numa função para facilitar a simulação, que é o que vai estar no script la em baixo, aqui a gente so estava dando uma olhada.

Mas o resultado disso, para um lag de 10 gerações fica o seguinte, sobre a dinâmica populacional.

01

A população consegue passar a capacidade suporte, no inicio, quando a responsabilidade bate, e vem o efeito dessa população alta, a taxa de crescimento fica negativa e ela cai, depois com a população abaixo da capacidade suporte, ela sobe denovo, mas esses altos e baixos vai diminuindo, até haver uma estabilização na capacidade suporte.

A gente pode olhar a taxa de crescimento de acordo com o tamanho populacional

02

Veja que legal, eu achei essa figura muito legal, a gente fica rodando perto taxa de crescimento zero, que é o que acontece quando a gente chega na capacidade suporte, veja que quando mais próximo da estabilidade, mais perto de parar em zero a gente fica rodando. Algo similar ao que acontece com o crescimento per capita, a contribuição que cada individuo da população pode fornecer.

03

Certo então temos o efeito dependente de densidade atrasado na dinamica populacional.

04

Que ela vai ter um momento inicial diferente do crescimento logístico, mas vai tender a uma estabilidade também, na capacidade suporte do ambiente.

Mas será que isso vale para todos os valores de atraso? Bem se a gente diminuir o lag para 3 por exemplo.

05

A gente vai se aproximar do crescimento logístico, o comportamento do crescimento logístico, onde a gente só diminui a taxa de crescimento per capita e temos um crescimento grande no começo e depois que chegamos a metade do tamanho populacional, começamos a ter uma queda na taxa de crescimento.

Mas e se a gente aumentar bastante o lag? Algo como lag = 20.

06

Que legal, a gente entra num ciclo, onde a população sobe muito acima da capacidade suporte, e depois cai bruscamente quando esse efeito vem a tona, mas depois como fica com um tamanho pequeno, ela vai e cresce bastante denovo, e fica girando nesse ciclo, diferente de um lag médio quando ela oscilava no começo mas depois estabilizava.

Bem é isso ai, eu estou vendo umas coisas legais em um curso do Edx chamado Quantitative Biology Workshop, esse modelo estava perdido no meio de um exercício, mas eu achei muito legal. No começo o curso estava meio chato, ficar fazendo exercícios em matlab (eu que nunca tenho grana sobrando, quiça para comprar o matlab) não me atrai muito, mas eles tem exercícios com simulações bem legais, e é apresentado aula dos professores que dão esse curso la no mit no final da semana, que é bem legal de ver também, mas então, 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.

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
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
#Post
#http://recologia.com.br/2015/01/crescimento-populacional-com-atraso-para-o-efeito-de-dependencia-da-densidade
 
#Parâmetros de entrada
N0 <- 1 #População inicial
r  <- 0.1 #Taxa de crescimento Populacional
K  <- 100 #Capacidade Suporte
t  <- 200 #Número de gerações
tau<- 20  #Atraso de resposta para o efeito dependente de densidade
 
simu_loglag<-function(N0=1,r=0.1,K=100,t=200,tau=10) {
    generation = seq(0,t,by=1)
    dNdt = rep(0,0,t)
    population = rep(0,0,t)
 
    for(i in 1:length(generation)) {
        #Aqui é o inicio da simulação, como não temos nada de população
        #ainda, usamos a população inicial
        if(i == 1){
            population[i] <- N0
            dNdt[i] <- r*population[i]*(1-population[i]/K)
        } else {
            if(i <= tau) {
                #Aqui, é o caso onde o lag ainda é maior que o tamanho da
                #populacao, lembre-se que se a gente tem um lag de 10, até o
                #tempo 10, não da para olhar 10 tempos atrás, diminuir 10 do
                #i, então a gente tem que usar o crescimento sem lag no
                #começo
                population[i] <- population[i-1]+dNdt[i-1]
                dNdt[i] = r*population[i]*(1-population[i]/K)
             } else {
                 population[i] = population[i-1]+dNdt[i-1]
                 dNdt[i] = r*population[i]*(1-population[i-tau]/K)
             }
        }
    }
    saida<-data.frame(generation,population,dNdt)
    return(saida)
}
 
 
simu_lag10<-simu_loglag(tau=10)
 
#
plot(population~generation,xlab="Tempo",ylab="Tamanho da população",
     data=simu_lag10,type="l",lwd=2,frame=F)
 
#
plot(dNdt~population,xlab="Tamanho da população",
     ylab="Crescimento populacional",data=simu_lag10,type="l",lwd=2,frame=F)
 
#
plot((dNdt/population)~population,data=simu_lag10,type="l",lwd=2,frame=F,
     xlab="Tamanho da população",ylab="Crescimento populacional per capita")
 
#Lag = 10
layout(matrix(1:3,ncol=1,nrow=3))
plot(population~generation,xlab="Tempo",ylab="Tamanho da população",
     data=simu_lag10,type="l",lwd=1,frame=F,
     main="Tamanho populacional ao longo do tempo, Lag = 10")
plot(dNdt~population,xlab="Tamanho da população",
     ylab="Crescimento populacional",data=simu_lag10,type="l",lwd=1,frame=F,
     main="Crescimento vs Tamanho populacional")
plot((dNdt/population)~population,data=simu_lag10,type="l",lwd=1,frame=F,
     xlab="Tamanho da população",ylab="Crescimento populacional per capita"
     ,main="Crescimento per capita vs Tamanho populacional")
 
#Lag = 3
simu_lag3<-simu_loglag(tau=3)
layout(matrix(1:3,ncol=1,nrow=3))
plot(population~generation,xlab="Tempo",ylab="Tamanho da população",
     data=simu_lag3,type="l",lwd=1,frame=F,
     main="Tamanho populacional ao longo do tempo, Lag = 3")
plot(dNdt~population,xlab="Tamanho da população",
     ylab="Crescimento populacional",data=simu_lag3,type="l",lwd=1,frame=F,
     main="Crescimento vs Tamanho populacional")
plot((dNdt/population)~population,data=simu_lag3,type="l",lwd=1,frame=F,
     xlab="Tamanho da população",ylab="Crescimento populacional per capita"
     ,main="Crescimento per capita vs Tamanho populacional")
 
#Lag = 20
simu_lag20<-simu_loglag(tau=20)
layout(matrix(1:3,ncol=1,nrow=3))
plot(population~generation,xlab="Tempo",ylab="Tamanho da população",
     data=simu_lag20,type="l",lwd=1,frame=F,
     main="Tamanho populacional ao longo do tempo, Lag = 20")
plot(dNdt~population,xlab="Tamanho da população",
     ylab="Crescimento populacional",data=simu_lag20,type="l",lwd=1,frame=F,
     main="Crescimento vs Tamanho populacional")
plot((dNdt/population)~population,data=simu_lag20,type="l",lwd=1,frame=F,
     xlab="Tamanho da população",ylab="Crescimento populacional per capita"
     ,main="Crescimento per capita vs Tamanho populacional")

R Orientado a Objetos e um exemplo S4

A gente já falou de programação orientado a objetos aqui, usando como exemplo classes do tipo S3.

Mas existe outra forma de definir classe no R que é a S4. Essa é uma forma que foi desenvolvida posteriormente ao método S3, e tem mais uma cara de orientação a objetos, no sentido que ela é mais estrita nas definições dela.

Bem vamos começar redefinindo nossa classe armadilha de sementes.

01

Então, primeiro definimos nossas classes, setando quais os seus atributos

1
2
3
4
#definindo classe
setClass("armadilhaDeSementes",representation(distancias = "numeric",
                                       contagem.sementes = "numeric",
                                       area.armadilha = "numeric"))

Depois nos podemos fazer um construtor, uma função que vai construir objetos dessa classe. Lembrando que essa parte é interessante e importante, pois é a hora que podemos validar os objetos, para que os atributos não tenham valores que não façam sentido.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
#construtor
setMethod("initialize","armadilhaDeSementes",
          function(.Object,
                   distancias = numeric(0),
                   contagem.sementes = numeric(0),
                   area.armadilha = numeric(0)) {
    if (length(distancias) != length(contagem.sementes)) {
        stop("Quantidade de distâncias e contagens é diferente.")
    }
    if (length(area.armadilha) != 1) {
        stop("Área da armadilha é ambigua.")
    }
    .Object@distancias <- distancias
    .Object@contagem.sementes <- contagem.sementes
    .Object@area.armadilha <- area.armadilha
    .Object
})

Agora a primeira diferença entre as classes s3 e s4 visível, é que objetos da classes s4 so podem ser construídos pela função new.

> s1 <- new("armadilhaDeSementes",distancias = 1:4, + contagem.sementes = c(4, 3, 2, 0), + area.armadilha = 0.0001) >

Mas veja que para ficar com uma cara do que vemos no R, basta fazer uma função que chama a função new.

1
2
3
4
#criando uma função para criar objetos s4
armadilhaDeSementes <- function(...) {
    new("armadilhaDeSementes",...)
}

Agora a gente pode criar objetos da forma convencional, que sempre vemos.

> s2<-armadilhaDeSementes(distancias=1:4,contagem.sementes=c(4, 3, 2, 0), + area.armadilha=0.0001) > class(s2) [1] "armadilhaDeSementes" attr(,"package") [1] ".GlobalEnv" >

Veja que esse é um bom exemplo do uso do …, quando uma função tem … como argumento, o que ela faz é pegar tudo que está como argumento, depois dos definidos nela, e passa pra frente, aqui a função armadilhaDeSementes pega tudo que tem de argumento dela e manda direto para o new(“armadilhaDeSementes”,…)

Outra coisa, é que para ver quais atributos uma classe tem, a gente usa a função slotNames

> slotNames(s1) [1] "distancias" "contagem.sementes" "area.armadilha"

E acessa os valores usando o arroba, ao invés do $

> #acessando os atributos > s1@distancias [1] 1 2 3 4

A gente vai definir métodos usando setMethod, colocando como assinatura qual a classe de que ela se trata.

1
2
3
4
5
6
7
8
9
10
11
12
#definindo métodos
setMethod("show",signature(object = "armadilhaDeSementes"),
          function(object) {
              str(object)
          }
          )
#
setMethod("mean",signature(x = "armadilhaDeSementes"),
          function(x, ...) {
              weighted.mean(x@distancias,w = x@contagem.sementes)
          }
          )

A essa funções vão ser atribuídas a essa classe direto, sem ter usar o esquema do ponto como nas classes S3. Aqui também existe mais uma nuance. Toda função deveria retornar algo, mesmo que retorne nada, um null, mas as funções no R, quanto nada é definido para um return, ela retorna a “coisa” da ultima linha da função, por isso o str(object) e weighted.mean(x@distancias,w = x@contagem.sementes) não estão dentro de um return(), mas se estivessem, teríamos a mesma coisa no fim das contas.

> s1 Formal class 'armadilhaDeSementes' [package ".GlobalEnv"] with 3 slots ..@ distancias : int [1:4] 1 2 3 4 ..@ contagem.sementes: num [1:4] 4 3 2 0 ..@ area.armadilha : num 1e-04 > mean(s1) [1] 1.777778

Cuidado, pois dessa forma, methods mostra apenas funções S3 (espero que não esteja falando besteira aqui).

> methods(mean) [1] mean.Date mean.default mean.difftime mean.POSIXct mean.POSIXlt

Precisamos usar o showMethods para ver os métodos S4

> showMethods("mean") Function: mean (package base) x="ANY" x="armadilhaDeSementes"

O showMethods também pode mostrar os métodos associados a uma classe, com uma determinada assinatura.

> showMethods(classes = "armadilhaDeSementes") Function: initialize (package methods) .Object="armadilhaDeSementes" Function: mean (package base) x="armadilhaDeSementes" Function: show (package methods) object="armadilhaDeSementes"

E a função getMethod mostra o código fonte da função, para classes S3 basta a gente digitar o nome da função sem parenteses, mas para funções S4 a gente usa o getMethod.

> getMethod("mean", "armadilhaDeSementes") Method Definition: function (x, ...) { weighted.mean(x@distancias, w = x@contagem.sementes) } Signatures: x target "armadilhaDeSementes" defined "armadilhaDeSementes"

Bem é isso ai, um tutorial legal para ver mais coisas de classes S4, é o tutorial do Hadley Wickham aqui (alias esse cara tem excelente tutoriais, vale a pena navegar na página dele aqui, e hey, ele fez o ggplot2), 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:
Owen Jones, Robert Maillardet & Andrew Robinson 2009 – Introduction to Scientific Programming and Simulation Using R. CRC Press 499pp

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
#definindo classe
setClass("armadilhaDeSementes",representation(distancias = "numeric",
                                       contagem.sementes = "numeric",
                                       area.armadilha = "numeric"))
 
#construtor
setMethod("initialize","armadilhaDeSementes",
          function(.Object,
                   distancias = numeric(0),
                   contagem.sementes = numeric(0),
                   area.armadilha = numeric(0)) {
    if (length(distancias) != length(contagem.sementes)) {
        stop("Quantidade de distâncias e contagens é diferente.")
    }
    if (length(area.armadilha) != 1) {
        stop("Área da armadilha é ambigua.")
    }
    .Object@distancias <- distancias
    .Object@contagem.sementes <- contagem.sementes
    .Object@area.armadilha <- area.armadilha
    .Object
})
 
#instanciando um objeto
s1 <- new("armadilhaDeSementes",distancias = 1:4,
          contagem.sementes = c(4, 3, 2, 0),
          area.armadilha = 0.0001)
 
#criando uma função para criar objetos s4
armadilhaDeSementes <- function(...) {
    new("armadilhaDeSementes",...)
}
 
#instanciando um objeto
s2<-armadilhaDeSementes(distancias=1:4,contagem.sementes=c(4, 3, 2, 0),
                 area.armadilha=0.0001)
class(s2)
 
#atributos da classe
slotNames(s1)
 
#acessando os atributos
s1@distancias
 
 
#definindo métodos
setMethod("show",signature(object = "armadilhaDeSementes"),
          function(object) {
              str(object)
          }
          )
#
setMethod("mean",signature(x = "armadilhaDeSementes"),
          function(x, ...) {
              weighted.mean(x@distancias,w = x@contagem.sementes)
          }
          )
 
#usando método show
s1
 
#método generico
mean(s1)
 
#métodos s3 de mean
methods(mean)
 
#métodos com assinatura s4
showMethods("mean")
 
#métodos da classe
showMethods(classes = "armadilhaDeSementes")
 
#vendo o código do método mean da classe armadilhaDeSementes
getMethod("mean", "armadilhaDeSementes")

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))