* [[:02_tutoriais:tutorial7b:start|Tutorial]] * [[:01_curso_atual:exercicios8| Exercícios]] ====== 7b. Modelos Lineares Múltiplos ====== 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. {{youtube>89hcuYp9oWo}} 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. **__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 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)''; ===== 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 ===== 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. {{youtube>Mx9skekN6e8}} 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. **__Delineamentos Experimentais__** ^ 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| /* ==== 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. */ ===== Simplificando Modelos ===== 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. Videoaula do curso **Planejamento e Análise de Dados**, os conceitos abordados são os mesmos, desconsidere referências à disciplina. {{youtube>5av4ffv89A0}} 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. ==== 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 * verifique se a ordem da retirada de termos de mesmo nível de complexidade influencia a retirada ou manutenção dos termos finais. ==== 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__ 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]]. * 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 Para simplificar nosso tutorial vamos usar apenas as preditoras: tempo de gestação, idade da mãe e se ela é fumante ou não. 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) 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 ===== 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. 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)), não retiramos os termos das variáveis isoladas quando ela está em algum termo de interação. ===== Interpretação do modelo ===== O ''summary'' nos fornece as principais informações sobre o modelo mínimo adequado. summary(ml02) 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 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'': (coefml02 <- coef(ml02)) confint(ml02) ==== 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: anova(ml02) 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 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. ml02b <- lm(bwt ~ age + smoke + gestation + gestation:smoke + gestation:age , data = bebes) anova(ml02b, ml02) anova(ml02b) 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 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'' ml05 <-lm(bwt ~ age + smoke + gestation, data = bebes) ml06 <-lm(bwt ~ age + smoke + gestation + gestation:smoke, data = bebes) anova(ml05, ml06) 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: 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) Estando tudo certo com nosso modelo podemos passar para outras fases como preparar gráficos e interpretar os resultados. ====== Videoaula Síncrona ====== Aula síncrona gravada pelo Google Meet em 05 de outubro de 2020 {{youtube>F2TVNNZscmM}} /* ===== 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 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: 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! */