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

O que acontece se o modelo ou formula da análise que você usa esta errado? O exemplo do paradoxo de Simpson.

Essa é uma pergunta que eu sempre pensei, o que acontece quando a gente usa a análise errada? Lembrando que Anova, Regressão, geral ou não, linear ou não. Você esta informando não so os seus dados, mas também o que você acha que acontece, como o mundo funciona.
E o que acontece quando a gente tem a informação correta, mas informa a formula errada?

Um exemplo interessante é o paradoxo de Simpson. Vamos tentar inventar um exemplo aqui para entender melhor esse caso.

Vamos supor que estamos interessados em saber como a altitude (variável preditora) influencia no DAP (diâmetro a altura do peito) de uma espécie de árvore. Talvez alguém queria saber onde plantar árvores para produzir mais madeira. Vamos inventar alguns dados apenas para ilustrar a situação. Então temos uma montanha, coletamos em três lugares, base, lateral e topo da montanha, e para cada planta medimos o DPA e registramos com o GPS a sua altura dessa árvore. No final temos várias plantas com suas respectivas alturas e DAP.

01

Legal, parece que quanto maior a altitude, maior o DAP, aplicamos uma regressão linear e como resultado vemos:

> summary(modelo) Call: lm(formula = resposta ~ altitude, data = dados) Residuals: Min 1Q Median 3Q Max -2.4560 -0.7646 0.2089 0.7864 3.5954 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.4484 0.5035 0.890 0.38084 altitude 0.8883 0.3024 2.937 0.00656 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.375 on 28 degrees of freedom Multiple R-squared: 0.2355, Adjusted R-squared: 0.2082 F-statistic: 8.626 on 1 and 28 DF, p-value: 0.00656

Legal, vemos uma estimativa positiva para a altitude (0.8883), p significativo (0.00656 **) com muitas estrelinha. Ai analisamos os resíduos para ver se existe muito problemas no modelo…

02

E fora alguns outliers, parece que os resíduos tem uma distribuição normal, e não apresentam uma tendência. Talvez esse não seja o melhor exemplo, mas o problema aqui é ao simplificar a forma como olhamos os dados, não considerando onde coletamos os dados, somente aquilo que era o interesse, altura e DAP, há uma inversão no padrão que é o mais importante da nossa questão, a inclinação da reta, a relação entre altitude e DAP.
Vamos olhar novamente a figura, agora separando por cores cada população coletada, a da base, entorno e pico da montanha.

03

Agora vemos que para cada população a tendência é que quanto maior a altitude, menor o DAP.

> summary(modelo2) Call: lm(formula = resposta ~ altitude + pop, data = dados) Residuals: Min 1Q Median 3Q Max -1.6268 -0.6942 -0.1737 0.6092 1.8863 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.7965 0.4470 4.019 0.000445 *** altitude -2.2413 0.5928 -3.781 0.000826 *** poppop2 2.7382 0.6269 4.368 0.000178 *** poppop3 6.7670 1.1999 5.640 6.27e-06 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.9565 on 26 degrees of freedom Multiple R-squared: 0.6566, Adjusted R-squared: 0.617 F-statistic: 16.57 on 3 and 26 DF, p-value: 3.184e-06

O que acontece é que cada população tem um intercepto diferente, ele começa em um ponto diferente. Podemos pensar por exemplo, que com o aumento da altitude, espécies vão sumindo, e que isso deixa a espécie de nosso interesse em melhor situação, crescendo mais, mas o padrão é que quanto maior a altitude, mais o DPA diminui. Outra explicação é que altitude também altera os nutrientes disponíveis, o que pode alterar o DAP, mas ainda sim quanto maior a altitude, menor o DPA, essa relação altitude altitude vs DPA é negativa. O ponto é que na tentantiva de simplificar o mundo com menos parâmetros, invertemos a relação que observamos. A conclusão ficou totalmente diferente. Apesar que os dois casos tudo é significativo, a conclusão é totalmente diferente. Pense nos conselhos que daríamos para “onde plantar árvores” para cada caso.
Esse caso particular é conhecido como paradoxo de Simpson, que diz respeito a como um padrão se altera ao não considerarmos uma diferença entre interceptos. Eu tentei exemplificar um caso que acho mais comum a biologia, mas no artigo do Wikipédia tem excelentes exemplo.

Agora uma coisa que podemos fazer é olhar os likelihood. Nesse caso se olharmos os likelihoods, vemos que:

> logLik(modelo) 'log Lik.' -51.09479 (df=3) > logLik(modelo2) 'log Lik.' -39.0885 (df=5) > logLik(modelo)/logLik(modelo2) 'log Lik.' 1.307157 (df=3)

O “modelo” que não considera a diferença dos interceptos tem um loglikelihood maior (mais negativo) que o do “modelo2”. Podemos dividir um pelo outro e vemos que somos mais inclinados a acreditar no modelo2, se o resultado fosse 1, seria porque os likelehoods são iguais, então os modelos seriam igualmente verosimilhantes, teriam uma capacidade igual de previsão, mas desvios dos 1 dizem respeito a probabilidade de previsão melhor de um em relação ao outro. So pensar que sempre o que for valores mais da “ponta do gradiente” sempre se pareceram com outliers no “modelo” mais simples, se você não vê esses valores, cada vez vai ficar mais parecido que é inútil olhar a diferença entre as populações. Ou seja quanto menos a amostra conter todo o gradiente, pior vai ser para descobrir que caimos dentro desse paradoxo. Alias aqui se a amostra não tiver uma distribuição uniforme de altitudes, e ficar concentrado numa altitude, maior o perigo de ficar enroscado nesse tipo de paradoxo, qualidade da informação acaba sendo mais importante que quantidade.

Mas podemos ainda testar esse desvio da razão dos likelihood usando a distribuição F. Assim temos.

> anova(modelo,modelo2) Analysis of Variance Table Model 1: resposta ~ altitude Model 2: resposta ~ altitude + pop Res.Df RSS Df Sum of Sq F Pr(>F) 1 28 52.965 2 26 23.789 2 29.177 15.944 3.027e-05 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Uma diferença entre os modelos, ou seja eles são diferentes, então deveríamos considerar nesse caso aquele que tem o maior likelihood. (lembrando que o maior likelihood é o valor do modelo 2 -39.09).
Esse paradoxo eu achei legal principalmente para ter idéia de que nem sempre o modelo mais simples é o melhor. Apesar de tendermos sempre a simplificar o máximo possível os modelos do mundo seguinte a teoria da Navalha de Occan, como disse Einstein “A scientific theory should be as simple as possible, but no simpler.”

Uma teoria cientifica tem que ser o mais simples possível, mas não mais simples.

#simpson paradox
set.seed(51)
pop1<-runif(10,0,1)
pop2<-runif(10,1,2)
pop3<-runif(10,2,3)
 
resp1<-rnorm(10,2-2*pop1)
resp2<-rnorm(10,4-2*pop2)
resp3<-rnorm(10,8-2*pop3)
 
dados<-data.frame(resposta=c(resp1,resp2,resp3),altitude=c(pop1,pop2,pop3),
pop=rep(c("pop1","pop2","pop3"),each=10))
head(dados)
 
plot(resposta~altitude,data=dados,pch=19,xlab="Altitude",ylab="DAP",frame.plot=F)
abline(lm(dados$resposta~dados$altitude))
 
modelo<-lm(resposta~altitude,data=dados)
summary(modelo)
 
par(mfrow=c(2,2))
plot(modelo)
 
plot(resposta~altitude,data=dados,pch=19,col=as.numeric(dados$pop)+1,xlab="Altitude",ylab="DAP",frame.plot=F)
abline(modelo)
abline(lm(dados[dados$pop=="pop1","resposta"]~
dados[dados$pop=="pop1","altitude"]),lty=2,col=2)
abline(lm(dados[dados$pop=="pop2","resposta"]~
dados[dados$pop=="pop2","altitude"]),lty=2,col=3)
abline(lm(dados[dados$pop=="pop3","resposta"]~
dados[dados$pop=="pop3","altitude"]),lty=2,col=4)
 
modelo2<-lm(resposta~altitude+pop,data=dados)
summary(modelo2)
 
extractAIC(modelo)
extractAIC(modelo2)
 
logLik(modelo)
logLik(modelo2)
 
anova(modelo,modelo2)

Parâmetros são mais importantes que valores p

Parâmetros são mais importantes que valores p, ou pelo menos deviam ser observados com carinho.

Os valores p muitas vezes são cruéis. E vou ilustrar com uma análise que vi durante uma disciplina para a graduação em biologia, em que alunos fazem projetos rápidos como treinamento para pesquisa.

Mudanças na diversidade, riqueza ou composição de espécies são sempre um projeto batata, diferenças entre borda e centro de matas por exemplo sempre ocorrem. Essas relações entre um gradiente abiótico e mudanças na riqueza, diversidade ou composição de espécies acabam sendo comum como projetos de treinamento desse tipo.

Pois nesse caso o projeto era o seguinte, avaliar a relação da profundidade de um local de uma lagoa com a diversidade de plantas aquáticas ali. Salvo outras considerações, tudo era simples e direto. Coletar uma área, um quadrado de 1 metro quadrado, pegando todas as plantas, quantificando a abundância de cada espécie nesse quadrado e medir a profundidade da agua nesse local.

Tudo bem, coletas feitas, muitas amostras, dados foram tabulados e análises produzidas.

E la saiu uma figura muito bonita.

Feito a análise, tudo significativo e com muitas estrelinhas:

 
> summary(modelo) Call: lm(formula = diversidade ~ profundidade) Residuals: Min 1Q Median 3Q Max -0.018415 -0.006439 -0.001135 0.004694 0.017791 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.997717 0.002693 1113.041 < 2e-16 *** profundidade 0.020549 0.002307 8.907 9.72e-12 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.009179 on 48 degrees of freedom Multiple R-squared: 0.623, Adjusted R-squared: 0.6152 F-statistic: 79.34 on 1 and 48 DF, p-value: 9.717e-12

Tudo bem e todo mundo feliz. Mas havia um problema. A inclinação do modelo era um valor muito baixo. Na verdade eu não lembro os valores reais. Mas era algo como isso, de acordo com o modelo, precisaria de 5 metros a mais de profundidade pra aumentar uma espécie, ou mesmo alterações na equitabilidade que alterassem a diversidade, seria impossível com essa inclinação realmente observar algo dentro do limite físico de profundidade de lagoas que existem.

Esse gráfico refletiria melhor o que acontece, sempre tudo igual. Existe um limite físico para profundidade das lagoas e nunca veríamos uma mudança na diversidade nesse ritmo de alteração.
Mas como o gráfico anterior ficou bonito e “significativo”, pouca atenção foi dada aos parâmetros estimados. Mas o ponto é que não havia la grandes problemas na análise , e sim no fato que mesmo os parâmetros sendo significativos eles não faziam muito sentido biológico, não naquela magnitude de mudança.
Mas as pessoas acabam por odiar quem tira valores p significativos delas.

set.seed(13)
profundidade<-runif(50,0,2)
diversidade<-rnorm(50,3+0.02*profundidade,0.01)
 
#Figura 1
plot(diversidade~profundidade,xlab="Profundidade local",
ylab="Diversidade de biológica",pch=19)
abline(lm(diversidade~profundidade))
 
#Regressão linear
modelo<-lm(diversidade~profundidade)
summary(modelo)
 
#Figura 2
#preparando para desenhar o intervalo de confiança, fazendo predições
novosdados<-data.frame(profundidade=seq(0,5,by=0.01))
predição<-predict(modelo,newdata=novosdados,interval = c("prediction"))
 
plot(diversidade~profundidade,ylim=c(2.5,3.5),xlim=c(0,5),
frame.plot=F,type="n",xlab="Profundidade local",
ylab="Diversidade biológica")
 
polygon(x=c(novosdados[,1],rev(novosdados[,1])),
y=c(predição[,3],rev(predição[,2])),
col="gray90",border="gray90")
 
points(predição[,1]~novosdados[,1],type="l",lwd=2)
points(diversidade~profundidade,type="p",pch=19)
 
legend("topleft",legend=c("Modelo","Intervalo de predição"),lty=c(1,1),
lwd=c(2,20),col=c("black","gray90"))