Diferença entre correlação e Regressão

Recentemente uma amiga minha me perguntou sobre qual a diferença entre regressão e correlação, para fazer inferência sobre relação entre duas variáveis. Pode parecer uma pergunta simples mas até eu parei para pensar nessa.

A regressão linear é uma velha conhecida na biologia.

01

Supondo que a gente tem duas variáveis, um x preditor e um y como resposta.
Na correlação, com os mesmo dados a gente também calcula a relação entre esses x e y, então o que é diferente e o que é igual?

Bem a regressão, como sabemos, é baseada nos mínimos quadrados.

02

Ou seja, minimizar a distância, no eixo y (veja que não é a distância diagonal, eu já pensei isso errado), mas a distância ao quadrado, para que não tenhamos que lidar com números negativos, todo número ao quadrado da um número da positivo.

Então a gente quer o seguinte:

d=y-\hat{y}

Aqui, y é o valor que observamos e \hat{y} é o valor previsto pelo modelo, seja la qual for, isso para os mínimos quadrado, por isso falei que é a diferença no eixo y apenas. Agora na regressão linear a gente tem a famosa formula \hat{y}=a+bx, então a gente pode substituir essa fórmula na equação dos mínimos quadrados anteriores.

Então :
d=y-\hat{y}
Vira
d=y-(a+bx)
E distribuindo o sinal negativo ali
d=y-a-bx

Ok, agora como temos muitas amostras, a nossa medida de ajuste são as distâncias, então são as somas quadradas desses valores
\sum\limits_{i=1}^n d^2= \sum\limits_{i=1}^n (y-a-bx)^2

Se a gente chutar uma fórmula qualquer, bem vamos tentar alguns valores para a (intercepto) e b (inclinação)

03_1

Temos essa superfície que é o ajuste, mas agora presta atenção nessa figura, veja que o melhor ajuste para a inclinação é sempre o mesmo independente do intercepto, ele só está “deitado” para a melhor inclinação, mas se a gente fizer alguns cortes nessa superfície, tipo pegando um valor fixo de intercepto, a gente vê o seguinte.

03_2

Veja que a inclinação vai cair sempre no mesmo lugar, e isso é uma função do segundo grau, que a gente consegue achar o mínimo com cálculo, de uma olhada aqui e aqui.

Como comenta o Crawley no livro dele, no R a gente vai usar os famosos 5 no R, que são \sum x,\sum x^2,\sum y,\sum y^2 e \sum xy

> sum(x);sum(x^2);sum(y);sum(y^2);sum(x*y) [1] 170.1989 [1] 1152.791 [1] 170.9315 [1] 995.8592 [1] 973.15

Não é a toa que álgebra linear está em todo lugar, aqui podemos calcular essas quantidades usando álgebra linear também

> XY <- cbind(1,x,y) > t(XY) %*% XY x y 30.0000 170.1989 170.9315 x 170.1989 1152.7910 973.1500 y 170.9315 973.1500 995.8592

Agora para chegar inclinação que produz o melhor ajuste, precisamos do soma dos quadrados corrigido e da soma dos produtos, isso vai ser:

SSX = \sum x^2 - \frac{(\sum x)^2}{n}
SSXY = \sum xy - \frac{(\sum x)(\sum y)}{n}

E a inclinação vai ser

b=\frac{SSXY}{SSX}

Repetindo isso no R temos:

> SSX<-sum(x^2)-sum(x)^2/n > SSXY <- sum(x*y)-(sum(x)*sum(y))/n > inclinacao<-SSXY/SSX > inclinacao [1] -1.103325

Agora a partir da melhor inclinação, podemos calcular o intercepto que vai ser

a=\bar{y}-b\bar{x}

Veja que é só a formula da regressão, isolando o a, e a reta de regressão sempre vai cruzar o ponto formado pela média de x e a média de y.

> intercepto<-mean(y)-inclinacao*mean(x) > intercepto [1] 5.659964

E fazendo a regressão usando o lm do R temos:

> lm(y~x) Call: lm(formula = y ~ x) Coefficients: (Intercept) x 5.660 -1.103

Mesmos valores, salvo o arredondamento, que o lm so mostra 3 casas.

> round(intercepto,digits=3) [1] 5.66 > round(inclinacao,digits=3) [1] -1.103

Agora qual a relação disso com correlação?
Bem na correlação a gente vai ver o quanto uma variável varia junto com a outra.

04

Veja que temos mais pontos nos quadradinhos azuis e vermelhos que qualquer outro. A correlação vai ser basicamente a inclinação normalizada, independente dos valores das variáveis, vai ser a inclinação variando em unidades de desvio padrão.

O Coeficiente de correlação de pearson vai ser dado por:

r=\frac{1}{(n-1)}\sum\limits_1^n \frac{(x_i-\bar{x})}{sd(x)} \frac{(y_i-\bar{y})}{sd(y)}

Só reproduzir isso direto no R

> 1/(n-1)*sum(((x-mean(x))/sd(x))*((y-mean(y))/sd(y))) [1] -0.9480939

E podemos conferir com a função cor dele

> cor(x,y) [1] -0.9480939

So que esse número então é proporcional a inclinação b

b=cor(y,x) \frac{sd(y)}{sd(x)}

Podemos testar no nosso exemplo

> inclinacao [1] -1.103325 > cor(y,x)*(sd(y)/sd(x)) [1] -1.103325

Por isso, por inferência estatística, se você conclui que y esta correlacionado com x, você também encontra significância na regressão.

Resumidamente

O coeficiente de regressão normalizado é o mesmo que o coeficiente de correlação de Pearson
O quadrado do coeficiente de Pearson é o mesmo que o R^2 na regressão linear

> summary(lm(y~x)) Call: lm(formula = y ~ x) Residuals: Min 1Q Median 3Q Max -2.14294 -0.72029 0.01953 0.70563 1.88121 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 5.65996 0.46646 12.13 1.15e-12 *** x -1.10333 0.06993 -15.78 1.83e-15 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.9876 on 28 degrees of freedom Multiple R-squared: 0.8989, Adjusted R-squared: 0.8953 F-statistic: 248.9 on 1 and 28 DF, p-value: 1.835e-15 > round(cor(y,x)^2,digits=4) [1] 0.8989

Nem a regressão linear nem a correlação respondem a questão sobre casualidade diretamente. Mas a equação da regressão pode ser usada para calcular valores previstos em y baseados em x, mas logicamente, ambos os casos estamos prevendo que se uma coisa muda em x, o y também muda, mas na correlação, como os valores são em desvios padrão, unidades de desvio padrão, tanto a correlação de x e y como de y e x são o mesmo valor, enquanto para a regressão, precisamos do preditor no eixo x.

Ou você tem muito de um e do outro (ou o contrário, muito de um pouco do outro) ou nada. A alternativa a isso, a não correlação é algo assim.

05

Então tanto a correlação com a regressão podem ser usadas para fazer uma previsão, mas com a regressão a gente termina com uma equação para usar diretamente nos valores de x, enquanto com a correlação, a gente tem a variação em unidades de desvio padrão. Que pode ser interessante para comparações quando temos experimentos com unidades diferentes, não comparáveis diretamente. Então o que usar acho que depende muito do contexto, mas no final das contas, podemos caminhar da inclinação da regressão para o coeficiente de correlação de Pearson e vice versa.

Historicamente, correlação foi inventada pelo Karl Pearson, que foi orientado do Francis Galton, que começou com a ideia da regresão a média. Mas o Karl Pearson que começou a usar os valores p para seu teste de chi-quadrado, depois o Ronald Fisher que populariou o uso do valor p, mas isso começou la no Laplace, que olhando orbitas de planetas já usava Verossimilhança e a ideia do valor p.

Eu imagino, se para a correlação é mais fácil conseguir chegar ao valor p, do que na regressão linear, lembrando que todo mundo fazia as coisas na mão no passado, não é a tranquilidade que é hoje com ajuda dos computadores. Mas então no fim, acho que a escolha entre correlação e regressão, para certas perguntas pode ser mais cultural, de o que você aprendeu, com o que você se sente mais seguro de expor num manuscrito, de comunicar as pessoas, já que no final das contas, a inclinação e a correlação são proporcionais.

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

E antes que eu mês esqueça, a nature liberou vários artigos dela sobre bioestatística nesse link:
http://www.nature.com/collections/qghhqm/
Achei bem legal.

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
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
109
110
111
112
#post:
#http://recologia.com.br/2015/02/diferenca-entre-correlacao-e-regressao
 
#Gerando um exemplo
set.seed(123)
n<-30
x<-runif(n,1,10)
y<-rnorm(n,5-x)
 
#figura 1
plot(y~x,frame=F)
 
 
#figura 2
plot(y~x,pch=19,frame=F)
modelo<-lm(y~x)
abline(modelo)
for(i in 1:n) { lines(cbind(x,x)[i,],cbind(predict(modelo),y)[i,]) }
 
 
#Figuras 3.1 e 3.2
inclinacoes <- seq(-10,10,0.1)
interceptos <- seq(3,7,0.1)
 
superficie<-matrix(0,nrow=length(interceptos),ncol=length(inclinacoes))
for(i in 1:length(interceptos)) {
    for(j in 1:length(inclinacoes)) {
        superficie[i,j]<-sum((y - interceptos[i] - inclinacoes[j]*x)^2)
    }
}
#3.1
persp(interceptos,inclinacoes,superficie,theta = 20, phi = 20)
 
SSE <- function(i) sum((y - 5 - inclinacoes[i]*x)^2)
#3.2
par(mfrow=c(2,2))
for(i in 3:6) {
plot(inclinacoes,sapply(1:length(inclinacoes),SSE),type="l",
     main=paste("Intercepto = ",i),xlab="inclinação b",
     ylab="soma dos resíduos quadrados",lty=3,cex=3)
}
 
 
#Famosos 5
sum(x);sum(x^2);sum(y);sum(y^2);sum(x*y)
XY <- cbind(1,x,y)
t(XY) %*% XY
 
SSX<-sum(x^2)-sum(x)^2/n
 
SSY<-sum(y^2)-sum(y)^2/n
 
SSXY <- sum(x*y)-(sum(x)*sum(y))/n
 
#Definindo inclinação e intercepto
inclinacao<-SSXY/SSX
inclinacao
 
intercepto<-mean(y)-inclinacao*mean(x)
intercepto
 
#Inclinação a partir da correlação
inclinacao
cor(y,x)*(sd(y)/sd(x))
 
#R^2
summary(lm(y~x))
round(cor(y,x)^2,digits=4)
 
 
#Comparando valores encontrados
lm(y~x)
round(intercepto,digits=3)
round(inclinacao,digits=3)
 
#figura 4
plot(y~x,frame=F,type="n")
abline(v=mean(x),lty=3,cex=2)
abline(h=mean(y),lty=3,cex=2)
 
indices<-x>mean(x) & y>mean(y)
points(x[indices],y[indices],col=1,pch=19)
indices<-x<=mean(x) & y>mean(y)
points(x[indices],y[indices],col=2,pch=19)
indices<-x<=mean(x) & y<=mean(y)
points(x[indices],y[indices],col=3,pch=19)
indices<-x>mean(x) & y<=mean(y)
points(x[indices],y[indices],col=4,pch=19)
dev.off()
 
#modelo sem correlação
x<-runif(n,1,10)
y<-rnorm(n,mean(x))
 
#figura 5
plot(y~x,frame=F,type="n")
abline(v=mean(x),lty=3,cex=2)
abline(h=mean(y),lty=3,cex=2)
 
indices<-x>mean(x) & y>mean(y)
points(x[indices],y[indices],col=1,pch=19)
indices<-x<=mean(x) & y>mean(y)
points(x[indices],y[indices],col=2,pch=19)
indices<-x<=mean(x) & y<=mean(y)
points(x[indices],y[indices],col=3,pch=19)
indices<-x>mean(x) & y<=mean(y)
points(x[indices],y[indices],col=4,pch=19)
dev.off()
 
#Calculando a correlação
cor(y,x)
1/(n-1)*sum(((y-mean(y))/sd(y))*((x-mean(x))/sd(x)))

Evolução – Fisherian Optimality Models – Parte 1

Com certeza, uma disciplina muito legal dentro da biologia é a evolução.

Como diz o celebre trabalho do Theodosius Dobzhansky de 1973:

“Nothing in Biology Makes Sense Except in the Light of Evolution”

“Nada na biologia faz sentido a não ser sobre a luz da evolução”

Mas nem tudo é alegria, o que essa frase anima a gente, as vezes outras são meio perturbadoras, como a do John Maynard Smith:

“If you can’t stand algebra, keep out of evolutionary biology”

“Se você não suporta álgebra, fique fora da biologia evolutiva”

Essa frase estava no inicio do livro Algebraic Statistics for Computational Biology do Lior Pachter e Bernd Sturmfels. Eu nunca dei conta de ler esse livro inteiro, porque é meio tedioso as vezes. Mas recentemente eu comecei a ler um livro muito legal chamado Modeling Evolution – an introduction to numerical methods do Derik a. Roff.

É o tipo de livro que eu acho legal, porque no próprio prefácio, o autor fala que o livro veio de uma disciplina que ele dava em Harvard, mas é o livro que ele queria ter lido na graduação. Essa são palavras chave pra gente ver um livro que o cara escreveu com gosto, e que não pega nem tão pesado na biologia como na matemática.

O livro coloca bem como começar do zero modelos em evolução, e como ir procedendo, desde os mais simples ao mais complexos. Eu pelo menos nunca tinha estudado muito a fundo como é a análise de modelos em evolução, o que é ruim porque a gente acaba simplesmente acreditando em tudo que le, ja que a análise passa a ser uma barreira.

Mas já que estamos aqui, vamos dar uma olhada em como são os modelos mais básicos, da forma como o Ronald Fisher começou a analisar as coisas em evolução, depois das explicações do livro, deu para dar uma melhorada. O próprio Darwin deu uma baqueada no fim da sua vida, pensando que sua teoria da evolução tinha talvez muitas coisas que não conseguia explicar, mas depois veio o Fisher e companhia com a teoria neodarwinista e começaram a matar a charada.

Os Fisherian Optimality Models são baseados na abordagem do Fisher para estudar questões evolutivas nas quais ele pegou o parâmetro malthusiano, r, como a medida de fitness. Essa medida pode ser ruim quando uma distribuição etária estável não pode ser assumida, existe uma dependência da densidade ou a frequência é um fator importante na evolução das características em questão, a população é caótica ou o fitness depende de interações sociais. Mas mesmo com tantas restrições, esse tipo de modelo é muito interessante para a análise de questões evolutivas e na geração de predições testáveis, alias essa é uma critica comum quando a pesquisa no Brasil, trabalhos que não tem predições a priore, não testam hipóteses, é coletar dados e ver o que da pra fazer. Agora da para tentar começar a usar esse tipo de abordagem para formular hipóteses, predições a priore.

Mas continuando, três dos parâmetros populacionais mais comumente usados nesse tipo de modelo são a taxa de crescimento populacional r, sucesso reprodutivo ao longo da vida (ou sucesso reprodutivo liquido) R0 e um componente do fitness o qual na sua maximização também maximiza o fitness. Dos três primeiros, o primeiro é o mais comumente usado. O parâmetro malthusiano é a taxa de crescimento da população que está com uma distribuição de idades estável (vou ficar devendo um post sobre isso), e é dado pela equação:

\int_0^{\inf} e^{-rt}l(x)m(x)dx=1

Onde o t é a tempo, l(x) é a taxa de sobrevivência no tempo x, e m(x) é o número de nascimentos por fêmea na idade x (normalmente a razão sexual é 1:1, mas isso não é obrigatório). Na ausência da dependência de densidade, é meio obvio que mutações que aumentem r devem aumentar sua frequência como a gente ja viu aqui. Se existe uma dependência da densidade, então a medida apropriada de fitness é o número relativo de indivíduos que passam do período no qual a dependência da densidade age. Por exemplo, suponha que a dependência da densidade ocorra na fase imatura e que nos queremos calcular o fitness de dois tipos de indivíduos, A e B. Se o número de adultos de A excede o número de adultos de B então A é necessariamente o com melhor fitness. A análise de tal situação pode se tornar mais complicada se o sucesso nesse período não depende simplesmente da densidade, mas também na frequência dos tipos na população.
Se a população não aumenta em tamanho, uma medida plausível do fitness é a taxa reprodutiva liquida.

R_0 = \int_0^{\inf} l(x)m(x)dx

Essa medida e algumas vezes é chamado de expectativa de sucesso reprodutivo ao longo da vida (“expected lifetime reproductive success”). A premissa aqui é que se A tem uma R_0 maior que B então A vai ter um fitness maior. Dependência da densidade é assumida implicitamente como não afetando as características de interesse. Quando existem flutuações estocásticas que movem a população para fora de uma estrutura estável de idade, nem r nem o R_0 podem ser assumidos como medidas apropriadas para medir o fitness.

Em alguns casos, particularmente estudos de comportamento, nosso interesse é em uma característica em particular e para simplificar a análise nos assumimos que a maximização de alguns componentes da historia de vida é equivalente a maximização do fitness. Por exemplo, nos podemos querer determinar a alocação ótima de tempo de forrageamento entre manchas que variam quanto a qualidade de seus recursos, a gente comentou sobre isso aqui, aqui e aqui. Mas uma premissa comum é que a maximização da obtenção de recursos por unidade de tempo é equivalente a maximização do fitness, mas precisamos de cuidado, porque ignorar outros componentes da historia de vida pode produzir conclusões erronias. Dois exemplos famosos são o paradoxo de Cole e a hipótese de Lack.

No seu clássico artigo de 1954, Cole sugeriu que havia um aparente paradoxo porque

“For an annual species, the absolute gain in intrinsic population growth which could
be achieved by changing to the perennial reproductive habit would be exactly equivalent to adding one individual to the average litter size”

“Para uma espécie anual, o ganho absoluto intrínseco no crescimento populacional que poderia ser obtido alterando o habito reprodutivo perenial seria exatamente equivalente a adicionar um individuo ao tamanho médio da ninhada”

Se a conclusão de Cole fosse verdadeira, nos esperaríamos que plantas perenes deveriam ser raras, o que certamente não é o que observamos. O problema foi que cole não incluiu diferenças na sobrevivência entre adultos e jovens, sua análise assumia implicitamente taxas de sobrevivência de 1 para ambos os estágios.

Ja Lack(1947) teve a ideia que em aves…

“the average clutch-size is ultimately determined by the average maximum number of young which the parents can successfully raise in the region and season in question”

“o tamanho médio da ninhada é ultimamente determinado pelo numero médio de jovens que os pais podem criar com sucesso na região em questão”

A hipótese de Lack assumia que apenas interações negativas dependentes de densidade entre irmãos dentro de um ninho eram importante para predizer o tamanho ótimo de filhos a se criar em um ninho para os pais, logo a quantidade de ovos que vemos em um ninho é o máximo que os pais podem criar. Lack não observou a possibilidade que a sobrevivência na fase a adulta dos filhotes da ninhada pode ser afetada pelo numero de irmãos com o qual cresceram. Experimentos manipulando o tamanho de ninhadas em aves, mamíferos, repteis, peixes e plantas demostrou efeitos negativos de uma ninhada maior na sobrevivência futura dos adultos e no seu sucesso reprodutivo. A incorporação de tais efeitos prediz que o tamanho de ninhada ótima vai ser menor que o valor proposto por Lack.

Métodos de análises:

O foco da análise é em condições de equilíbrio e não na trajetória evolucionaria até ele (mas a trajetória pode ser analisada também, espero abordar isso depois quando ler todo o livro de Evolução). Aqui a dependência da densidade não é explicitamente considerada. Dependência da frequência é também assumida como ausente. A formulação do modelo que gostaríamos de chegar é

W = f( \theta_1, \theta_2, \dots, \theta_k, x_1, x_2, \dots , x_n)

onde W é o fitness, \theta_1,\theta_2,\dots,\theta_k são os parâmetros e x_1 ,x_2 , \dots , x_{n} são as características.
Essa função é para ser lida como Fitness é uma função de k parâmetros e n características (“Fitness is a function of k parameters and n traits.”).

Uma regra geral é manter o número de parâmetros e características no mínimo. Suponha que tenhamos 5 parâmetros e queremos avaliar a performance de todas as combinações, dividindo cada parâmetro em 10 partes, que é uma número razoável de divisões, isso nos da 10^5, 100000 combinações para examinar. Apesar de ser possível (usando o computador), essa pode não ser o caminho preferível, se puder ser evitado.

A função de fitness vai invariavelmente ser feita de um número de funções componentes, como veremos no exemplo daqui a pouco, no qual o fitness é o produto da fecundidade e da sobrevivência e ambas são função do tamanho do corpo. Essas funções podem produzir curvas suaves de fitness que podem ser analisadas usando calculo (integral e diferencial) ou a função pode ter descontinuidades ou ser difícil de ser analisadas por algum outro motivo (uma superfície irregular) usando o calculo.

É preciso cuidado também examinando modelos com descontinuidades e lugares no qual os componentes do modelo podem tomar valores fisicamente impossíveis (irrelevantes do ponto de vista biológico). No exemplo que a gente vai olhar aqui é dado por uma função linear negativa do tamanho do corpo, o que significa que abaixo de um tamanho particular a sobrevivência vai ficar negativa, o que não é possível no mundo real. No exemplo isso não é um problema porque o fitness será negativo nesse caso. Entretanto, suponha que o modelo contenha o produto de duas funções de sobrevivência que podem ser matematicamente menores que zero. Agora é possível que os dois valores negativos produzam um valor positivo e o valor de fitness não pode ser visto como incorreto.

No fim, existem três casos gerais em que acabamos sempre caindo.

Primeiro: fitness pode ser escrito como uma função de características de interesse e a diferenciação não é um problema.

Segundo: fitness pode ser escrito como uma função de características de interesse mas a diferenciação é problemática.

Terceiro: fitness não pode ser escrito como uma função na qual o fitness pode ser isolado.

Métodos de análises:
W=f(\theta_1,\theta_2,\dots,\theta_k,x_1,x_2,\dots,x_n) é bem comportado.

Assumindo uma função de fitness diferenciável (uma superfície suave de fitness), nos achamos o valores ótimos para as características diferenciando a função com respeito as características e igualando o tudo a zero. Por exemplo, suponha que a função de fitness é uma função quadrada.

W=f(\theta_1,\theta_2,\theta_3,x) = \theta_1 + \theta_2x + \theta_3 x^2

Então:

\frac{dW}{dx} = \theta_2 + 2 \theta_3 x

\frac{dW}{dx} =  0 quando  x = \frac{-\theta_2}{2 \theta_3}

Note que um fitness intermediário requer que  \theta_3 \textless 0 e \theta_2 > 0 . Podemos ficar tentados a esperar apenas que \frac{-\theta_2}{2 \theta_3} \textless 0, entretanto, a condição \theta_3 > 0 e \theta_2 \textless 0 define uma função convexa de forma que o ponto de flexão seja um mínimo e não uma máxima. Sempre que possível, a função deve ser plotada numa figura para garantir que um ponto apropriado de flexão exista. Se a função é complexa e contem muitos termos, a diferenciação pode se tornar muito tediosa e cuidado é necessário para que termos não sejam acidentalmente omitidos ou sinais trocados.

Código de R para derivadas:

A gente já viu sobre derivadas aqui, mas vamos olhar o código do livro para seguir a linha de raciocínio.

Então suponha que a gente quer derivar uma equação de segundo grau

a+bx+cx^2

No R a gente declara ela como uma formula usando o ~(tio) e usa a função deriv

> y <- deriv(~a+b*x+c*x^2,"x") > y expression({ .value <- a + b * x + c * x^2 .grad <- array(0, c(length(.value), 1L), list(NULL, c("x"))) .grad[, "x"] <- b + c * (2 * x) attr(.value, "gradient") <- .grad .value })

A derivada é dada na linha .grad[, “x”] <- b + c * (2 * x). Para calcular o valor em algum gradiente, a gente precisa de um pouco mais de código. Em alguns casos a saída é dividida em expressões separadas. Por exemplo ao derivar a expressão [latex]e^{ax}+bx+cx^2[/latex] A sintaxe no R é a mesma usando a função deriv

> y <- deriv(~exp(a*x)+b*x+c*x^2,"x") > y expression({ .expr2 <- exp(a * x) .value <- .expr2 + b * x + c * x^2 .grad <- array(0, c(length(.value), 1L), list(NULL, c("x"))) .grad[, "x"] <- .expr2 * a + b + c * (2 * x) attr(.value, "gradient") <- .grad .value })

Mas aqui a derivada é composta, veja que a derivada

.grad[, "x"] <- .expr2 * a + b + c * (2 * x)

precisa da expr2 que está aqui

.expr2 <- exp(a * x)

O primeiro cenário:

Esse cenário ilustra a análise no caso em que não existe uma estrutura de idade e o fitness é simplesmente uma função quadrada da caraterística sendo estudada. Por causa da sua simplicidade, a análise e visualização são facilmente completadas. Apesar da simplicidade, é valido começar com algo simples e ir aprimorando, onde podemos adicionar complexidade aos poucos. Essa abordagem permite visualizar a importância de algumas premissas em particular, por exemplo aqui consideramos o caso de organismos semelparos, que a característica sendo analisada varia positivamente com um componente do fitness e negativamente com o outro componente. A priori, podemos ser tentados a assumir que isso vai necessariamente resultar em uma função que tem um máximo em algum valor intermediário da característica. Mas isso não é necessariamente verdade, suponha por exemplo duas funções de fitness

y_1=e^{ax}
y_2=e^{-bx}

Onde y_1 e y_2 são características como fecundidade e sobrevivência, a e b são constantes e x é a característica sendo estudada, nesse caso o tamanho do corpo. Vamos supor que o fitness seja o produto de y_1 e y_2, então o fitness W é dado por

W = y_1y_2=e^{(a-b)x}

Que não tem um valor ótimo intermediário. No primeiro cenário, um valor ótimo não é garantido, mas ocorre para um particular set de valores de parâmetros.

Pressupostos gerais:

1. O organismo é semelparo.
2. Fecundidade, F, aumenta com o tamanho do corpo x.
3. Sobrevivência, S, diminui com o tamanho do corpo x.
4. O fitness, W, é uma função da fecundidade e da sobrevivência.

Os pressupostos acima descrevem uma situação bem geral. Para fazer uma previsão, nos precisamos converter esses pressupostos em expressões matemáticas. Mas mesmo antes nos podemos fazer uma afirmação geral. Como um componente do fitness, fecundidade, aumenta com o tamanho do corpo, enquanto o outro componente, sobrevivência, diminui com o tamanho do corpo, em muitas circunstancias, mas não todas, nos esperamos que um tamanho ótimo do corpo maximize o fitness. Para investigar isso, nos vamos converter esses pressupostos em expressões matemáticas.

Pressupostos matemáticos

1. Fecundação aumenta linearmente com o tamanho do corpo:

F=a_F + b_F x

onde a_F e b_F são constantes.

2. Sobrevivência diminui linearmente com o tamanho do corpo:

S= a_S - b_S x

3. O fitness, W, é a expectativa reprodutiva ao longo da vida, R_0, dado como o produto da fecundidade e da sobrevivência.

W = R_0 = FS
W = (a_F + b_F x)(a_S - b_S x)
W = a_F a_S -b_F b_S x^2 + (a_S b_F - a_F b_S)x

A equação acima descreve uma parábola que é concava para baixo, isso é, tem um valor máximo em algum tamanho de corpo, digamos x^*. Apesar de sabermos de antemão que existe um valor ótimo de tamanho de corpo, nos não sabemos se ele ocorre em um valor plausível de tamanho de corpo. O primeiro passo na análise é plotar W com x para confirmar visualmente que o máximo existe e mostrar que ele ocorre para a espécie em estudo.

Plotando a função de fitness

Normalmente é uma boa ideia plotar a função para examinar visualmente o seu comportamento. Enquanto nos sabemos que neste caso a função é quadrática e concava para baixo, nos não sabemos se o valor ótimo para a característica avaliada é plausível biologicammente dados os parâmetros usados. Assim nos definimos valores dos parâmetros, vamos supor que valores razoáveis são a_F = 0, b_F=4, a_S = 1 e b_S = 0.5.
A equação do fitness pode então ser escrita como:

W = 0 \cdot 1 - 4 \cdot 0.5 x^2 + (1 \cdot 4 - 0 \cdot 0.5)x
W = -2x^2 + 4x

E no gráfico:

01

A partir da figura nos podemos ver que é um máximo em torno de 1.


Achando o máximo usando calculo

Para obter x^* usando calculo nos diferenciamos a equação, igualamos a zero e achamos o valor de x que satisfaz a condição.
Nesse caso a diferenciação é muito simples

\frac{dW}{dx} = -4x+4

\frac{dW}{dx} = 0 quando -4x+4 =0

logo

x=1

Diferenciação simbólica no R

Diferenciação simbólica, ou também chamada de diferenciação algorítmica (porque tem algorítimo que pensa nisso pra você) no R pode ser feita usando deriv

> y <- deriv(~-2*x^2+4*x,"x") > y expression({ .value <- -2 * x^2 + 4 * x .grad <- array(0, c(length(.value), 1L), list(NULL, c("x"))) .grad[, "x"] <- 4 - 2 * (2 * x) attr(.value, "gradient") <- .grad .value })

O output é meio “bagunçado”, mas isso porque ele não é pensado como um sumário, mas como uma entrada funcional para outras funções.
Com essa saída, a gente pode determinar em que valor a derivativa é zero usando a função uniroot. Nos setamos a derivativa como uma função separada e chamada pela função uniroot.

1
2
3
4
5
6
7
8
# 
FUNC <- function(w) {
    y <- deriv(~-2*x^2+4*x,"x")
    x <- w
    z <- eval(y)
    d <- attr(z,"gradient")
    return(d)
}
> B <- uniroot(FUNC, interval= c(-2,4)) > B$root [1] 1

Achando o máximo usando uma aproximação por métodos númericos.

Em muitos casos, a função pode não ser assim tão facilmente diferenciada, por exemplo, a função pode consistir de uma modelo de simulação ou ela pode ter discontinuidades, nos vamos encontrar comumente, mas esse cenário pode servir de exemplo para ilustrar como é o procedimento.
As rotinas disponíveis tipicamente localizão o valor mínimo da função. No nosso caso queremos o máximo. Para usar rotinas de minimização nos podemos simplismente pegar o valor negativo da equação. Logo, ao invés de procurar o máximo de 4x-2x^2, nos procuramos o mínimo de -4x+2x^2. Existe muitas rotinas prontas no R para achar o mínimo, aqui a gente vai usar o nlm (uma altenativa ao optimize).

1
2
3
FITNESS <- function(x) {
    return(2*x^2-4*x)
}
> nlm(FITNESS, p=-2) $minimum [1] -2 $estimate [1] 0.9999995 $gradient [1] 2.220446e-10 $code [1] 1 $iterations [1] 1 > nlm(FITNESS, p=-2) $minimum [1] -2 $estimate [1] 0.9999995 $gradient [1] 2.220446e-10 $code [1] 1 $iterations [1] 1

Agora que temos o resultado, é só sair medindo indivíduos dessa espécie na natureza e ver se eles medem algo em torno de 1 para saber que a teoria está correta, senão alterar os pressupostos e começar denovo. Mas o legal é que temos agora um resultado esperado a priori. Tudo bem que faltou definir unidades direito, precisamos pensar também como são os parâmetros da espécie em questão, mas da pra começar a pensar na metodologia para tentar aplicar. Gostaria de ter tido uma cabeça melhor para ter estudado esse tipo de coisa na graduação. Mas Modeling Evolution – an introduction to numerical methods é um livro muito legal, ele nem existia antes mesmo, então não tinha como ler, mas conforme for lendo mais, vou ir postando aqui porque se eu não escrevo, eu não aprendo.

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:
Derek A. Roff 2010 – Modeling Evolution – an introduction to numerical methods. Oxford University Press 352pp.

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
y <- deriv(~a+b*x+c*x^2,"x")
y
 
#
y <- deriv(~exp(a*x)+b*x+c*x^2,"x")
y
 
#
x <- seq(0,2,length=1000)
W <- (-2*x^2 + 4*x)
jpeg("01.jpg")
plot(x,W,type="l",xlab="Tamanho do corpo, x",ylab="Fitness, W",las=1,lwd=3,
     frame=F)
dev.off()
 
 
y <- deriv(~-2*x^2+4*x,"x")
y
 
# 
FUNC <- function(w) {
    y <- deriv(~-2*x^2+4*x,"x")
    x <- w
    z <- eval(y)
    d <- attr(z,"gradient")
    return(d)
}
 
#
B <- uniroot(FUNC, interval= c(-2,4))
B$root
 
 
FITNESS <- function(x) {
    return(2*x^2-4*x)
}
 
nlm(FITNESS, p=-2)