Diferenças

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

--- 02_tutoriais:tutorial7:start [2020/10/01 15:24]
adalardo [Variável Indicadora]
+++ 02_tutoriais:tutorial7:start [2020/10/02 18:26]
adalardo [7a. Regressão Linear Simples]
@@ Linha 7: / Linha 7: @@
 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 =====
@@ Linha 12: / Linha 16: @@
 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 40%>
+<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 do trabalho é fundamental!
+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!
@@ Linha 47: / Linha 54: @@
-Parece complicado, mas é simples gerar dados aleatórios com essa estrutura do R. Vamos definir primeiro qual 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$.
+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 proveniente 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:
+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)
-x <- round(seq(12, 220, len = 15), 1)
+x1 <- round(seq(12, 220, len = 15), 1)
-y0 <- 10.3 + 0.12 * x
+y0 <- 10.3 + 0.12 * x1
-res <- rnorm(length(x), 0, 5)
+res <- rnorm(length(x1), 0, 5)
 y1 <- y0 + res
+xm <- mean(x1)
+ym <- mean(y1)
 </code>
@@ Linha 67: / Linha 76: @@
 <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))
+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)
@@ Linha 74: / Linha 83: @@
 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(x, y0, pch = 16, cex = 0.8, col = cores[1] )
+points(x1, y0, pch = 16, cex = 0.8, col = cores[1] )
-points(x, y1, pch = 19, col = cores[2])
+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>
@@ Linha 84: / Linha 93: @@
 ==== 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''.
@@ Linha 119: / Linha 130: @@
 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.
+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.
@@ Linha 289: / Linha 319: @@
   *1. Cria valores aleatórios de uma normal com média zero e desvio padrão igual a 5;
-  *2. Somo esses valores ao valor ''y0 '', proveniente da relação ''y ~ x'' da população e cria os dados ''ySim'';
+  *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'';
@@ Linha 403: / Linha 433: @@
 ===== Tabela de Anova de uma Regressão =====
-<WRAP center round box 60%>
+<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>
@@ Linha 691: / Linha 725: @@
 ===== 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]]:
+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>

ecoR

Ferramentas do usuário

Ferramentas do site

Diferenças

Ferramentas da página