Ferramentas do usuário

Ferramentas do site


02_tutoriais:tutorial7:start

Diferenças

Aqui você vê as diferenças entre duas revisões dessa página.

Link para esta página de comparações

02_tutoriais:tutorial7:start [2022/06/21 11:55]
adalardo [R² Ajustado]
02_tutoriais:tutorial7:start [2023/09/11 15:58]
Linha 1: Linha 1:
-<WRAP tabs> 
-  * [[02_tutoriais:​tutorial7:​start|Tutorial]] 
-  * [[01_curso_atual:​exercicios7| Exercícios]] 
-  * [[03_apostila:​06-modelos| Apostila]] ​ 
-</​WRAP>​ 
-====== 7a. Regressão Linear Simples ​ ====== 
- 
-A primeira parte desse tutorial é baseado no [[http://​labtrop.ib.usp.br/​doku.php?​id=cursos:​planeco:​roteiro:​08-lm_rcmdr| tutoria de modelos lineares da disciplina Princípios de Planejamento e Análise de Dados]], inclusive as vídeoaulas. Aqui iremos focar no código que estava subjacente ao tutorial. 
-<WRAP center round tip 60%> 
-A videoaula gravada no google meet no dia 02 de outubro de 2020 está ao final do tutorial. Dê preferência para as videoaulas do curso de **Principios de Planejamento e Análise de Dados** que estão colocadas ao longo do tutorial. Eles tratam o tema de modelos lineares de forma mais sucinta e tiveram alguma edição. Desconsiderem nesses vídeos as referências à disciplina. 
-</​WRAP>​ 
- 
-===== Modelos Lineares ===== 
- 
-{{:​02_tutoriais:​tutorial7:​lmComics.jpg?​400 ​ |}} 
- 
-Os modelos lineares são uma generalização dos testes de hipótese clássicos mais simples. Uma regressão linear, por exemplo, só pode ser aplicada para dados em que tanto a variável preditora quanto a resposta são contínuas, enquanto uma análise de variância é utilizada quando a variável preditora é categórica. Os modelos lineares não têm essa limitação,​ podemos usar variáveis contínuas ou categóricas indistintamente. 
-<WRAP center round box 80%> 
-__**Videoaula Modelo Linear I**__ 
-O vídeo é proveniente de outra disciplina, desconsidere qualquer referência a ela. 
-<WRAP center round tip 80%> 
-{{youtube>​b4VgLr6loGE}} 
- 
-</​WRAP>​ 
- 
- 
-</​WRAP>​ 
- 
-No quadro abaixo estão listados alguns dos testes clássicos frequentistas. Estes testes foram criados para diferentes naturezas de variáveis respostas e preditoras. Você já refletiu sobre a natureza das variáveis do seu estudo? Esse é um passo importante na tomada de decisão da análise adequada, assim como seu acoplamento com a hipótese de trabalho é fundamental! 
- 
- 
-<WRAP center round box 80%> 
-{{  :​02_tutoriais:​tutorial7:​tabTestes.png?​700 ​ |}} 
-</​WRAP>​ 
- 
-O modelo linear é uma unificação que incorpora os testes da tabela acima que tem a __** variável resposta contínua**__. Portanto, já não há mais necessidade dos nomes: //​teste-t//,​ //Anova//, //Anova Fatorial//, //​Regressão Simples//, //​Regressão Múltipla//,​ //Ancova// entre muitos outros nomes de testes que foram incorporados nos modelos lineares. ​   
-Isso não livra o bom usuário de estatística de entender a natureza das variáveis que está utilizando. Isso continua sendo imprescindível para tomar boas decisões ao longo do processo de análise dados e interpretação dos resultados. 
- 
- 
-===== Simulando dados ===== 
-  ​ 
-Vamos começar com um exemplo simples de regressão, mas de forma diferente da usual. Vamos inverter o procedimento para entender bem o que os modelos estatísticos nos dizem e como interpretar os resultados produzidos. Para isso, vamos inicialmente gerar dados fictícios a partir de um universo conhecido. Esses dados terão dois componentes:​ uma estrutura determinística e outra aleatória. A primeira está relacionada ao processo de interesse do estudo e relaciona a variável resposta à preditora. No caso, essa estrutura é uma relação linear e tem a seguinte forma: 
- 
- 
-$$ y = {\alpha} + {\beta} x$$ 
- 
-O componente aleatório é expresso por uma variável probabilística Gaussiana da seguinte forma: ​ 
- 
-$$ \epsilon = N(0, \sigma) $$ 
-  
-Portanto, nossos dados serão uma amostra de uma população com a seguinte estrutura: 
- 
-$$ y = {\alpha} + {\beta} x + \epsilon$$ 
- 
- 
-Parece complicado, mas é simples gerar dados aleatórios com essa estrutura do R. Vamos definir primeiro quais são os parâmetros que estão na nossa população,​ ou seja qual o valor de $\alpha$ e $\beta$ da relação entre ''​y''​ e ''​x''​ na população. Além disso, vamos definir também qual a variabilidade associada a essa relação, o nosso $\epsilon$. 
- 
-$$ y = 5.3 + 0.12 x + N(0, 5) $$  
- 
- 
-Antes de gerar os dados aleatórios,​ vamos utilizar uma ferramenta que define a raiz da semente aleatória que o R irá usar. Com isso, apesar dos dados gerados serem provenientes de uma amostra aleatória, todos que utilizarem a mesma semente terão os mesmo valores amostrados. Em seguida vamos criar uma sequência para representar a variável preditora ''​x''​ e, a partir da relação acima, calcular o ''​y0'',​ que são os valores associados a essa relação determinística com ''​x''​ e também criar um vetor ''​res''​ que define a variabilidade do nossos dados: 
- 
-<code rsplus> 
-set.seed(1) 
-x1 <- round(seq(12,​ 220, len = 15), 1) 
-y0 <- 10.3 + 0.12 * x1 
-res <- rnorm(length(x1),​ 0, 5) 
-y1 <- y0 + res 
-xm <- mean(x1) 
-ym <- mean(y1) 
-</​code>  ​ 
- 
- 
-Vamos olhar esses valores em um gráfico utilizando o procedimento que aprendemos no tutorial [[02_tutoriais:​tutorial5b:​start|]]:​ 
- 
-<code rsplus> 
-par(mar = c(4, 4, 2, 2), cex.lab = 1.5, cex.axis = 1.5, las = 1, bty = "​n"​) 
-plot(x1, y1, type = "​n",​ axes = FALSE, ann = FALSE, ylim = range(y1), xlim = range(x1)) 
-rect(par()$usr[1],​ par()$usr[3],​ par()$usr[2],​ par()$usr[4], ​ col = rgb(0, 0, 0, 0.15)) 
-axis(1) 
-mtext(text = "​Variável preditora (x1)", side = 1, line = 3, cex = 1.5) 
-axis(2) 
-mtext(text = "​Variável resposta (y1)", side = 2, line = 3, cex = 1.5, las =0) 
-cores <- c(rgb(1, 0, 0, 0.3), rgb(0, 0, 1, 0.3)) 
-points(x1, y0, pch = 16, cex = 0.8, col = cores[1] ) 
-points(x1, y1, pch = 19, col = cores[2]) 
-legend("​bottomright",​ legend = c("y0 = 10.3 + 0.12 x1", "y1 = y0 + N(0, 5)"), bty = "​n",​ col = cores, pch = 19) 
-</​code>​ 
- 
-<WRAP center round box 80%> 
-{{  :​02_tutoriais:​tutorial7:​plotxy01.png?​600 ​ |}} 
-</​WRAP>​ 
- 
-==== Estimando os parâmetros ==== 
- 
-{{:​02_tutoriais:​tutorial7:​convAxisLeg.png?​600 ​ |}} 
- 
-Os valores em ''​x1''​ e ''​y1'',​ a partir desse ponto, são os dados, provenientes de uma amostra aleatória de ''​y1''​ e de sua relação determinística com ''​x1''​. Vamos esquecer que geramos os dados. Vamos partir do ponto que temos os dados coletados e precisamos fazer as inferências relacionadas à relação ''​y1 ~ x1''​. ​ 
- 
-Com esses dados em mãos precisamos estimar os parâmetros que definem, se é que existe, a relação ''​y0 ~ x1''​. Como não temos ''​y0'',​ apenas a amostra com ''​y1'',​ eles são nossa melhor pista para verificar se essa relação existe. A partir deles então, buscamos estimar os parâmetros da população $\alpha$, $\beta$ e $\epsilon$. Há várias maneiras de estimarmos os parâmetros de uma relação linear. Podemos usar o clássico <wrap em>​método de mínimos quadrados</​wrap>​ que, nesse caso específico,​ também é a solução de <wrap em>​máxima verossimilhança</​wrap>​. Ambos são medidas que indicam o quanto o modelo, no caso uma relação linear, está ajustado aos dados. Buscamos então os parâmetros que melhor descrevem a relação.  ​ 
- 
-Vamos ilustrar aqui como esses métodos funcionam de forma intuitiva, usando a solução numérica, também chamada de //força bruta//. Apesar de existir solução algébrica, esse método é muito útil em casos mais complexos e bastante intuitivo. Tudo começa com o ponto de fulcro, o centro de massa dos nossos dados, as medias de ''​x1''​ e ''​y1''​. Esse ponto, necessariamente está na reta que melhor representa a relação ''​y1 ~ x1''​. 
-Em seguida, podemos testar diferentes inclinações,​ passando pelo fulcro e buscar a inclinação que melhor representa a relação. A medida de melhor ajuste pode variar, podemos usar a soma dos desvios quadráticos ou a verossimilhança. Nesse caso, ambos irão chegar ao mesmo resultado. ​ 
- 
-Um modelo possível é a inclinação ''​b = 0'',​ abaixo o gráfico dessa relação e o cálculo da soma dos desvios quadráticos para esse modelo: 
- 
-<code rsplus> 
-plot(x1, y1, type = "​n",​ axes = FALSE, ann = FALSE, , ylim = range(y1), xlim = range(x1)) 
-rect(par()$usr[1],​ par()$usr[3],​ par()$usr[2],​ par()$usr[4], ​ col = rgb(0, 0, 0, 0.15)) 
-axis(1) 
-mtext(text = "​Variável preditora (x1)", side = 1, line = 3, cex = 1.5) 
-axis(2) 
-mtext(text = "​Variável resposta (y1)", side = 2, line = 3, cex = 1.5, las =0) 
-points(x1, y1, pch = 19, col = cores[2], cex = 1.5) 
-points(xm, ym, pch=19) 
-points(xm, ym, cex = 2.5) 
-text(xm, ym-2 , labels= "​fulcro",​ cex = 1.5) 
-abline(mean(y1),​ 0, col= "​white",​ lwd=2 ) 
-segments(x0 = x1, y0 = y1, x1= x1, y1= mean(y1) , lty=2) 
-text(70, 40, paste("​Desvio Quadratico =", round( sum((y1 - mean(y1))^2),​2)),​ cex= 1.2) 
-points(xm, ym, pch=19) 
-points(xm, ym, cex = 2.5) 
-</​code>​ 
- 
- 
-<WRAP center round box 80%> 
-{{  :​02_tutoriais:​tutorial7:​modNull.png?​600 ​ |}} 
-</​WRAP>​ 
- 
- 
-Agora, precisamos encontrar entre as muitas possibilidades de inclinação,​ aquela que <wrap em>​minimiza</​wrap>​ o valor de soma dos desvios quadráticos. ​ 
-Por sorte já sabemos fazer iterações e gráficos no R! Vamos utilizar o gráfico anterior como base e fazer muitos gráficos com diferentes inclinações e ao mesmo tempo calcular a soma dos desvios ao quadrado. Antes vamos calcular os valores para diferentes inclinações:​ 
- 
-<code rsplus> 
-bsim <- seq(-0.05, 0.29, 0.005) 
-asim <- ym - bsim * xm 
-nsim <- length(bsim) 
-ML <- rep(NA, nsim) 
-MQ <- rep(NA, nsim) 
-for (i in 1:nsim) 
-{ 
-    spred <- asim[i] +  bsim[i] * x1 
-    MQ[i] <- sum((spred - y1)^2) 
-    ML[i] <-  prod(dnorm(y1,​ mean = spred, sd = 5, log = FALSE)) 
-} 
-estima <- data.frame(alfa = asim, beta = bsim, ML = ML, MQ = MQ, logML = log10(ML)) 
-</​code>​ 
- 
-Acima, criamos uma sequencia de valores para representar diferentes inclinações ''​bsim'',​ estimamos o intercepto, ''​asim'',​ a partir dessas inclinações. Em seguinda utilizamos esses coeficientes para calcular o valor de soma quadrática,​ ''​MQ'',​ e verossimilhança,​ ''​ML''​. ​ 
- 
- 
-Abaixo um gráfico animado. A única função que ainda não conhecemos é ''​Sys.sleep''​((Consulte a documentação!)),​ que é uma forma de fazer o R dar um cochilada durante os ciclos, sem isso não seria possível ver os quadros.  ​ 
- 
-<code rsplus> 
-graphics.off() 
-x11(width = 12, heigh = 6) 
-par(mfrow = c(1,2),mar= c(5, 5, 2, 2), las=1, cex.axis= 1.2, cex.lab=1.2,​ bty= "​n"​) 
-for (i in 1: length(bsim)) 
-{ 
-    spred = asim[i] + bsim[i] * x1 
-    plot(x1, y1, type = "​n",​ axes = FALSE, ann = FALSE, , ylim = range(y1), xlim = range(x1)) 
-    rect(par()$usr[1],​ par()$usr[3],​ par()$usr[2],​ par()$usr[4], ​ col = rgb(0, 0, 0, 0.15)) 
-    axis(1) 
-    mtext(text = "​Variável preditora (x1)", side = 1, line = 3, cex = 1.5) 
-    axis(2) 
-    mtext(text = "​Variável resposta (y1)", side = 2, line = 3, cex = 1.5, las =0) 
-    points(x1, y1, pch = 19, col = cores[2], cex = 1.75) 
-    points(xm, ym, pch=19) 
-    points(xm, ym, cex = 2.5) 
-    abline(asim[i],​ bsim[i], col= cores[1], lwd = 4, lty = 1) 
-    segments(x0 = x1, y0 = y1, x1= x1, y1= spred, lty=2) 
-    points(xm, ym, pch=19) 
-    points(xm, ym, cex = 2.5) 
-    plot(NULL, NULL, type="​n",​ axes = FALSE, ann = FALSE, xlim =range(bsim),​ ylim = range(estima$MQ)) ​ 
-    rect(par()$usr[1],​par()$usr[3],​ par()$usr[2],​ par()$usr[4],​ col="​gray80"​) 
-    points(bsim[1:​i],​ estima$MQ[1:​i],​ pch = 19, cex = 1.75, col = rgb(1, 0.4, 0.8, 0.9)) 
-    axis(1) 
-    axis(2) 
-    mtext("​Inclinação",​ 1, at = mean(bsim), line=2.5, cex = 1.5) 
-    mtext("​Soma dos desvios quadráticos",​ 2, at = par()$usr[3]+ ((par()$usr[4] - par()$usr[3])/​2),​ line = 4, cex = 1.5, las =0) 
-    Sys.sleep(0.2) 
-} 
-</​code>​ 
- 
- 
-<WRAP center round box 100%> 
-{{  :​02_tutoriais:​tutorial7:​minquadSmall_plot.gif?​800 ​ |}} 
-</​WRAP>​ 
- 
-Com o resultado dessa simulação podemos buscar o valor de inclinação que melhor representa a relação: 
- 
-<code rsplus> 
-(bEst <- bsim[which.min(estima$MQ)]) 
- 
-</​code>​ 
- 
-Como temos um ponto que conhecemos do nosso modelo, o fulcro, podemos calcular o intercepto a partir dele: 
- 
-$$ym = aEst + bEst * xm$$ 
- 
-$$aEst = ym - (bEst * xm)$$ 
- 
-<code rsplus> 
-(aEst <- ym - (bEst * xm)) 
-</​code>​ 
- 
- 
-Nossas estimativas parecem muito boas!  
-Quando construímos modelos, não precisamos fazer a simulação para fazer estimativas,​ a função ''​lm''​ faz isso para nós. 
- 
-===== Modelos Lineares ===== 
- 
-Para definir o modelo linear na função ''​lm''​ precisamos utilizar e expressão em formula do R: 
- 
-''​lm(y ~ x, data = xy)''​ 
- 
-Onde ''​y''​ representa a variável resposta e ''​x''​ a preditora. Além disso, podemos utilizar o argumento ''​data''​ para indicar o objeto onde estão armazenadas as variáveis ''​x''​ e ''​y''​. Vamos criar nosso primeiro modelo e verificar o resultado: 
- 
-<code rsplus> 
-lmxy01 <- lm(y1 ~ x1) 
-class(lmxy01) 
-str(lmxy01) 
-</​code> ​ 
- 
-  ​ 
-O objeto de modelo é bastante complexo. Para acessar as informações utilizamos as funções extratoras que retiram do objeto as informações solicitadas,​ por exemplo, a classe em  ''​class(lmxy1)''​. Vamos avaliar as principais informações que podem ser extraídas de um objeto de modelo da classe ''​lm''​. 
- 
- 
-==== Coeficientes do Modelo ==== 
- 
-Os coeficientes estimados pelo modelo podem ser obtidos pela função ''​coef''​. A função ''​confint''​ retorna os intervalos de confiança associados à essas estimativas:​ 
- 
-<code rsplus> 
-coef(lmxy01) 
-confint(lmxy01) 
-</​code>​ 
- 
- 
-==== Predito pelo Modelo ==== 
- 
-O predito pelo modelo, ou sua esperança, é o valor que a reta do modelo prevê para cada observação em ''​x1''​. Podemos usar a estimativas dos coeficientes do modelo e a equação da reta para estimar esses valores, ou usar a função ''​predict'':​ 
- 
-<code rsplus> 
-(predxy01 <- coef(lmxy01)[1] + coef(lmxy01)[2] * x1) 
-predict(lmxy01) 
-</​code>​ 
- 
- 
- 
-==== Resíduos do Modelo ==== 
- 
-Os resíduos do modelo são o quanto cada observação se afasta do predito pelo modelo. Ou seja, o observado em ''​y1''​ menos o predito pelo modelo para cada observação em ''​x1'':​ 
- 
-<​code>​ 
-(resxy01 <- y1 - predxy01) 
-residuals(lmxy01) 
-</​code> ​ 
- 
-==== Resumo do Modelo ==== 
- 
-O ''​summary''​ organiza as principais informações do modelo. Uma de nossas metas é entender todos os elementos que são organizados pelo ''​summary''​.  ​ 
- 
-<code rsplus> 
-summary(lmxy01) 
-</​code>​ 
- 
-Abaixo a saída do ''​summary'':​ 
- 
-<​code>​ 
-Call: 
-lm(formula = y1 ~ x1) 
- 
-Residuals: 
-    Min      1Q  Median ​     3Q     ​Max ​ 
--11.449 ​ -3.645 ​  ​1.079 ​  ​2.792 ​  ​7.386 ​ 
- 
-Coefficients:​ 
-            Estimate Std. Error t value Pr(>​|t|) ​   ​ 
-(Intercept) 10.97202 ​   2.81569 ​  ​3.897 ​ 0.00184 **  
-x1           ​0.11855 ​   0.02124 ​  5.582 8.89e-05 *** 
---- 
-Signif. codes: ​ 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
- 
-Residual standard error: 5.28 on 13 degrees of freedom 
-Multiple R-squared: ​ 0.7056,​ Adjusted R-squared: ​ 0.683  
-F-statistic:​ 31.16 on 1 and 13 DF,  p-value: 8.892e-05 
-</​code>​ 
- 
- 
-Após a fórmula do modelo são apresentados os quartis associados aos resíduos e as estimativas dos coeficientes na coluna ''​Estimate'':​ 
- 
-<code rsplus> 
-quantile(resxy01) 
-coef(lmxy01) 
-</​code>​ 
- 
-Para entender o que são os outros valores que aparecem na tabela ''​Coefficients'',​ precisamos entender como a estatística clássica frequentista entende a imprecisão das estimativas. ​ 
- 
- 
-===== Múltiplos Experimentos ===== 
- 
- 
-Uma das bases da teoria da estatística frequentista é que uma amostra e seus resultados são apenas uma realização dentre os possíveis resultados provenientes de uma população real, a qual não temos acesso. Esta abstração da realização de múltiplos experimentos nos ajuda a entender as estimativas dos modelos. Vamos verificar o que acontece com nossas estimativas ao realizar muitas vezes amostras similares provenientes da mesma população,​ vamos guardar as estimativas e as somas quadráticas dos desvios a cada iteração. Para isso, vamos criar os vetores para armazenar as informações antes abrir o ''​for( )'':​ 
- 
-<code rsplus> 
-nSimula <- 1000 
-aEst <-  rep(NA, nSimula) 
-bEst <-  rep(NA, nSimula) 
-desQuad <-  rep(NA, nSimula) 
-</​code>​ 
- 
-Vamos resgatar((esse vetores foram construídos logo no início do tutorial!)) a relação determinística do nossa população:​ 
- 
-<​code>​ 
-x1 <- round(seq(12,​ 220, len = 15), 1) 
-y0 <- 10.3 + 0.12 * x1 
-</​code>​ 
- 
-Agora criamos a iteraçao com ''​nSimula''​ ciclos. A cada ciclo a sequência de operações é: 
- 
-  *1. Cria valores aleatórios de uma normal com média zero e desvio padrão igual a 5; 
-  *2. Soma esses valores ao valor ''​y0 '',​ proveniente da relação ''​y ~ x''​ da população e cria os dados ''​ySim''; ​ 
-  *3. Cria o modelo com o ''​ySim ~ x1'';​ 
-  *4. Guarda os coeficientes do modelo em ''​cSim'';​ 
-  *5. Armazena as estimativas dos coeficientes nos respectivos vetores, ''​aEst''​ e ''​bEst'',​ na posição ''​i'';​ 
-  *6. Calcula e armazena a soma dos desvios quadráticos do modelo em ''​desQuad''​ na posição ''​i''​. 
- 
-<code rsplus> 
-for(i in 1:nSimula) 
-{ 
-rSim <- rnorm(length(x1),​ 0, 5) 
-ySim <- y0 + rSim 
-lmSim <- lm(ySim ~ x1) 
-cSim <- coef(lmSim) 
-aEst[i] <- cSim[1] 
-bEst[i] <- cSim[2] 
-desQuad[i] <- sum(residuals(lmSim)^2) 
-} 
-</​code>​ 
- 
- 
-Vamos olhar agora as médias e a distribuição dos valores estimados em um histograma: 
- 
-<code rsplus> 
-(med_aEst <- mean(aEst)) 
-(med_bEst <- mean(bEst)) 
-## 
-hist(aEst, main = "​Estimativa do intercepto",​ ylab = ""​) 
-abline(v = med_aEst, col = "​red",​ lwd = 4, lty = 2) 
-mtext("​Frequência",​ side = 2, line = 3.5, cex = 1.5, las = 0) 
-## 
-hist(bEst, main = "​Estimativa da inclinação",​ ylab = ""​) 
-abline(v = med_bEst, col = "​red",​ lwd = 4, lty = 2) 
-mtext("​Frequência",​ side = 2, line = 3.5, cex = 1.5, las = 0) 
-</​code>​ 
- 
- 
-<WRAP center round box 100%> 
-{{  :​02_tutoriais:​tutorial7:​LmEstHist.png?​800 ​ |}} 
-</​WRAP>​ 
- 
- 
-Apesar de nossa primeira estimativa dos parâmetros ter sido muito boa, poderíamos ter feito uma amostra com valores mais diferentes. Não é depressível a chance de nossa amostra gerar valores de intercepto menores que ''​8''​ ou de inclinação maiores que ''​0.14''​. Essa imprecisão associada a estimativa dos parâmetros é o <wrap em>erro padrão</​wrap>​. O erro padrão é o desvio padrão dessa distribuição de valores das estimativas se pudéssemos refazer a amostra muitas vezes na população,​ como acabamos de fazer com a nossa simulação! 
- 
-Vamos calcular o desvio padrão dos múltiplos experimentos:​ 
- 
-<code rsplus> 
-(med_aEst <- mean(aEst)) 
-(med_bEst <- mean(bEst)) 
-</​code>​ 
-  
-Quando só temos um experimento,​ que é a situação real de um estudo, é possível calcular o erro padrão, por exemplo, da inclinação,​ com a fórmula: 
- 
-<WRAP center round box 60%> 
- 
-$$\beta_{se} = \sqrt{\frac{\frac{1}{n-2}\sum_{i=1}^{n}\hat{\epsilon}_{i}^{2}}{\sum_{n=1}^{n}(x_i - \bar{x})}} $$ 
- 
-</​WRAP>​ 
- 
-O valor calculado a partir dessa formula é apresentado no tabela do ''​summary''​ em ''​Std. Error''​. ​ Representa a estimativa do desvio padrão da distribuição dos coeficientes estimados em múltiplos experimentos,​ análogo ao que fizemos a partir das simulações. ​ Em resumo, o erro padrão representa a imprecisão na estimativa dos coeficientes,​ baseado no conceito que se repetíssemos nosso experimento muitas vezes o desvio padrão da distribuição dos valores obtidos para o coeficiente é o erro padrão.  ​ 
- 
-Agora que já sabemos o que é o <wrap em>erro padrão</​wrap>​ das estimativas dos coeficientes,​ vamos avaliar os outros valores da tabela ''​Coefficients''​ do summary: 
- 
-  * ''​t value'': ​ é a razão entre a estimativa apresentada na coluna ''​Estimate''​ dividida pelo ''​Std. Error''​. 
-  * ''​Pr(>​|t|)'':​ probabilidade do teste t bicaudal para esse valor. ​   
- 
-Este último valor está relacionada à probabilidade de incorremos em erro ao afirmar que a estimativa do coeficiente em questão é diferente de zero. 
- 
-===== Erro Padrão do Modelo ===== 
- 
-O próximo valor temos no ''​summary''​ é o ''​Resídual Standard Error''​. Ele nada mais é do que a estimativa dos desvios padrão dos resíduos das amostras. Ou seja, nossa estimativa do desvio padrão da população que modelamos, no caso,  ''​sd = 5''​. ​ 
- 
- 
-Nesse ponto já deve estar capacitado para fazer a leitura de grande parte do ''​summary''​ de um modelo linear. Confira se entendeu os elementos que aparecem até ''​Residual standard error'',​ inclusive. 
- 
-<code rsplus> 
-summary(lmxy01) 
- 
-</​code>​ 
- 
- 
- 
-<​code>​ 
-Call: 
-lm(formula = y1 ~ x1) 
- 
-Residuals: 
-    Min      1Q  Median ​     3Q     ​Max ​ 
--11.449 ​ -3.645 ​  ​1.079 ​  ​2.792 ​  ​7.386 ​ 
- 
-Coefficients:​ 
-            Estimate Std. Error t value Pr(>​|t|) ​   ​ 
-(Intercept) 10.97202 ​   2.81569 ​  ​3.897 ​ 0.00184 **  
-x1           ​0.11855 ​   0.02124 ​  5.582 8.89e-05 *** 
---- 
-Signif. codes: ​ 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
- 
-Residual standard error: 5.28 on 13 degrees of freedom 
-Multiple R-squared: ​ 0.7056,​ Adjusted R-squared: ​ 0.683  
-F-statistic:​ 31.16 on 1 and 13 DF,  p-value: 8.892e-05 
- 
-</​code>​ 
- 
- 
- 
- 
- 
- 
- 
- 
- 
- 
- 
-===== Tabela de Anova de uma Regressão ===== 
- 
-<WRAP center round box 100%> 
-<WRAP center round tip 80%> 
-Video na disciplina de Princípios de Planejamento e Análise de Dados. Desconsidere qualquer referência à disciplina. O tema tratado é a partição de variação dos dados. ​ 
-{{youtube>​C4urUFRGDvo}} 
- 
-</​WRAP>​ 
-</​WRAP>​ 
- 
-Os modelos, independente da natureza da variável preditora, podem ser analisados através do método de partição de variância que aprendemos no roteiro de [[02_tutoriais:​tutorial6b:​start|]]. Vamos utilizar os dados simulados no tópico [[#​simulando_dados|Simulando dados]] para exemplificar essa propriedade. Garanta que têm na sua sessão do R os objetos ''​x1''​ e ''​y1'',​ caso contrário retorne ao tópico e recrie os objetos a partir do ''​set.seed(1)''​. 
- 
-Veja o que acontece quando a função ''​anova''​ é aplicada ao objeto de modelo ''​lm'':​ 
- 
-<code rsplus> 
-anova(lmxy01) ​ 
-</​code>​ 
- 
-O resultado é uma tabela de anova com a partição da variância dos dados do modelo: 
- 
-<​code>​ 
-Analysis of Variance Table 
- 
-Response: y1 
-          Df Sum Sq Mean Sq F value    Pr(>​F) ​   ​ 
-x1         1 868.51 ​ 868.51 ​ 31.158 8.892e-05 *** 
-Residuals 13 362.37 ​  ​27.87 ​                     ​ 
---- 
-Signif. codes: ​ 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
-</​code>​ 
- 
-Nessa tabela temos todos os componentes da partição da variabilidade dos dados, da maneira clássica apresentada no roteiro de [[02_tutoriais:​tutorial6b:​start|Anova]]. Note que estamos usando a função ''​anova''​ no R, não para fazer uma //análise de variância//​ clássica, mas para construir uma tabela da partição da variação dos dados. A primeira diferença que notamos nessa tabela é que a variabilidade total não é apresentada. Nem precisaria, tendo em vista a propriedade aditiva da partição das variabilidades. A primeira linha é relacionada a variabilidade explicada pela variável preditora. A segunda linha é a variabilidade associada ao resíduo, ou seja, a variabilidade não explicada pela preditora. A razão das médias quadráticas,​ ou variância, ​ é a estatística <wrap em>F de Fisher</​wrap>​. Fazemos a leitura da mesma maneira que fazemos para a análise de variância, inclusive para calcular o <wrap em>​coeficiente de determinação</​wrap>​ que é a variação explicada sobre a variação total, utilizando as somas quadráticas. 
-Vamos representar os gráficos associados a representação dessas duas variações nos dados: 
- 
-<code rsplus> 
-cores <- c(rgb(1, 0, 0, 0.3), rgb(0, 0, 1, 0.3)) 
-x11(width = 12, height = 14) 
-par(mar = c(4, 5, 2, 2), las = 1, cex = 1.5, cex.lab = 1.3, cex.axis = 1.3) 
-layout(matrix(c(1,​ 3, 2, 3), ncol = 2), heights = c(4, 6)) 
-layout.show(3) 
-plot(x1, y1, type = "​n",​ axes = FALSE, ann = FALSE, , ylim = range(y1), xlim = range(x1)) 
-rect(par()$usr[1],​ par()$usr[3],​ par()$usr[2],​ par()$usr[4], ​ col = rgb(0, 0, 0, 0.15)) 
-axis(1) 
-mtext(text = "​Variável preditora (x1)", side = 1, line = 3, cex = 1.5) 
-axis(2) 
-mtext(text = "​Variável resposta (y1)", side = 2, line = 3, cex = 1.5, las =0) 
-segments(x0 = x1, y0 = mean(y1), y1 = y1, col = "​white",​ lty = 2, lwd = 2.5) 
-points(x1, y1, pch = 19, col = cores[2], cex = 1.5) 
-abline(h = mean(y1), lty = 2, lwd = 3) 
-## 
-plot(x1, y1, type = "​n",​ axes = FALSE, ann = FALSE, , ylim = range(y1), xlim = range(x1)) 
-rect(par()$usr[1],​ par()$usr[3],​ par()$usr[2],​ par()$usr[4], ​ col = rgb(0, 0, 0, 0.15)) 
-axis(1) 
-mtext(text = "​Variável preditora (x1)", side = 1, line = 3, cex = 1.5) 
-axis(2) 
-mtext(text = "​Variável resposta (y1)", side = 2, line = 3, cex = 1.5, las =0) 
-segments(x0 = x1, y0 = y1, y1 = predict(lmxy01),​ col = "​white",​ lty = 2, lwd = 2.5) 
-points(x1, y1, pch = 19, col = cores[2], cex = 1.5) 
-abline( lmxy01, lty = 1, lwd = 3,  col = cores[1]) 
-## 
-plot(x1, y1, type = "​n",​ axes = FALSE, ann = FALSE, , ylim = range(y1), xlim = range(x1)) 
-rect(par()$usr[1],​ par()$usr[3],​ par()$usr[2],​ par()$usr[4], ​ col = rgb(0, 0, 0, 0.15)) 
-axis(1) 
-mtext(text = "​Variável preditora (x1)", side = 1, line = 3, cex = 1.5) 
-axis(2) 
-mtext(text = "​Variável resposta (y1)", side = 2, line = 3, cex = 1.5, las =0) 
-segments(x0 = x1, y0 = y1, y1 = predict(lmxy01),​ col = "​white",​ lty = 2, lwd = 3.5) 
-points(x1, y1, pch = 19, col = cores[2], cex = 2.5) 
-abline( lmxy01, lty = 1, lwd = 5,  col = cores[1]) 
-abline(h =  mean(y1), lty = 2, lwd = 3)#,  col = cores[1]) 
-predxy <- predict(lmxy01) 
-devExp <​-rep(NA,​ length(y1)) ​ 
-devExp[c(1, 3, 6, 8, 9, 11, 12, 15)] <- predxy[c(1, 3, 6, 8, 9, 11, 12, 15)]  ​ 
-devExp[c(2, 4, 5, 7, 10, 13, 14)] <- y1[c(2, 4, 5, 7, 10, 13, 14)]  
-segments(x0 = x1, y0 = mean(y1), y1 = devExp, col = rgb(0, 1, 0, 0.7), lty = 2, lwd = 3.5) 
-legend(x = 170, y = 20, legend = c("​grande média",​ "​modelo linear",​ "​resíduo",​ "​variação explicada"​),​ lty = c(2, 1, 2, 2), col = c(1, cores[1], "​white",​ "​green"​),​ bty = "​n",​ cex = 1.5, lwd = 3) 
-</​code>​ 
- 
-<WRAP center round box 80%> 
-{{  :​02_tutoriais:​tutorial7:​lmPartVar.png?​500 ​ |}} 
-</​WRAP>​ 
- 
-==== Coeficiente de Determinação (R²) ==== 
- 
-Assim como na **análise de variância** podemos utilizar a partição da variação dos dados para calcular a proporção da variação explicada. Vamos calcular esse valor: 
- 
-<code rsplus> 
-dqTotal <- sum((y1 - mean(y1))^2) 
-dqRes <- sum(residuals(lmxy01)^2) 
-R2 <- (dqTotal - dqRes)/​(dqTotal) 
-R2 
-</​code>​ 
- 
-==== R² Ajustado ==== 
- 
-O R² ajustado está relacionado a uma maior precisão na estimativa do R², que depende do tamanho amostral. Existem vários tipos de ajustes, no caso do ''​summary''​ de um objeto ''​lm''​ a formula é: 
- 
-$$R^{2}_{adj} = 1 - (1- R^2) \frac{n - 1}{n - p - 1}$$ 
- 
-<code rsplus> 
-R2adj <- 1 - (1- R2) * ((length(y1) - 1)/​(length(y1) - 1 - 1)) 
-R2adj 
-</​code>​ 
- 
-Com os R², R² ajustado e o teste da partição da variação da tabela da função ''​anova'',​ fechamos o ''​summary''​ do modelo linear simples. Novamente, reconheça esses valores e interprete os valores: 
- 
-<code rsplus> 
-summary(lmxy01) 
-</​code>​ 
- 
-<​code>​ 
-Call: 
-lm(formula = y1 ~ x1) 
- 
-Residuals: 
-    Min      1Q  Median ​     3Q     ​Max ​ 
--11.449 ​ -3.645 ​  ​1.079 ​  ​2.792 ​  ​7.386 ​ 
- 
-Coefficients:​ 
-            Estimate Std. Error t value Pr(>​|t|) ​   ​ 
-(Intercept) 10.97202 ​   2.81569 ​  ​3.897 ​ 0.00184 **  
-x1           ​0.11855 ​   0.02124 ​  5.582 8.89e-05 *** 
---- 
-Signif. codes: ​ 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
- 
-Residual standard error: 5.28 on 13 degrees of freedom 
-Multiple R-squared: ​ 0.7056,​ Adjusted R-squared: ​ 0.683  
-F-statistic:​ 31.16 on 1 and 13 DF,  p-value: 8.892e-05 
-</​code> ​ 
- 
- 
-===== Comparando Modelos ===== 
- 
-A partição da variância é uma ferramenta poderosa para compararmos modelos aninhados. Ou seja, no caso do modelo mais simples estar contido no modelo mais complexo. De fato, a ''​anova''​ já estava fazendo isso quando ​ 
-utilizamos a partição da variação acima. Vamos montar dois modelos: um com e um sem a variável preditora ''​x1'':​ 
- 
-<code rsplus> 
-lmxy01 <- lm(y1 ~ x1) 
-lmNull <- lm(y1 ~ 1) 
-</​code> ​   ​ 
- 
-Acima, montamos novamente o modelo ''​lmxy01''​ e em seguida criamos outro modelo ''​lmNull''​. Para indicar que não há nenhuma preditora, colocamos o numeral ''​1'',​ que indica que o modelo é definido pela grande média de ''​y1''​. ​ 
- 
-No caso estamos comparando as seguintes representações gráficas: 
- 
-<code rsplus> 
-cores <- c(rgb(1, 0, 0, 0.3), rgb(0, 0, 1, 0.3)) 
-x11(width = 12, height = 6, xpos = 1100) 
-par(mfrow = c(1, 2), mar = c(4, 5, 2, 2), las = 1, cex = 1.5, cex.lab = 1.3, cex.axis = 1) 
-plot(x1, y1, type = "​n",​ axes = FALSE, ann = FALSE, , ylim = range(y1), xlim = range(x1)) 
-rect(par()$usr[1],​ par()$usr[3],​ par()$usr[2],​ par()$usr[4], ​ col = rgb(0, 0, 0, 0.15)) 
-axis(1) 
-mtext(text = "​Variável preditora (x1)", side = 1, line = 3, cex = 1.5 ) 
-axis(2) 
-mtext(text = "​Variável resposta (y1)", side = 2, line = 3, cex = 1.5, las =0) 
-segments(x0 = x1, y0 = mean(y1), y1 = y1, col = "​white",​ lty = 2, lwd = 2.5) 
-points(x1, y1, pch = 19, col = cores[2], cex = 1.5) 
-abline(h = mean(y1), lty = 1, lwd = 3, col = cores[1]) 
-text(100, 42, "​Modelo:​ y1 ~ 1", cex = 1.2) 
-## 
-plot(x1, y1, type = "​n",​ axes = FALSE, ann = FALSE, , ylim = range(y1), xlim = range(x1)) 
-rect(par()$usr[1],​ par()$usr[3],​ par()$usr[2],​ par()$usr[4], ​ col = rgb(0, 0, 0, 0.15)) 
-axis(1) 
-mtext(text = "​Variável preditora (x1)", side = 1, line = 3, cex = 1.5) 
-axis(2) 
-mtext(text = "​Variável resposta (y1)", side = 2, line = 3, cex = 1.5, las =0) 
-segments(x0 = x1, y0 = y1, y1 = predict(lmxy01),​ col = "​white",​ lty = 2, lwd = 2.5) 
-points(x1, y1, pch = 19, col = cores[2], cex = 1.5) 
-abline( lmxy01, lty = 1, lwd = 3,  col = cores[1]) 
-text(100, 42, "​Modelo:​ y1 ~ x1", cex = 1.2) 
-</​code>​ 
- 
- 
-<WRAP center round box 80%> 
-{{  :​02_tutoriais:​tutorial7:​m1mNull.png?​600 ​ |}} 
-</​WRAP>​ 
- 
-==== Comparando Modelos com Partição da Variação ==== 
- 
-Para fazer a comparação com esse método é necessário que os modelos sejam <wrap em>​aninhados</​wrap>​ como já foi dito anteriormente. Como o modelo mais complexo contém o mais simples, a única possibilidade é que o mais complexo explique a mesma ou mais variação que o mais simples, não há como explicar menos, já que ele contém todas as variáveis preditoras que o mais simples. Neste caso, a variação explicada é o quanto o modelo mais complexo se ajusta melhor que o mais simples. Enquanto, o que correspondia à variação total, neste caso, é o que não é explicado pelo modelo mais simples. 
-Veja como funciona com os dois modelos que criamos: 
- 
-<code rsplus> 
-anova(lmNull,​ lmxy01) 
-</​code>​ 
- 
-<​code>​ 
-Analysis of Variance Table 
- 
-Model 1: y1 ~ 1 
-Model 2: y1 ~ x1 
-  Res.Df ​    RSS Df Sum of Sq      F    Pr(>​F) ​   ​ 
-1     14 1230.87 ​                                 ​ 
-2     ​13 ​ 362.37 ​ 1    868.51 31.158 8.892e-05 *** 
---- 
-Signif. codes: ​ 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
-</​code>​ 
- 
- 
-Como no nosso caso, o modelo mais simples é o modelo onde a variação de ''​y1''​ não é explicada por ''​x1'',​ temos a mesma partição do que quando fizemos a ''​anova''​ para o modelo ''​y1 ~ x1''​ isoladamente,​ os elementos apenas estão rearranjados de forma diferente. Reconheça todos os elementos dessa partição da variação e sua correspondência com a tabela de partição anterior. 
- 
-Nesse ponto, é desejável que tenha entendido que a partição da variância de um modelo é correspondente a compará-lo com o modelo nulo, ou seja, quanta variância o modelo é capaz de explicar em relação ao modelo nulo. Esse modelo nulo representa o modelo mais simples representado pelos dados, com a variação total no resíduo, e é representado por apenas um parâmetro, a média da variável resposta.  ​ 
- 
- 
-===== As Lagartas ===== 
- 
-<WRAP center round box 60%> 
-**__Dieta de Lagarta__** 
-  
- 
-{{  :​02_tutoriais:​tutorial7:​lagartaTrans.png?​400 ​ |}}  
- 
- 
-O nosso próximo exercício usa os {{ :​02_tutoriais:​tutorial7:​lagarta.txt |dados de crescimento de lagartas}} submetidas a dietas de folhas com diferentes concentrações de taninos. São apenas duas variáveis, ''​growth'',​ o crescimento da lagarta, e ''​tannins'',​ a concentração de taninos. O objetivo é verificar se há relação entre o crescimento da lagarta e a concentração de taninos da dieta. ​ Baixe o arquivo no seu diretório de trabalho e faça a leitura dos dados. 
- 
-</​WRAP>​ 
- 
-Neste caso temos dois modelos plausíveis,​ um em que o crescimento está relacionado à dieta e o outro onde o crescimento independe da dieta de taninos: 
- 
-<code rsplus> 
-lag <- read.table("​lagarta.txt",​ header = TRUE, sep = "​\t"​) 
-str(lag) 
-summary(lag) 
-lmlag00 <- lm(growth ~ 1, data = lag) 
-summary(lmlag00) 
-lmlag01 <- lm(growth ~ tannin, data = lag) 
-summary(lmlag01) 
-</​code>​ 
- 
-Sempre é bom visualizar os modelos em um gráfico, vamos fazê-los: 
- 
-<code rsplus> 
-par(mfrow = c(1, 2), mar = c(4, 4, 2, 1), las = 1, cex = 1.5, cex.lab = 1.3, cex.axis = 1) 
-plot(NULL, NULL, type = "​n",​ axes = FALSE, ann = FALSE, , ylim = range(lag$growth),​ xlim = range(lag$tannin)) 
-rect(par()$usr[1],​ par()$usr[3],​ par()$usr[2],​ par()$usr[4], ​ col = rgb(0, 0, 0, 0.15)) 
-axis(1) 
-mtext(text = "​Taninos",​ side = 1, line = 2.5, cex = 1.75 ) 
-axis(2) 
-mtext(text = "​Crescimento",​ side = 2, line = 2.5, cex = 1.75, las =0) 
-segments(x0 = lag$tannin, y0 = lag$growth, y1 = predict(lmlag00),​ col = "​white",​ lty = 2, lwd = 2.5) 
-points(lag$tannin,​ lag$growth, pch = 19, col = cores[2], cex = 1.5) 
-abline(lmlag00,​ lty = 1, lwd = 3, col = cores[1]) 
-## 
-plot(NULL, NULL, type = "​n",​ axes = FALSE, ann = FALSE, , ylim = range(lag$growth),​ xlim = range(lag$tannin)) 
-rect(par()$usr[1],​ par()$usr[3],​ par()$usr[2],​ par()$usr[4], ​ col = rgb(0, 0, 0, 0.15)) 
-axis(1) 
-mtext(text = "​Taninos",​ side = 1, line = 2.5, cex = 1.75 ) 
-axis(2) 
-mtext(text = "​Crescimento",​ side = 2, line = 2.5, cex = 1.75, las =0) 
-segments(x0 = lag$tannin, y0 = lag$growth, y1 = predict(lmlag01),​ col = "​white",​ lty = 2, lwd = 2.5) 
-points(lag$tannin,​ lag$growth, pch = 19, col = cores[2], cex = 1.5) 
-abline(lmlag01,​ lty = 1, lwd = 3, col = cores[1]) 
-</​code>​ 
- 
-<WRAP center round box 80%> 
-{{  :​02_tutoriais:​tutorial7:​lag2models.png?​600 ​ |}} 
-</​WRAP>​ 
- 
- 
-==== Comparando Modelos de Dieta de Lagartas ==== 
- 
-Com a anova podemos comparar os dois modelos para ver qual o mais plausível: 
- 
-<code rsplus> 
-anova(lmlag00,​ lmlag01) 
-</​code>​ 
- 
-O teste ''​F''​ desta comparação é similar ao da análise de variância. Neste caso a hipótese estatística nula subjacente é que os modelos não diferem na sua capacidade em explicar a variação dos dados. A hipótese alternativa é que o modelo mais complexo explica mais variação do que o mais simples. Neste caso,  o ''​F''​ e o ''​p-value''​ associado nos permite afirmar que o modelo mais complexo explica melhor, pois a probabilidade de incorrermos em erro com essa proposição é muito pequena. 
- 
-Nesse caso, podemos partir para a interpretação do modelo selecionado:​ 
- 
-<code rsplus> 
-summary(lmlag01) 
-</​code>  ​ 
- 
-E novamente, percebam que a partição do modelo isoladamente é a mesma partição de variância da comparação com o modelo nulo, ou aquele que é representado apenas pela média da variável resposta: 
- 
- 
-<code rsplus> 
-anova(lmlag01) 
-</​code>  ​ 
- 
-==== Diagnóstico do Modelo ==== 
- 
-Ao encontrar um modelo adequado para representar os dados, antes de interpretá-lo,​ precisamos fazer o diagnósticos das premissas do modelo. Neste caso simples, podemos fazer isso apenas pelos gráficos padrão do R. São quatro gráficos básicos que nos permitem diagnosticar:​ (1) se os residuos desviam de uma distribuição normal, (2) se a relação entre as variáveis apresenta alguma curvatura, (3) dados influentes e limítrofes. 
- 
-<code rsplus> ​ 
-par(mfrow = c(2, 2)) 
-plot(lmlag01) 
-</​code>​ 
- 
-Não está no escopo deste curso ensinar esses gráficos diagnósticos. Para isso veja o tutorial de [[http://​labtrop.ib.usp.br/​doku.php?​id=cursos:​planeco:​roteiro:​07a-clasrcmdr| Regressão]] da disciplina [[http://​labtrop.ib.usp.br/​doku.php?​id=cursos:​planeco| Princípios de Planejamento e Análise de Dados ]] 
-===== Variável Indicadora ===== 
- 
-Nos modelos lineares as variáveis contínuas e categóricas são colocadas como preditoras, indistintamente. Entretanto, no modelo elas são tratadas diferentemente. As variáveis categóricas são transformadas em <wrap em>​variáveis indicadoras ou dummy</​wrap>​. Entender essa transformação é importante para a interpretação dos resultados do modelo. Vamos aqui analisar o caso bem simples da produtividade de culturas em diferentes solos, os dados que utilizamos no tutorial de [[02_tutoriais:​tutorial6b:​start|Partição de Variação dos Dados]]. Vamos então refazer os primeiros passos de leitura dos dados: 
- 
-<code rsplus> 
-are <- c(6,​10,​8,​6,​14,​17,​ 9, 11, 7, 11) 
-arg <-  c(17, 15, 3, 11, 14, 12, 12, 8, 10, 13) 
-hum <- c(13, 16, 9, 12, 15, 16, 17, 13, 18, 14) 
-solo <- rep(c("​arenoso",​ "​argiloso",​ "​humico"​),​ each = 10) 
-cultivar <- data.frame(producao = c(are, arg, hum), solo = solo) 
-str(cultivar) 
-</​code>​ 
- 
-Vamos fazer uma inspeção nos dados na forma de um ''​boxplot'':​ 
- 
-<​code>​ 
-cols <- c(rgb(0, 0, 0, 0.1), rgb(1, 0,0, 0.3), rgb(0,0,1, 0.3))    
-par(mar = c(4,4,2,1), las = 1, cex = 1.5) 
-boxplot(producao ~ solo, data = cultivar, col = cols, xlab = "Tipo de solos",​ ylab = "​Produção (ton/​ha)",​ range = 0) 
-points(x = jitter(rep(1:​3,​ each = 10)), y = cultivar$producao,​ bg = rep(cols, each = 10), pch = 21) 
-</​code>​ 
- 
-<WRAP center round box 70%> 
- 
-{{  :​02_tutoriais:​tutorial6:​cultFig.png?​600 ​ |}} 
- 
-</​WRAP>​ 
- 
- 
-Podemos construir o modelo com esses dados da maneira convencional no ''​lm'':​ 
- 
-<code rsplus> 
-lmcult <- lm(producao ~ solo, data = cultivar) 
-summary(lmcult) 
-</​code>​ 
- 
-O primeiro passo é reconhecer no resumo do modelo os mesmos resultados do teste de **análise de variáncia** que fizemos no tutorial anterior. Apenas com um pouco mais de informação associado aos coeficientes do modelo e as imprecisões nas suas estimativas. 
- 
-<​code>​ 
-Call: 
-lm(formula = producao ~ solo, data = cultivar) 
- 
-Residuals: 
-   ​Min ​    1Q Median ​    ​3Q ​   Max  
-  -8.5   ​-1.8 ​   0.3    1.7    7.1  
- 
-Coefficients:​ 
-             ​Estimate Std. Error t value Pr(>​|t|) ​   ​ 
-(Intercept) ​    ​9.900 ​     1.081   9.158 9.04e-10 *** 
-soloargiloso ​   1.600      1.529   ​1.047 ​ 0.30456 ​   ​ 
-solohumico ​     4.400      1.529   ​2.878 ​ 0.00773 **  
---- 
-Signif. codes: ​ 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
- 
-Residual standard error: 3.418 on 27 degrees of freedom 
-Multiple R-squared: ​ 0.2392,​ Adjusted R-squared: ​ 0.1829 ​ 
-F-statistic:​ 4.245 on 2 and 27 DF,  p-value: 0.02495 
-</​code>​ 
- 
-Outro ponto interessante no resumo é que, apesar de termos apenas uma variável preditora, tipo de ''​solo'',​ temos três coeficientes estimados, sendo um deles um intercepto, que não faz muito sentido para variáveis preditoras categóricas. 
- 
-Qual o significado desses coeficientes?​ Vamos ver a representação gráfica deste modelo, que já fizemos anteriormente:​ 
- 
-<code rsplus> 
-mSolos <- tapply(cultivar$producao,​ cultivar$solo,​ mean) 
-mSolosVetor <- rep(mSolos, each = 10) 
-colvector <- rep(cols, each = 10) 
-par(mar = c(4,4,2,1), las = 1, cex = 1.5) 
-plot(x = 1:30, y = cultivar$producao , ylim = c(0,20), xlim = c(0, 30), pch=(rep(c(15,​16,​17),​ each=10)), col = colvector, ylab = "​Produtividade (ton/​ha)",​ xlab = "​Observações",​ cex = 1.5) 
-segments(x0 = 1:30, y0 = cultivar$producao,​ y1= mSolosVetor,​ col = colvector, lwd = 2) 
-segments(x0 = c(1,11, 21), y0 = mSolos, x1 = c(10, 20, 30), col = cols, lwd =1.5) 
-legend("​bottomright",​ legend = c("​arenoso",​ "​argiloso",​ "​humico"​),​ pch = 15:17 ,col = cols, title = "​Solos",​ bty = "​n"​) 
-</​code>​ 
- 
-<WRAP center round box 80%> 
-{{  :​02_tutoriais:​tutorial6b:​resdSolo.png?​600 ​ |}} 
-</​WRAP>​ 
- 
- 
-Nessa representação,​ o modelo é representado pelas três linhas horizontais,​ que representam as médias de cada tipo de solo, que calculamos acima no objeto ''​mSolos''​. Então os coeficientes são as estimativas da produtividade média dos solos? Quase isso! 
- 
-<code rsplus> 
-tapply(cultivar$producao,​ cultivar$solo,​ mean) 
-coef(lmcult) 
-</​code>​ 
- 
-Os coeficientes são a diferenças entres as médias do nível em relação ao ''​intercepto'',​ que no caso das variáveis categóricas é sempre o primeiro nível. Caso os níveis não tenham sido colocados na ordem utilizando a função ''​factor''​ o ''​lm''​ automaticamente ordena os níveis em ordem alfabética e coloca no intercepto o primeiro deles. Portanto, uma variável preditora categórica gera o número de níveis de coeficientes,​ sendo que um deles é deslocado para o intercepto. 
- 
-Veja as soma dos coeficientes com o intercepto: 
- 
-<code rsplus> 
-tapply(cultivar$producao,​ cultivar$solo,​ mean) 
-coef(lmcult) 
-coef(lmcult)[1] + coef(lmcult)[2:​3] 
-</​code>​ 
- 
-Vamos demonstrar o que acontece dentro do modelo fazendo a transformação de variável indicadora, antes de criar o modelo. A partir da variável solo criamos duas outras, ''​arg''​ e ''​hum'',​ que indicam qual solo pertence cada observação. As variáveis indicadoras tem ''​1'',​ se a observação pertence ao solo em questão, ou ''​0'',​ caso pertença a outro tipo de solo: 
- 
-<code rsplus> 
-cultivar$arg <- rep(0, 30) 
-cultivar$hum <- rep(0, 30) 
-cultivar$arg[cultivar$solo == "​argiloso"​] <- 1 
-cultivar$arg 
-cultivar$hum[cultivar$solo == "​humico"​] <- 1 
-cultivar$hum 
-head(cultivar) 
-cultivar[11:​16,​] 
-tail(cultivar) 
-</​code> ​ 
- 
-Agora, vamos construir nosso modelo com as variáveis indicadoras:​ 
- 
-<code rsplus> 
-lmcultInd <- lm(producao ~ arg + hum, data = cultivar) 
-summary(lmcultInd) 
-</​code>​ 
- 
-Reconheça o resultado idêntico dos modelos, ​ veja o que acontece se comparamos os modelos por ''​anova'':​ 
- 
-<code rsplus> 
-anova(lmcultInd,​ lmcult) 
-</​code>​ 
- 
-<WRAP center round box 90%> 
-**__O que devo ter entendido até aqui:__** 
-<WRAP center round tip 100%> 
- 
-  * 1. Todos os valores apresentados no ''​summary''​ de um modelo. 
-  * 2. A partição da variação no ''​lm''​. 
-  * 3. O que a ''​anova''​ apresenta ao comparar dois modelos. 
-  * 4. Como as variáveis categóricas são tratadas nos modelos no ''​lm''​. 
-</​WRAP>​ 
- 
-São muitos conceitos não triviais. Em um primeiro momento, busque alguma compreensão,​ mesmo que ainda pareça nebuloso. 
- 
-</​WRAP>​ 
- 
-<WRAP center round box 100%> 
-Aula síncrona da disciplina no google meet, gravada em 01 de outubro de 2020. Nela abordo a construção e interpretação de modelos lineares simples no ambiente de programação R, focando no resumo (''​summary''​) com as principais informações do modelo. Uma bom entendimento do resumo do modelo é essencial para interpretação correta do resultado. Veja curso completo em: 
-http://​ecor.ib.usp.br 
-<WRAP center round tip 80%> 
-{{youtube>​VRrJ487k5qY}} 
-</​WRAP>​ 
- 
-</​WRAP>​ 
- 
- 
- 
-<WRAP center round box 50%> 
-Siga para: 
-  * Capítulo da apostila [[03_apostila:​06-modelos|]]:​ do início até o tópico [[03_apostila:​06-modelos#​variavel_factor_como_variavel_indicadora_dummy| Variável fator como variável indicadora...]] 
-  * Exercícios [[01_curso_atual:​exercicios7|]] 
-</​WRAP>​ 
  
02_tutoriais/tutorial7/start.txt · Última modificação: 2023/09/11 15:58 (edição externa)