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.

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

**__Formulas Estatísticas__**

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 utilizar
a função '''I()''':
  * ''' 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))**//;

**__Delineamentos Experimentais__**
^ Expressão ^ Significado ^
| //Y ~ X// | Modele //Y// como função estatística de //X// |
| //A + B// | inclui ambos os fatores //A// e //B// |
| //A - B// | inclui todos os efeitos em //A//, exceto os que estão em //B// |
| //A * B// | //A + B + A:B// |
| //A / B// | //A + B %in% (A)// modelos hierárquicos |
| //A:B// | efeito da interação entre os fatores //A// e //B// |
| //B %in% A// | efeitos de //B// dentro dos níveis de //A// |
| //A^m// | todos os termos de //A// cruzados até à ordem //m// | 




</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 ((no exercício terão que usar os dados brutos e todas as variáveis)).

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

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

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

<code rsplus>
summary(ml02)
confint(ml02)
anova(ml02)

</code>

===== Diagnóstico do 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>

<WRAP center round important 60%>
Os dados desse estudo serão usados também no exercício, porém lá, vamos partir dos dados brutos com mais variáveis
</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!
*/
ecoR

Ferramentas do usuário

Ferramentas do site

Ferramentas da página
