\documentclass[a4paper, oneside, 10pt]{article} \usepackage[english]{babel} \usepackage[unicode]{hyperref} \usepackage[utf8x]{inputenc} \usepackage{listings} \usepackage{graphicx} \graphicspath{{media/}} \date{\today} \title{} \author{} \begin{document} \begin{itemize} \item \href{http://ecor.ib.usp.br/doku.php?id=02_tutoriais:tutorial6b:start}{Tutorial} \item \href{http://ecor.ib.usp.br/doku.php?id=01_curso_atual:exercicios6b}{ Exercícios} \item Apostila \end{itemize} \section{\texorpdfstring{6b. Partição da Variação dos Dados}{6b Partiao da Variaao dos Dados}} \label{sec:6b_partiao_da_variaao_dos_dados} \subparagraph{\texorpdfstring{Vídeo gravado pelo Google Meet em aula síncrona no dia 30 de setembro de 2020. Sem edição.}{Video gravado pelo Google Meet em aula sincrona no dia 30 de setembro de 2020 Sem ediao}} \label{sec:video_gravado_pelo_google_meet_em_aula_sincrona_no_dia_30_de_setembro_de_2020_sem_ediao} \centering\includegraphics[keepaspectratio=true,width=0.8\textwidth]{02_tutoriais/tutorial6b/phdStatistics} O teste t, apresentado no \href{http://ecor.ib.usp.br/doku.php?id=02_tutoriais:tutorial6:start}{tutorial 6a}, é usado apenas para o caso de termos uma variável resposta numérica contínua e uma preditora categórica com \textbf{dois níveis}. Caso a preditora tenha mais do que dois níveis, precisamos usar um outro teste que é uma generalização do teste t, o teste de Análise de Variância ou ANOVA. O teste está baseado no princípio de partição da variação dos dados. A variação total dos dados é particionada nos componentes do que é explicado e aquele que não é explicado pela variável preditora categórica. Esse conceito é aplicado de maneira mais ampla na estatística, utilizado em outros tipos de estatística e para a tomada de decisão do modelo que melhor explica a variação nos dados. Por isso, vamos focar este tutorial no conceito da partição da variação. Para exemplificar a partição da variância associada à ANOVA, vamos usar o exemplo de dados de colheita de um cultivar em diferentes tipos de solos, apresentado no livro de Robert Crawley, \href{http://www.bio.ic.ac.uk/research/mjcraw/therbook/index.htm/}{The R Book}, como segue abaixo: \textbf{\underline{Tradução livre da descrição do livro ,,\emph{The R Book}" (\href{http://www.bio.ic.ac.uk/research/mjcraw/therbook/index.htm/}{Crawley, 2007}) }} \raggedleft\includegraphics[keepaspectratio=true,width=0.8\textwidth]{02_tutoriais/tutorial6/crawley} ,,... a melhor forma de entender o que está acontecendo é trabalharmos um exemplo. Temos um experimento em que a produção agrícola, por unidade de área, é medida em 10 campos de cultivo, selecionados aleatoriamente em cada um de três tipos diferentes de solo. Todos os campos foram semeados com a mesma variedade de semente e manejados com as mesmas técnicas (fertilizantes, controle de pragas). O objetivo é verificar se o tipo de solo afeta significativamente o rendimento de culturas, e caso afete, quanto."\footnote{The best way to see what is happening is to work through a simple example. We have an experiment in which crop yields per unit area were measured from 10 randomly selected fields on each of three soil types. All fields were sown with the same variety of seed and provided with the same fertilizer and pest control inputs. The question is whether soil type significantly affects crop yield, and if so, to what extent.} Vamos organizar os dados apresentados no livro em um objeto no R diretamente:\lstset{frame=single, language=rsplus} \begin{lstlisting} 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) \end{lstlisting} Primeiro vamos fazer uma inspeção nos dados na forma de um \texttt{boxplot}:\lstset{frame=single} \begin{lstlisting} 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) \end{lstlisting} \centering\includegraphics[keepaspectratio=true,width=0.8\textwidth]{02_tutoriais/tutorial6/cultFig} A pergunta aqui é: Há variação na produtividade média entre solos? \subsection{\texorpdfstring{Variação Total dos Dados}{Variaao Total dos Dados}} \label{sec:variaao_total_dos_dados} A partição da variação inicia-se pelo reconhecimento e cálculo da variação total associada aos dados. A variação total dos dados é baseada nos desvios das observações em relação à grande média, que no caso, é a média de todos os campos de cultivo. Vamos representar esses desvios:\lstset{frame=single, language=rsplus} \begin{lstlisting} par(mar = c(4,4,2,1), las = 1, cex = 1.5) colvector <- rep(cols, each= 10) 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 = "Variável Resposta", xlab = "Observações", cex = 1.5) for(i in 1:30) { lines(c(i,i),c(cultivar$producao[i], mean(cultivar$producao)), col = colvector[i]) } abline(h = mean(cultivar$producao), lty = 2) legend("bottomright", legend = c("arenoso", "argiloso", "humico"), pch = 15:17 ,col = cols, title = "Solos", bty = "n") \end{lstlisting} \textbf{\underline{ Representação dos desvios dos dados:} } em relação à média geral. \includegraphics[keepaspectratio=true,width=0.8\textwidth]{02_tutoriais/tutorial6b/desCult} No gráfico esta variação é representada pelos segmentos verticais coloridos. A grande média é definida como a média de produtividade de todos os campos de cultivo (n=30), independente do tipo de solo, e é representada pela linha preta horizontal tracejada. Medimos essa variação total pela \texttt{soma quadrática}: os valores dos desvios dos dados em relação à grande média (segmentos verticais no gráfico) elevados ao quadrado e posteriormente somados. Essa soma quadrática total é nossa medida de variação. $$ SQ_{"total"} = \sum_{i=1}^k\sum_{j=1}^n (y_{ij} - \bar{\bar{y}})^2 $$ \\ Aqui calculamos o valor de somatória quadrática que é a base da análise de ANOVA.\lstset{frame=single, language=rsplus} \begin{lstlisting} (mGeral <- mean(cultivar$producao)) (dGeral <- cultivar$producao - mGeral) (dqGeral <- dGeral^2) (sqGeral <- dqGeral) (sdqGeral <- sum(sqGeral)) \end{lstlisting} Fizemos acima todos os passos isoladamente, pois iremos utilizar mais à frente alguns desses valores intermediários. Vamos iniciar a construção da nossa tabela de ANOVA, incluindo a medida de variação total na sua posição: \centering\includegraphics[keepaspectratio=true,width=0.8\textwidth]{02_tutoriais/anova1} \subsection{\texorpdfstring{Particionando a Variação}{Particionando a Variaao}} \label{sec:particionando_a_variaao} Um ponto importante da soma quadrática é que esta variação é aditiva e pode ser decomposta. Uma parte dessa variação é explicada pelo tratamento, que são nossas variáveis preditoras \footnote{podemos ter mais de uma}. A porção não explicada pelas variáveis preditoras é o resíduo, algumas vezes chamado de \emph{erro}. A porção não explicada está associada à variação aleatória dos dados ou a algum, ou vários fatores que não foram incluídos ou controlados no nosso experimento. \subsubsection{\texorpdfstring{Variação Não Explicada}{Variaao Nao Explicada}} \label{sec:variaao_nao_explicada} Vamos primeiro representar essa variação não explicada, ou variação interna aos níveis da variável categórica:\lstset{frame=single, language=rsplus} \begin{lstlisting} (mSolos <- tapply(cultivar$producao, cultivar$solo, mean)) mSolosVetor <- rep(mSolos, 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") \end{lstlisting} No código acima utilizamos a função \texttt{segments} e não foi necessário utilizar a iteração. Além disso, criamos um vetor de médias com a repetição das médias dos respectivos solos para cada observação no \texttt{mSolosVetor}, para facilitar a construção do gráfico. \textbf{\underline{ Representação dos resíduos dos modelo:} } desvios dos dados em relação às médias dos tratamentos. \includegraphics[keepaspectratio=true,width=0.8\textwidth]{02_tutoriais/tutorial6b/resdSolo} Para calcular a variação não explicada precisamos usar o mesmo procedimento da variação total: elevar os resíduos ao quadrado e somar, resultando também em uma somatória quadrática.\lstset{frame=single, language=rsplus} \begin{lstlisting} (sqIntra <- sum((cultivar$producao - mSolosVetor)^2)) \end{lstlisting} Incluindo esse valor na nossa tabela: \centering\includegraphics[keepaspectratio=true,width=0.8\textwidth]{02_tutoriais/anova3} \subsubsection{\texorpdfstring{Variação Explicada}{Variaao Explicada}} \label{sec:variaao_explicada} A variação explicada pelas variáveis preditoras é a diferença entre a variação total e a não explicada, devido à característica aditiva das somatórias quadráticas. Vamos representar a variação que foi explicada, ou variação entre os níveis da variável preditora categórica, em um gráfico:\lstset{frame=single, language=rsplus} \begin{lstlisting} 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(0, 1 ,2), each=10)), col = colvector, ylab = "Produtividade (ton/ha)", xlab = "Observações", cex = 1) points(x = 1:30, y = mSolosVetor, pch = rep(c(15,16,17), each=10), col = colvector, cex = 1.5) segments(x0 = 1, y0 = mGeral, x1= 30, col = 1, lty = 2, lwd = 1.5) segments(x0 = 1:30, y0 = mSolosVetor, y1 = rep(mGeral, 30), col = colvector, lwd =1.5) legend("bottomright", legend = c("arenoso", "argiloso", "humico"), pch = 15:17 ,col = cols, title = "Solos", bty = "n") \end{lstlisting} \textbf{\underline{ Representação da variação explicada dos dados} } \includegraphics[keepaspectratio=true,width=0.8\textwidth]{02_tutoriais/tutorial6b/soloExplica} Vamos calcular e incluir na nossa tabela essa variação, calculada pelo soma quadrática dos segmentos representados na figura:\lstset{frame=single, language=rsplus} \begin{lstlisting} (sqEntre <- sum((mSolosVetor - mGeral)^2)) \end{lstlisting} Incluindo esse valor na nossa tabela: \centering\includegraphics[keepaspectratio=true,width=0.8\textwidth]{02_tutoriais/tutorial6b/tabAnova04} \subsubsection{\texorpdfstring{Graus de Liberdade}{Graus de Liberdade}} \label{sec:graus_de_liberdade} Precisamos agora calcular os \textbf{desvios quadráticos médios} que são as somas quadráticas dividido pelos graus de liberdade (\texttt{gl}). Vamos utilizar a soma quadrática total para exemplificar o cálculo dos graus de liberdade: para calcular essa soma utilizamos os \texttt{30} valores de produtividade e calculamos a média geral. No caso, por termos calculado a média, perdemos um grau de liberdade e ficamos com \texttt{29 gl}. Utilizando a mesma lógica para a soma quadrática não explicada, partimos das mesmas \texttt{30} informações e calculamos as \texttt{3} médias de produtividade dos solos, portanto, ficamos com \texttt{27 gl}. Por fim, na soma quadrática explicada temos \texttt{3} informações, as médias de cada solo, e calculamos uma estatística, a média geral, ficando com \texttt{2} graus de liberdade. Com esses informação podemos então, calcular esses desvios quadráticos médios:\lstset{frame=single, language=rsplus} \begin{lstlisting} (sq <- c(sdqGeral, sqIntra, sqEntre)) glib <- c(29, 27, 2) (msq <- sq/glib) \end{lstlisting} Agora nossa tabela está quase completa, se calcularmos a razão entre os desvios médios explicado pelo não explicado, temos o cálculo da estatística da ANOVA, o F-Fisher:\lstset{frame=single, language=rsplus} \begin{lstlisting} (fcultiva <- msq[3]/msq[2]) \end{lstlisting} \centering\includegraphics[keepaspectratio=true,width=0.8\textwidth]{02_tutoriais/anova5} Só falta agora o cálculo do p-valor associado ao \href{https://pt.wikipedia.org/wiki/Teste_F}{teste F} e à estatística F. O F-Fisher é uma distribuição probabilística que tem dois parâmetros: os graus de liberdade dos cálculos da (1) variação média entre e (2) intra grupos. \lstset{frame=single, language=rsplus} \begin{lstlisting} pf(q = fcultiva, df1 = 2, df2 = 27, lower.tail = FALSE) \end{lstlisting} Pronto! Nosso almejado p-valor e tabela completa! \centering\includegraphics[keepaspectratio=true,width=0.8\textwidth]{02_tutoriais/tutorial6b/tabAnova06} \paragraph{\texorpdfstring{Distribuição de F}{Distribuiao de F}} \label{sec:distribuiao_de_f} Os gráficos de outras aulas apresentaram a distribuição de densidade probabilística, onde a variável \texttt{y} é relacionada à probabilidade de cada valor em intervalos muito pequenos. O valor da probabilidade cumulativa é a área da curva até o valor fornecido, o que é retornado pela função \texttt{pf}. No caso, como utilizamos o argumento \texttt{lower.tail = FALSE}, a função retorna a outra área da curva, representada pela figura a seguir:\lstset{frame=single, language=rsplus} \begin{lstlisting} curve(expr=df(x, 2,27), main="Distribuição F de Fisher (df=2,27)", xlab="Valor FALSE",ylab="Densidade Probabilística (df)", xlim=c(0,10)) abline(v = fcultiva, col="red") abline(h = 0, lty = 2) xf <- seq (fcultiva, 10, 0.01) ydf <- df(xf, 2, 27) polygon(c(fcultiva, xf), c(0, ydf), col="red") text(6, 0.1,paste("pf(x) =",round(pf(fcultiva,2,27,lower.tail=F),4)), cex = 1.2, col="red") \end{lstlisting} \centering\includegraphics[keepaspectratio=true,width=0.8\textwidth]{02_tutoriais/tutorial6b/pvalor} \subsubsection{\texorpdfstring{Anova no R}{Anova no R}} \label{sec:anova_no_r} A tabela de anova é tão importante que a função \texttt{anova} no R retorna a tabela para objetos de modelos. Vamos ver isso no próximo tópico. Para o teste da análise de variância, propriamente dito, a função utilizada é \texttt{aov}. No caso do objeto produzido pela função \texttt{aov} a tabela aparece se aplicarmos a função \texttt{anova} ou \texttt{summary} ao objeto \texttt{aov}, que também é um objeto da classe modelo linear \texttt{lm}.\lstset{frame=single, language=rsplus} \begin{lstlisting} (cultAov <- aov(producao ~ solo, data = cultivar)) class(cultAov) anova(cultAov) summary(cultAov) \end{lstlisting} \end{document}