7b. Modelos Lineares Múltiplos [ecoR]

Essa página está em modo somente de leitura. Você pode visualizar a fonte, mas não alterá-la. Informe-se com o administrador do Wiki, caso você ache que isso está incorreto.
<WRAP tabs>
  * [[:02_tutoriais:tutorial7b:start|Tutorial]]
  * [[:01_curso_atual:exercicios8| Exercícios]] 
 </WRAP>

====== 7b. Modelos Lineares Múltiplos ======

<WRAP center round box 80%>
Videoaula do curso **Princípios de Planejamento e Análise de Dados**. Os conceitos abordados são os mesmos desse tutorial, desconsidere referências à disciplina. 
<WRAP center round tip 80%>
{{youtube>89hcuYp9oWo}}

</WRAP>
</WRAP>

Os modelos lineares permitem que sejam incluídas mais do que uma variável preditora, como fizemos até aqui. Nesse tutorial vamos aprender alguns princípios básicos desses modelos mais complexos. Em modelos com mais de uma variável preditora, precisamos tomar a decisão de quais variáveis devemos reter em nosso modelo. É desejável interpretar o modelo mais simples e que contém apenas as variáveis que explicam porções consideráveis da variação na variável resposta.

<WRAP center round box 90%>
**__Formulas Estatísticas__**

<WRAP center round tip 90%>
O argumento ''formula'' da função ''lm'' funciona de forma diferente das formulas matemáticas e deve-se ter cuidado com a inclusão de termos ou operações dentro dela.

Alguns aspectos básicos do argumento:

  * ''y ~ x'' indica: construa o modelo da variável resposta ''y'' como **função estatística linear** de ''x'';
  * ''y ~ x1 + x2'' indica: construa o modelo estatístico de ''y'' como **função linear** das variáveis ''x1'' e ''x2'' como tendo efeitos aditivos;

Se quisermos utilizar os símbolos matemáticos no sentido matemático usual **dentro** de uma fórmula estatística, temos que utilizara a função ''I()'' para que a operação seja realizada antes da construção do modelo:

  * ''y ~ I( x1^2 * x2^3 )'' indica: modele y como função estatística **da variável** ''x1^2 * x2^3'';
  * 
  * '' y ~ I( x1 / x2 )'' indica: modele y como função estatística **da variável** ''x1/x2'';

No caso de utilizarmos **funções matemáticas** específicas a função ''I()'' torna-se desnecessária:
  * ''log(y) ~ log(x)'' indica: modele o ''log(y)'' com função estatística da variável ''log(x))'';
  * ''log(y) ~ log(x1^2 * x2)'' indica: modele o **log(y)** com função estatística da variável ''log(x1^2 * x2)'';


</WRAP>


</WRAP>


===== Modelos Plausíveis =====

Quando temos uma hipótese onde há mais de uma variável preditora, precisamos avaliar, antes de iniciar as análises, quais modelos são plausíveis e relacionados a que hipótese alternativa. Os modelos só são construídos a partir dessa avaliação. Vamos usar o exemplo da videoaula para exemplificarmos os procedimentos e conceitos relacionados a esse tutorial. 


===== Interação entre preditoras =====

<WRAP center round box 80%>
Videoaula do curso **Princípios de Planejamento e Análise de Dados**. Os conceitos abordados são os mesmos desse tutorial, desconsidere referências à disciplina. 
<WRAP center round tip 80%>
{{youtube>Mx9skekN6e8}}

</WRAP>
</WRAP>
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.  

<WRAP center round box 90%>

**__Delineamentos Experimentais__**

<WRAP center round tip 100%>
^ Expressão ^ Significado ^
| ''y ~ x'' | Modele ''y'' como função estatística de ''x'' |
| ''y ~ x1 + x2'' | inclua as variáveis ''x1'' e ''x2'' como preditoras|
| ''y ~ x1 + x2 + x1:x2'' | inclua também a interação de ''x1'' com ''x2''|
| ''y ~ x1 * x2'' | mesmo que ''y ~ x1 + x2 + x1:x2''  |
| ''y ~ (x1 + x2 + x3)^2'' | Adiciona acima '' + x3 + x1:x3 + x2:x3''| 
| ''y ~ (x1 + x2 + x3)^3'' | Adiciona acima '' + x1:x2:x3''|
| ''y ~ (x1 + x2 + x3)^3 - x1:x2'' | Retira o termo ''x1:x2'' da fórmula acima|

</WRAP>

</WRAP>

/*
==== 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.

<WRAP center round box 90%>
  *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.

</WRAP>

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

</code>

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.

<code rsplus>

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

</code>

==== Modelo Cheio ====

Abaixo construímos o modelo cheio com as variáveis adubo e tipo de solo.

<code rsplus>

soloFull <- lm(prod ~ adubo + solo + solo:adubo, data = dadosolo)
summary(soloFull)

</code>


==== 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.  

<code rsplus>

solo01 <- lm(prod ~ adubo + solo , data = dadosolo)
anova(solo01, soloFull)

</code>


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.


<code rsplus>

solo00 <- lm(prod ~ 1 , data = dadosolo)
anova(solo00, soloFull)
anova(soloFull)

</code>

==== Diagnóstico ====

O passo final é investigar se o modelo cumpre com as premissas do modelo linear. 

<code rsplus>

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

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! 





<WRAP center round tip 100%>

**__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.

 

</WRAP>

*/


===== Simplificando Modelos =====

<WRAP center round box 80%>
Videoaula do curso **Princípios de Planejamento e Análise de Dados**. Os conceitos abordados são os mesmos desse tutorial, desconsidere referências à disciplina. 
<WRAP center round tip 80%>

 Videoaula do curso **Planejamento e Análise de Dados**, os conceitos abordados são os mesmos, desconsidere referências à disciplina. 
{{youtube>5av4ffv89A0}}

</WRAP>
</WRAP>

Um dos procedimento de simplificar modelos é partir do modelo cheio e ir simplificando, retirando variáveis preditoras que não ajudam na explicação da variabilidade dos dados.
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.

<WRAP center round box 90%>

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

</WRAP>

<WRAP center round box 90%>
==== 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
    * verifique se a ordem da retirada de termos de mesmo nível de complexidade influencia a retirada ou manutenção dos termos finais.
    

</WRAP>

<WRAP center round box 90%>
==== Tomada de decisão  ====
<WRAP center round tip 80%>
** A diferença não é significativa: **
    * retenha o modelo mais simples
    * continue simplificando

</WRAP>
<WRAP center round tip 80%>

**A difereça é significativa: **
    * retenha o modelo complexo 
    * este é o modelo __MINÍMO ADEQUADO__

</WRAP>

</WRAP>

Já utilizamos esse procedimento no tutorial [[02_tutoriais:tutorial7:start|]], quando comparamos o modelo linear com preditora e o modelo sem nenhuma variável preditora.





===== 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]]. 

<WRAP center round box 80%>
  * baixe o arquivo  {{::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
   
</WRAP>

Para simplificar nosso tutorial vamos usar apenas as preditoras: tempo de gestação, idade da mãe e se ela é fumante ou não.

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

</code>

<code>
Call:
lm(formula = bwt ~ gestation + age + smoke + gestation:age + 
    gestation:smoke + age:smoke + gestation:age:smoke, data = bebes)

Residuals:
    Min      1Q  Median      3Q     Max 
-51.433 -10.647   0.156   9.800  50.994 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)              1.843e+02  5.443e+01   3.385 0.000735 ***
gestation               -2.262e-01  1.938e-01  -1.167 0.243542    
age                     -6.010e+00  1.942e+00  -3.095 0.002014 ** 
smokeTRUE               -1.830e+02  8.188e+01  -2.235 0.025635 *  
gestation:age            2.177e-02  6.926e-03   3.143 0.001716 ** 
gestation:smokeTRUE      6.192e-01  2.934e-01   2.110 0.035056 *  
age:smokeTRUE            3.967e+00  2.956e+00   1.342 0.179915    
gestation:age:smokeTRUE -1.397e-02  1.061e-02  -1.317 0.187994    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 16.1 on 1166 degrees of freedom
Multiple R-squared:  0.233,	Adjusted R-squared:  0.2284 
F-statistic:  50.6 on 7 and 1166 DF,  p-value: < 2.2e-16

</code>
===== Interação Tripla =====

Vamos simplificar o modelo, retirando a interação ''gestation:age:smoke'' que aparenta não ser importante.

<code rsplus>

ml01 <- lm(bwt ~ gestation + age + smoke
           + gestation:age + gestation:smoke
         + age: smoke, data = bebes) 
anova(ml01, mlfull)
summary(ml01)

</code>

===== 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.

<code rsplus>
## 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)

</code>

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.
Poderíamos continuar simplificando para garantir que não retemos nenhum termo que não é relevante para explicar o peso do bebê ao nascer. Entretanto, a menos que se tenha um bom motivo ((desenhos experimentais aninhados podem incluir a variável aninhada apenas na interação)), <wrap em>não retiramos os termos das variáveis isoladas quando ela está em algum termo de interação</wrap>.

===== Interpretação do modelo =====

O ''summary'' nos fornece as principais informações sobre o modelo mínimo adequado.  

<code rsplus>
summary(ml02)
</code>

<code rsplus>
Call:
lm(formula = bwt ~ gestation + age + smoke + gestation:age + 
    gestation:smoke, data = bebes)

Residuals:
    Min      1Q  Median      3Q     Max 
-51.978 -10.769   0.108  10.027  50.599 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)         135.598062  41.406657   3.275 0.001088 ** 
gestation            -0.055381   0.147986  -0.374 0.708301    
age                  -4.248772   1.458653  -2.913 0.003650 ** 
smokeTRUE           -75.235972  17.213833  -4.371 1.35e-05 ***
gestation:age         0.015584   0.005224   2.983 0.002911 ** 
gestation:smokeTRUE   0.239947   0.061676   3.890 0.000106 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 16.1 on 1168 degrees of freedom
Multiple R-squared:  0.2317,	Adjusted R-squared:  0.2284 
F-statistic: 70.45 on 5 and 1168 DF,  p-value: < 2.2e-16
</code>

Uma interpretação importante é com relação a variável ''smoke''. Onde foi parar o nível ''smokeFALSE''? Como é uma variável categórica de dois níveis, ''smoke''foi transformada em  variáveis indicadoras e um dos níveis deslocado para o intercepto. O que está representado no intercepto?  É a estimativa do modelo para uma mulher que não é fumante com tempo de gestação ''zero'' e  idade ''zero''. O que não faz sentido biológico nenhum. 
 
O intervalo de confiança dos coeficientes é retornado pela função ''confint'':

<code rsplus>
(coefml02 <- coef(ml02))
confint(ml02)
</code>


==== Interpretação da tabela de Anova em Modelos Multiplos ====

A função ''anova'' aplicada a um único modelo com múltiplas preditoras, nos fornece a comparação de múltiplos modelos na ordem em que as variáveis foram colocadas na fórmula. Vamos interpretar a tabela de ''anova'' do nosso modelo:

<code rsplus>
anova(ml02)
</code>

<code rsplus>
Analysis of Variance Table

Response: bwt
                  Df Sum Sq Mean Sq  F value    Pr(>F)    
gestation          1  65450   65450 252.4963 < 2.2e-16 ***
age                1    939     939   3.6241 0.0571933 .  
smoke              1  19024   19024  73.3941 < 2.2e-16 ***
gestation:age      1   1964    1964   7.5776 0.0060012 ** 
gestation:smoke    1   3923    3923  15.1354 0.0001057 ***
Residuals       1168 302757     259                       
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
</code>

A segunda linha nos diz que o modelo com ''gestação'' ao adicionar ''age'' não explica muita variação a mais. Na terceira linha a comparação é entre os modelos ''bwt ~ gestation + age'' com o modelo ''bwt ~ gestation + age + smoke'' a quarta é a comparação deste último com ''bwt ~ gestation + age + smoke + gestation:age''  e assim por diante, sempre comparando o modelo com tedos os termos anteriores e o que inclui todos os termos anteriores mais o termo que está na linha da tabela. Portanto, se colocarmos termos em outra ordem, as comparações serão outras.

<code rsplus>
ml02b <- lm(bwt ~ age + smoke + gestation + gestation:smoke 
           + gestation:age , data = bebes) 
anova(ml02b, ml02)
anova(ml02b)
</code>

<code rsplus>
Analysis of Variance Table

Response: bwt
                  Df Sum Sq Mean Sq  F value    Pr(>F)    
age                1    287     287   1.1068 0.2929867    
smoke              1  23757   23757  91.6509 < 2.2e-16 ***
gestation          1  61370   61370 236.7568 < 2.2e-16 ***
smoke:gestation    1   3580    3580  13.8130 0.0002115 ***
age:gestation      1   2307    2307   8.9001 0.0029108 ** 
Residuals       1168 302757     259                       
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

</code>

Só para entendermos o que está apresentado nessa ''anova'', vamos comparar os modelos:

  * 1. ''bwt ~ age + smoke + gestation''
  * 2. ''bwt ~ age + smoke + gestation + smoke:gestation''

<code rsplus>
ml05 <-lm(bwt ~ age + smoke + gestation, data = bebes)
ml06 <-lm(bwt ~ age + smoke + gestation + gestation:smoke, data = bebes)
anova(ml05, ml06)
</code>

Pode haver pequenas variações nos valores por conta arredondamentos. O importante aqui é que um termo pode ser significativo ou não dependendo da ordem que for colocado, principalmente se há alguma colinearidade entre as variáveis incluídas. Ou seja, o termo que é colocado antes explica a variação que o termo que vem depois poderia explicar também!

===== Diagnóstico do modelo =====

O diagnóstico das premissas do modelo é importante, para mais informações veja o tutorial da disciplina [[http://labtrop.ib.usp.br/doku.php?id=cursos:planeco:roteiro:07a-clasrcmdr| Princípios de Planejamento e Análise de Dados]] sobre o assunto. O basico pode ser interpretado nos gráficos que são feitos por padrão se usamos a função ''plot'' no objeto de modelo:

<code rsplus>
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)
</code>


Estando tudo certo com nosso modelo podemos passar para outras fases como preparar gráficos e interpretar os resultados. 


====== Videoaula Síncrona ======

<WRAP center round box 80%>
Aula síncrona gravada pelo Google Meet em 05 de outubro de 2020
<WRAP center round tip 80%>
{{youtube>F2TVNNZscmM}}
</WRAP>

</WRAP>



/*

===== 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:
<code rsplus>
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))
</code>

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:
<code rsplus>
plot(Pollution~Industry, data=poluicao)
abline(pol.m1, col="blue")
</code>

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:
<code rsplus>
 abline(lm(Pollution~Industry, data=poluicao, subset=Industry<max(Industry)), col="red")
</code>

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:
<code rsplus>
pol.m0 <- lm(Pollution~1, data=poluicao)
##compare:
anova(pol.m0,pol.m1)
anova(pol.m1)
</code>
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:
<code rsplus>
mean(poluicao$Pollution)
sd(poluicao$Pollution)
</code>

Um gráfico com as linhas dos dois modelos concorrentes:
<code rsplus>
plot(Pollution~Industry, data=poluicao)
abline(pol.m1, col="blue")
abline(h=mean(poluicao$Pollution),col="red")
</code>

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:
<code rsplus>
pol.m3 <- update(pol.m2,.~.+Wind)
anova(pol.m2,pol.m3)
anova(pol.m3)
</code>

Não houve ganho importante de explicação. Ficamos com o modelo anterior. Vamos inspecionar os gráficos de diagnóstico:

<code rsplus>
par(mfrow=c(2,2))
plot(pol.m2)
par(mfrow=c(1,1))
</code>

Há um ponto discrepante (41). vamos identificá-lo nos dados e nos gráficos:
<code rsplus>
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))
</code>

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:
<code rsplus>
anova(lm(Pollution~Temp+Wind+Industry, data=poluicao))
anova(lm(Pollution~Industry+Wind+Temp, data=poluicao))
</code>

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//: 
<code rsplus>
library(lattice)
temp.shingle <- equal.count(poluicao$Temp,number=4)
summary(temp.shingle)
</code>

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'')
<code rsplus>
xyplot(Pollution~Industry|temp.shingle,data=poluicao,
   panel=function(x,y,...){
     panel.xyplot(x,y,...)
     panel.abline(lm(y~x),...)
   }
       )
</code>

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:
<code rsplus>
xyplot(Pollution~Industry|equal.count(poluicao$Wind,number=4),data=poluicao,
   panel=function(x,y,...){
     panel.xyplot(x,y,...)
     panel.abline(lm(y~x),...)
   }
       )
</code>

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:
<code rsplus>
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)
</code>

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:
<code rsplus>
##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)
</code>

Em seguida fazemos três graficos com os valores selecionados po faixa de temperatura
<code rsplus>
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))
</code> 

*/

/*
===== O Modelo por trás da ANOVA =====

Neste tópico 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:
<code rsplus>
mudas <- read.csv("altura-mudas.csv")
mudas$bloco <- as.factor(mudas$bloco)
mudas$substrato <- as.factor(mudas$substrato)
</code>

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:
<code rsplus>
summary(mudas.m1)
## Coeficientes do modelo
cf.mud1 <- coef(mudas.m1)
cf.mud1
</code>

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'':
<code rsplus>
head(model.matrix(mudas.m1)%*%cf.mud1)
tail(model.matrix(mudas.m1)%*%cf.mud1)
</code>

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)

<box red| CUIDADO!> 
Prossiga apenas após estar certo(a) do significado dos coeficientes neste primeiro modelo 
</box>

Agora vamos criar um novo modelo, que terá também o efeito do substrato:

     mudas.m2 <- update(mudas.m1,.~.+substrato)

Inspecione os coeficientes:
<code rsplus>
cf.mud2 <- coef(mudas.m2)
cf.mud2
</code>

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!
*/