* [[:cursos:ecor:02_tutoriais:tutorial7b:start|Tutorial]]
* [[:cursos:ecor:01_curso_atual:exercicios8| Exercícios]]
====== 7b. Modelos Lineares Múltiplos ======
====== Modelos Lineares Múltiplos ======
===== Simplificando Modelos =====
Durante o curso usaremos o procedimento de simplificar o modelo a partir do modelo cheio.
O procedimento consiste em comparar modelos aninhados, dois a dois, retendo o que está mais acoplado aos dados. Caso os modelos não seja diferentes no seu poder explicativo, retemos o modelo mais simples, apoiados no princípio da parcimônia.
==== Princípio da parcimônia (Navalha de Occam) ====
* número de parâmetros menor possível
* linear é melhor que não-linear
* reter menos pressupostos
* simplificar ao mínimo adequado
* explicações mais simples são preferíveis
==== Método do modelo cheio ao mínimo adequado ====
- ajuste o modelo máximo (cheio)
- simplifique o modelo:
* inspecione os coeficientes (summary)
* remova termos não significativos
- ordem de remoção de termos:
* interação não significativos (maior ordem)
* termos quadráticos ou não lineares
* variáveis explicativas não significativas
* agrupe níveis de fatores sem diferença
* ANCOVA: intercepto não significativos -> 0
==== Tomada de decisão ====
** A diferença não é significativa: **
* retenha o modelo mais simples
* continue simplificando
**A difereça é significativa: **
* retenha o modelo complexo
* este é o modelo __MINÍMO ADEQUADO__
===== Interação entre preditoras =====
A interação é um elemento muito importante quando temos mais de uma preditora, pois desconsiderá-la pode limitar o entendimento dos processos envolvidos. Um exemplo cotidiano da interação é visto no uso de medicamentos e o alerta da bula sobre interação medicamentosa ou efeitos colaterais para pessoas portadoras de doenças crônicas. Dizemos que um medicamento tem interação com outra substância quando o seu efeito é modificado pela presença de outra substância, como por exemplo a ingestão de álcool junto com muitos medicamentos. Nos modelos, a interação tem uma interpretação similar, a resposta pelo efeito de uma variável preditora se altera com a presença de outra preditora.
==== Simulando um experimento plausível ====
Vimos que existe um efeito do tipo de solo na produção de um cultivar no exemplo de ANOVA. Uma expectativa plausível é que a adição de adubo também tenha efeito na produtividade e modifique o efeito do solo. Esse é nosso próximo exemplo. Para ele vamos usar uma simulação de dados similar ao que fizemos no modelo linear simples.
Nos dados originais do exercício de ANOVA a produtividade média nos solos foi de:
* arenoso: 9.9
* argiloso: 11.5
* humico: 14.3
Vamos, a partir dessa informação, criar um experimento onde, além da diferença do solo, metade dos cultivos foram tratados com adubo orgânico.
*1. Criamos vetores para representar as variáveis solo e adubo.
*2. Para cada observação incluímos o efeito médio de produtividade de cada solo (10 réplicas para cada solo)
*3. Associamos um valor de efeito do tratamento adubo, como:
* arenoso: + 2.7
* argiloso: + 0.7
* humico: + 0.2
*4. Em seguida somamos um efeito aleatório na resposta para criar um data frame com as variáveis preditoras e resposta.
set.seed(42)
solo <- rep(c("are", "arg", "hum"), each=20)
adubo <- as.factor(rep(rep(c(0, 1), each=10), 3))
meansolo <- rep(c(9.9, 11.5, 14.3), each=20)
efeitoadubo <- rep(c(0, 2.4, 0, 0.9, 0, 0.2), each=10)
residuo <- rnorm(60, 0, 0.5)
dadosolo <- data.frame(solo, adubo, prod = meansolo + efeitoadubo + residuo)
str(dadosolo)
Confira se o objeto ''dadosolo'' foi organizado corretamente
==== Gráfico dos dados ====
Agora um gráfico simples. Busque entender todos os argumentos das funções abaixo.
par( mar=c(4,4,2,2), cex.lab=1.5, cex.axis=1.2, las=1, bty="n")
boxplot(prod ~ adubo + solo, data = dadosolo, ann= FALSE, xaxt= "n", outline= FALSE, col = rep(c(rgb(0,0,0, 0.1),rgb(0,0,0, 0.5)), 3) )
mtext(c("arenoso", "argiloso", "húmico"), side = 1, at= c(1.5, 3.5, 5.5), line = 1, cex = 2)
legend("bottomright", legend= c("sem", "com"),title = "Adubo", bty= "n", pch = 15, cex = 1.5,col = c(rgb(0,0,0, 0.1),rgb(0,0,0, 0.5)))
==== Modelo Cheio ====
Abaixo construímos o modelo cheio com as variáveis adubo e tipo de solo.
soloFull <- lm(prod ~ adubo + solo + solo:adubo, data = dadosolo)
summary(soloFull)
==== Modelo sem Interação ====
A primeira simplificação possível é retirar o efeito da interação entre as preditoras e comparar com o modelo cheio.
solo01 <- lm(prod ~ adubo + solo , data = dadosolo)
anova(solo01, soloFull)
O resultado nos indica que o modelo cheio é o modelo mínimo adequado. Ou seja, explica uma porção considerável da variação dos dados a mais que o modelo mais simples, sem a interação entre tipo de solo e adubo. Para completar, vamos fazer a comparação com o modelo nulo. Essa comparação pode ser feito de duas maneiras: (1) construindo o modelo nulo e comparando por anova, ou (2) interpretando a tabela de anova do modelo mínimo adequado.
solo00 <- lm(prod ~ 1 , data = dadosolo)
anova(solo00, soloFull)
anova(soloFull)
==== Diagnóstico ====
O passo final é investigar se o modelo cumpre com as premissas do modelo linear.
par(mfrow = c(2,2), mar=c(4,4,2,2), cex.lab=1.2,
cex.axis=1.2, las=1, bty="n")
plot(soloFull)
Não poderia ser mais comportado. Isso significa que criamos os dados corretamente!!
Agora é a parte mais difícil e interessante de qualquer análise de dados, a interpretação biológica suscita do resultado!
**__Interpretando Variáveis Indicadoras (Dummy)__**
As variáveis indicadoras devem ser interpretadas com cuidado. No exemplo acima, o modelo pode ser descrito da seguinte forma:
$$ y_{tr} = \alpha + \beta_1 * arg + \beta_2 * hum + \beta_3 * adubo + \beta_4 * arg * adubo + \beta_5 * hum * adubo $$
As variáveis __arg__, __hum__ e __adubo__ são dummy ou indicadoras, representadas por 1 quando presente e 0 quando ausentes. $\alpha, \beta_i$ representam as estimativas do modelo e estão relacionados, nesse caso, ao efeito de cada tratamento.
Para calcular o valor predito para o tratamento no solo arenoso com adubo, temos:
$$ y_{arenAdubo} = \alpha + \beta_3 * adubo $$
Isso em decorrência do tratamento **arenoso sem adubo** estar representado pelo intercepto ($\alpha$) do modelo.
Para o tratamento de solo **argiloso com adubo** o predito é:
$$ y_{argAdubo} = \alpha + \beta_1 * arg + \beta_3 * adubo + \beta_4 * arg * adubo $$
E assim por diante, usando as variáveis indicadoras e os coeficientes estimados para o cálculo do predito pelo modelo.
===== Peso de bebês ao nascer =====
Vamos analisar o dado de peso dos bebês ao nascer e como isso se relaciona às características da mãe. Esses dados pode ser consultados em [[https://www.stat.berkeley.edu/users/statlabs/labs.html]].
* baixe o arquivo {{::cursos:ecor:02_tutoriais:tutorial7c:babies.csv|}} no seu diretório de trabalho
* Vamos selecionar o modelo mínimo adequado a partir das variáveis:
* resposta **bwt** : peso do bebê ao nascer em onças(oz)
* preditoras:
* gestation: tempo de gestação (dias)
* age: idade
* weight: peso da mãe
* smoke: 0 não fumante; 1 fumante
Para simplificar nosso tutorial vamos usar apenas as preditoras: tempo de gestação, idade da mãe e se ela é fumante ou não ((no exercício terão que usar os dados brutos e todas as variáveis)).
bebes <- read.table("babies.csv", header= TRUE, as.is = TRUE, sep= "\t")
str(bebes)
mlfull <- lm(bwt ~ gestation + age + smoke
+ gestation:age + gestation:smoke
+ age: smoke + gestation:age:smoke, data = bebes)
summary(mlfull)
===== Interação Tripla =====
Vamos simplificar o modelo, retirando a interação ''gestation:age:smoke'' que aparenta não ser importante.
ml01 <- lm(bwt ~ gestation + age + smoke
+ gestation:age + gestation:smoke
+ age: smoke, data = bebes)
anova(ml01, mlfull)
summary(ml01)
===== Interações Dupla =====
Continuamos a simplificação, retirando as interações duplas uma a uma para avaliar quais delas devem ser mantidas. Os testes parciais das variáveis no ''summary'' nos dá uma indicação de quais devem ser mantidas, mas uma boa prática é fazer o processo completo, já que um elemento no modelo pode mudar o efetividade de outro, principalmente quando compartilham alguma porção de variação explicada.
## sem age:smoke
ml02 <- lm(bwt ~ gestation + age + smoke
+ gestation:age + gestation:smoke, data = bebes)
anova(ml01, ml02)
## sem gestation:smoke
ml03 <- lm(bwt ~ gestation + age + smoke
+ gestation:age + age:smoke, data = bebes)
anova(ml01, ml03)
## sem gestation:age
ml04 <- lm(bwt ~ gestation + age + smoke
+ gestation:smoke + age: smoke, data = bebes)
anova(ml01, ml04)
A única interação dupla que não parece fazer diferença quando retiramos do modelo é a ''age:smoke'', as outras explicam uma porção razoável da variação dos dados.
===== Interpretação do modelo =====
O ''summary'' nos fornece as principais informações sobre o modelo mínimo adequado.
summary(ml02)
confint(ml02)
anova(ml02)
===== Diagnóstico do modelo =====
par(mfrow = c(2,2), mar=c(4,4,2,2), cex.lab=1.2,
cex.axis=1.2, las=1, bty="n")
plot(ml02)
Os dados desse estudo serão usados também no exercício, porém lá, vamos partir dos dados brutos com mais variáveis
/*
===== Regressão Múltipla passo a passo =====
Neste tutorial vamos usar o conjunto de dados //Pollute//, que tem dados do nível de poluição em algumas cidades e atributos destas cidades que podem servir como variáveis preditoras. Baixe o arquivo com os dados [[http://www.bio.ic.ac.uk/research/mjcraw/therbook/data/Pollute.txt|daqui]], e crie um objeto de dados no R com ele:
poluicao <- read.table("Pollute.txt", header=T, sep="\t")
==== Regressão Múltipla e Seleção de Modelos ====
Vamos usar como preditoras as variáveis número de indústrias, a temperatura média e velocidade média do vento. Inspecione a relação da variável-resposta com cada preditora:
par(mfrow=c(1,3))
plot(Pollution~Industry, data=poluicao)
plot(Pollution~Temp, data=poluicao)
plot(Pollution~Wind, data=poluicao)
par(mfrow=c(1,1))
Nosso primeiro modelo de regressão linear terá apenas o efeito do número de indústrias:
pol.m1 <- lm(Pollution~Industry, data=poluicao)
Um gráfico da linha de valores esperados para este modelo sobre os dados:
plot(Pollution~Industry, data=poluicao)
abline(pol.m1, col="blue")
Há um ponto extremo, correspondente à cidade mais industrializada, que pode ter uma alavancagem forte. Vamos verificar acrescentando a reta de um modelo ajustado sem este ponto:
abline(lm(Pollution~Industry, data=poluicao, subset=Industry
A alavancagem não é forte. Vamos seguir em frente, com um teste do modelo
anova(pol.m1)
Este teste é na verdade uma comparação com o modelo mais simples possível, que é ausência de efeito de qualquer preditora. Neste caso, o valor previsto para a variável resposta é simplesmente sua média:
pol.m0 <- lm(Pollution~1, data=poluicao)
##compare:
anova(pol.m0,pol.m1)
anova(pol.m1)
Já entendeu? Veja o resumo do modelo mais simples:
summary(pol.m0)
E compare com a média e o desvio-padrão da variável-resposta:
mean(poluicao$Pollution)
sd(poluicao$Pollution)
Um gráfico com as linhas dos dois modelos concorrentes:
plot(Pollution~Industry, data=poluicao)
abline(pol.m1, col="blue")
abline(h=mean(poluicao$Pollution),col="red")
Agora vamos complicar o modelo, com o acréscimo de mais variáveis preditoras. Acrescentamos o efeito de temperatura:
pol.m2 <- lm(Pollution~Industry+Temp, data=poluicao)
Comparação dos dois modelos:
anova(pol.m1,pol.m2)
Há um ganho importante de variação explicada com a inclusão da temperatura, para o qual o teste F é significativo.
Prosseguindo, vamos adicionar o efeito do vento, e comparar com o modelo anterior:
pol.m3 <- update(pol.m2,.~.+Wind)
anova(pol.m2,pol.m3)
anova(pol.m3)
Não houve ganho importante de explicação. Ficamos com o modelo anterior. Vamos inspecionar os gráficos de diagnóstico:
par(mfrow=c(2,2))
plot(pol.m2)
par(mfrow=c(1,1))
Há um ponto discrepante (41). vamos identificá-lo nos dados e nos gráficos:
poluicao[41,c(1:3,5)]
par(mfrow=c(1,2))
plot(Pollution~Industry, data=poluicao)
points(Pollution~Industry, data=poluicao[41,], col="red", pch=19)
plot(Pollution~Temp, data=poluicao)
points(Pollution~Temp, data=poluicao[41,], col="red", pch=19)
par(mfrow=c(1,1))
Apesar deste ponto, o modelo parece razoável. Vamos continuar na próxima seção investigando interações entre as preditoras. Mas antes é muito importante que você entenda a diferença entre os testes de cada fator no resumo do modelo obtido com
summary(pol.m2)
E o resultado do comando:
anova(pol.m2)
O comando ''anova'' sempre compara modelos. Se o argumento é um único modelo, ele retorna o resultado de um **teste sequencial**, isto é, o aumento de variação explicada pelo acréscimo de uma preditora em relação a um modelo com as preditoras **acima dele na tabela** A sequência das variáveis preditoras nesta tabela é a ordem delas na fórmula do modelo. Por isso, os dois comandos abaixo não retornam exatamente os mesmos resultados, compare os valores das somas dos quadrados:
anova(lm(Pollution~Temp+Wind+Industry, data=poluicao))
anova(lm(Pollution~Industry+Wind+Temp, data=poluicao))
Já o teste t de cada preditora no resumo obtido com a função ''summary'' é um **teste marginal**, pois testa a hipótese de que o efeito seja zero, **descontados os efeitos das demais preditoras**. Para a variável //Industry// este teste **equivale**((não é o mesmo teste)) ao seguinte comando, com a função ''anova'':
anova(lm(Pollution~Temp+Industry, data=poluicao),lm(Pollution~Temp, data=poluicao))
E para a variável Temp corresponde a:
anova(lm(Pollution~Temp+Industry, data=poluicao),lm(Pollution~Industry, data=poluicao))
Note que neste caso o teste é sempre uma comparação de um modelo com as demais variáveis exceto a que se pretende testar, com outro modelo que tem todas as preditoras, e também a que se pretende testar. É essencial que você compreenda isto para interpretar corretamente os resultados de qualquer modelo linear.
==== Interação ====
Começamos com uma análise gráfica, para descobrir sinais de interação entre as duas variáveis preditoras incluídas no modelo.
Primeiro definimos 4 faixas de temperatura parcialmente sobrepostas com a função ''equal.count'', do pacote //lattice//:
library(lattice)
temp.shingle <- equal.count(poluicao$Temp,number=4)
summary(temp.shingle)
E agora condicionamos o gráfico por esta variável categórica. O argumento ''panel'' é uma função((se você ainda não chegou à unidade sobre funções, não se preocupe, apenas analise o gráfico resultante)) a ser aplicada em cada painel, no caso plotar os pontos
e fazer uma regressão apenas com estes pontos (detalhes na ajuda do ''xyplot'')
xyplot(Pollution~Industry|temp.shingle,data=poluicao,
panel=function(x,y,...){
panel.xyplot(x,y,...)
panel.abline(lm(y~x),...)
}
)
Fazemos o mesmo tipo de gráfico para a variável //Wind//. Embora não tenha sido incluída no modelo, sua interação com outras preditoras pode ocorrer:
xyplot(Pollution~Industry|equal.count(poluicao$Wind,number=4),data=poluicao,
panel=function(x,y,...){
panel.xyplot(x,y,...)
panel.abline(lm(y~x),...)
}
)
Os dois gráficos mostram que, para algumas faixas de temperatura e de vento as inclinações das retas de regressão são muito diferentes, o que sugere interações entre as variáveis climáticas e o número de indústrias. Vamos criar modelos com estas interações, e avaliar se há um ganho com isto:
pol.m4 <- update(pol.m2,.~.+Industry:Temp)## acrescenta apenas a interação
pol.m5 <- update(pol.m2,.~.+Wind+Industry:Wind)##acrescenta a interação e seus componentes
anova(pol.m2,pol.m4)
anova(pol.m2,pol.m5)
Nenhuma das interações melhora o modelo significativamente. Continuamos com o modelo anterior.
==== Um Gráfico do Modelo ====
Comparar os os dados com os valores previstos pelo modelo que selecionamos não é fácil, pois agora temos duas variáveis preditoras, o que define uma **superfície de resposta**, e não mais uma linha. Uma possibilidade é dividir os dados em conjuntos definidos por uma das variáveis, e investigar a relação com a outra variável, como nos gráficos abaixo.
Primeiro criamos um fator por faixa da variável //Temp//, usando seus intervalos de tercis (quantis a 1/3, 2/3 e 3/3) para separar as observações em terços:
##Tercis da variavel Temp
Temp.tercis <- quantile(poluicao$Temp,c(0,1/3,2/3,1))
Temp.tercis
##Criando um fator com estes breakpoints, com a funcao cut
poluicao$Temp.c <- cut(poluicao$Temp, breaks=Temp.tercis)
## Verificando: três fatores ordenados de faixas de temperatura
levels(poluicao$Temp.c)
Em seguida fazemos três graficos com os valores selecionados po faixa de temperatura
cf.m2 <- coef(pol.m2)
par(mfrow=c(1,3))
##Primeiro Grafico: primeiro nivel do fator faixa de temperatura
plot(Pollution~Industry, data=poluicao, subset=Temp.c==levels(poluicao$Temp.c)[1],
main=paste("Temp: ",Temp.tercis[1]," a ",Temp.tercis[2],sep=""),
xlim=c(min(poluicao$Industry),max(poluicao$Industry)),
ylim=c(min(poluicao$Pollution),max(poluicao$Pollution)) )
##Linhas do previsto para ponto medio da faixa de temperatura
abline(cf.m2[1]+cf.m2[3]*mean(Temp.tercis[1:2]),cf.m2[2], col="blue", lwd=2)
##Segundo Grafico
plot(Pollution~Industry, data=poluicao, subset=Temp.c==levels(poluicao$Temp.c)[2],
main=paste("Temp: ",round(Temp.tercis[2],1)," a ",round(Temp.tercis[3],1),sep=""),
xlim=c(min(poluicao$Industry),max(poluicao$Industry)),
ylim=c(min(poluicao$Pollution),max(poluicao$Pollution)) )
abline(cf.m2[1]+cf.m2[3]*mean(Temp.tercis[2:3]),cf.m2[2], col="blue", lwd=2)
##Terceiro grafico
plot(Pollution~Industry, data=poluicao, subset=Temp.c==levels(poluicao$Temp.c)[3],
main=paste("Temp: ",round(Temp.tercis[3],1)," a ",round(Temp.tercis[4],1),sep=""),
xlim=c(min(poluicao$Industry),max(poluicao$Industry)),
ylim=c(min(poluicao$Pollution),max(poluicao$Pollution)) )
abline(cf.m2[1]+cf.m2[3]*mean(Temp.tercis[3:4]),cf.m2[2], col="blue", lwd=2)
par(mfrow=c(1,1))
*/
===== O Modelo por trás da ANOVA =====
Neste tutorial usaremos dados de um [[http://en.wikipedia.org/wiki/Randomized_block_design|experimento em blocos]] para comparar o crescimento de mudas de diferentes espécies de árvores em diferentes substratos. Baixe o conjunto de dados [[dados:dados-mudas|aqui]].
Neste tutorial vamos desconsiderar os blocos, e usar uma análise de variância de dois fatores para testar diferenças entre espécies das mudas e entre os substratos usados.
Veremos que por trás da ANOVA há um modelo linear, que estabelece uma relação entre uma resposta contínua (crescimento) e variáveis preditoras categóricas, que são os fatores (espécie de planta e substrato usado).
Comece com a leitura dos dados e conversão da variaveis de bloco((mesmo não usando, é sempre bom converter fatores representados por números em categorias))e substrato, que são numericas, em fatores:
mudas <- read.csv("altura-mudas.csv")
mudas$bloco <- as.factor(mudas$bloco)
mudas$substrato <- as.factor(mudas$substrato)
Ajustamos um primeiro modelo, no qual a altura das mudas depende da especie
mudas.m1 <- lm(altura~especie, data=mudas)
Faça uma avaliação deste modelo, e inspecione seus coeficientes:
summary(mudas.m1)
## Coeficientes do modelo
cf.mud1 <- coef(mudas.m1)
cf.mud1
O que significam estes coeficientes? A resposta está na matriz do modelo:
head(model.matrix(mudas.m1))
Vamos usá-la para o cálculo dos valores esperados pelo modelo. Como são muitos valores, inspecione o começo e o fim do vetor resultante com as funções ''head'' e ''tail'':
head(model.matrix(mudas.m1)%*%cf.mud1)
tail(model.matrix(mudas.m1)%*%cf.mud1)
Já entendeu o que são os coeficientes? Isto pode ajudar: veja as médias de alturas por espécies:
tapply(mudas$altura,mudas$especie,mean)
E compare com os coeficientes no resumo do modelo:
summary(mudas.m1)
Prossiga apenas após estar certo(a) do significado dos coeficientes neste primeiro modelo
Agora vamos criar um novo modelo, que terá também o efeito do substrato:
mudas.m2 <- update(mudas.m1,.~.+substrato)
Inspecione os coeficientes:
cf.mud2 <- coef(mudas.m2)
cf.mud2
E a matriz do modelo:
head(model.matrix(mudas.m2))
Repita a multiplicação da matriz do modelo pelos coeficientes e entenda como isto gera os valores previstos.
Avalie a adição do fator substrato ao modelo com o comando:
anova(mudas.m2)
**MORAL DA HISTÓRIA**: uma análise de variância é o teste de significância marginal para um modelo linear com variáveis categóricas!