* [[bie5782:02_tutoriais:tutorial7:start|Tutorial]] * [[bie5782:01_curso_atual:exercicios7| Exercícios]] * [[bie5782:03_apostila:06-modelos| Apostila]] ====== 7a. Regressão Linear Simples ====== ===== Ajuste e diagnóstico de uma regressão ===== Vamos usar o objeto ''Animals'', do pacote //MASS// library(MASS) data(Animals) str(Animals) Nele estão as massas corporais e do cérebro de 28 espécies de vertebrados, em gramas. Vamos explorar a relação entre tamanho do cérebro e do corpo em um gráfico de dispersão: plot(brain~body, data=Animals) A relação claramente não é linear, então podemos tentar em escala //log//: plot(brain~body,data=Animals, log="xy") Ou transformar os dados em seus logaritmos, dentro da função ''plot'': plot(log(brain)~log(body), data=Animals) Agora a relação parece linear, mas há 3 pontos bem discrepantes, de animais muito grandes mas com cérebros pequenos. Quem são? Animals[log(Animals$body) > 8 & log(Animals$brain) < 6, ] São os dinossauros! Ainda assim vamos ajustar uma regressão linear com todos os dados, transformados para seu logaritmos. Note que o argumento que damos para a função ''lm'' é o mesmo que demos para a função ''plot'': anim.m1 <- lm(log(brain)~log(body),data=Animals) E vamos adicionar a linha de regressão ao gráfico: abline(anim.m1, col="blue") Antes de prosserguir, vamos inspecionar os gráficos de diagnóstico (um a um, aperte ''Enter'' para passar para o próximo): plot(anim.m1) Há vários problemas. Você os identifica? Vamos agora ajustar uma regressão sem os dinossauros, retirados com o argumento ''subset'': anim.m2 <- lm(log(brain)~log(body),data=Animals, subset=!(log(Animals$body)>8&log(Animals$brain)<6)) Um gráfico com as retas dos dois modelos, para comparação; veja como a exclusão destes três pontos mudou bastante a inclinação da reta: plot(log(brain)~log(body), data=Animals) abline(anim.m1, col="blue", lty=2) abline(anim.m2, col="red") Novos Diagnósticos (4 gráficos de uma vez) mostram que agora ainda temos alguns problemas (quais?): par(mfrow=c(2,2)) plot(anim.m2) par(mfrow=c(1,1)) Mas vamos prosseguir. Fazemos a ANOVA e o resumo do modelo: anova(anim.m2) summary(anim.m2) Calculamos os coeficientes estimados e os intervalos de confiança: coef(anim.m2) confint(anim.m2) Como a relação entre os logaritmos dos valores é linear, na escala original a relação é de uma equação de potência (//power law//), muito usada para expressar relações alométricas: $$y \ = \ Cx^b $$ Onde $b =$ inclinação e $C = e^a$, sendo $a$ o intercepto na escala log-log, se usamos logaritmos naturais. Agora vamos acrescentar à planilha os valores previstos, já na escala original: cf <- as.numeric(coef(anim.m2)) Animals$prev <- round(exp(cf[1])*Animals$body^cf[2],1) Animals Para quais animais a estimativa foi melhor? E para quais foi pior? Você consegue relacionar isto com os diagnósticos que fizemos acima? ===== Simulação com dados Não Normais ===== Neste tutorial vamos investigar a robustez do modelo de regressão linear simples (Gaussiana) à violação de seus pressupostos. Para isto, vamos criar uma variável resposta cujo valor médio é uma função linear de uma variável preditora, mas com uma distribuição muito diferente da Normal. Antes de mais nada, defina uma semente de números aleatórios, para reproduzir a mesma simulação, e verificar os resultados relatados abaixo. set.seed(42) Vamos então simular 30 valores da variável preditora, sorteados de uma distribuição uniforme: x <- runif(30,0,10) Em seguida criamos uma variável resposta, vinda de uma população com [[http://en.wikipedia.org/wiki/Exponential_distribution|distribuição exponencial]], cuja a média equivale a (2 + 5 * preditora). A distribuição exponencial tem apenas um parâmetro, λ, que corresponde ao inverso da média. y <- rexp(30,rate=1/(2+5*x)) Inspecionando o gráfico de dispersão: plot(y~x) Ajustando a regressão linear: m1 <- lm(y~x) Os gráficos de diagnóstico mostram desvios importantes em relação às premissas do modelo de regressão, como esperado: os resíduos não são normais e não têm variância constante. Há alguns pontos com alavancagem (//leverage//) forte. par(mfrow=c(2,2)) plot(m1) par(mfrow=c(1,1)) Quais os coeficientes estimados pela regressão? coef(m1) São bem diferentes dos valores reais de intercepto (a = 2) e inclinação (b = 5). Será que este padrão se mantém? Podemos investigar isto repetindo várias vezes a tomada de uma amostra de y e o ajuste do modelo, e para cada uma delas guardar os valores dos coeficientes. Vamos fazer isso, e também guardar os limites dos intervalos de confiança a 5% dos coeficientes, para verificar quantas vezes eles de fato incluem os valores verdadeiros. Primeiro criamos uma matriz para guardar os resultados. Cada linha guarda o resultado de uma aleatorização, e deixamos seis colunas, para as estimativas do intercepto e inclinação, e para os limites de seus intervalos de confiança: result <- matrix(ncol=6,nrow=2000) Em seguida usamos um //loop// (função ''for'') para repetir as simulações 2000 vezes. A cada vez guardamos os coeficientes e limites de seus intervalos de confiança: for(i in 1:2000){ yr <- rexp(30,rate=1/(2+5*x)) m <- lm(yr~x) CIs <- confint(m) result[i,1:2] <- as.numeric(coef(m))# guarda os coeficientes nas duas primeiras colunas result[i,3:4] <- CIs[1,]#guarda limites dos ICs do intercepto result[i,5:6] <- CIs[2,]#guarda limites dos ICs da inclinacao } Agora vamos verificar a distribuição dos valores dos interceptos estimados, que estão guardadas na primeira coluna da matriz de resultados: hist(result[,1]) Acrescentamos a média destes valores: abline(v=mean(result[,1]), lty=2, col="blue") E o valor real de a = 2 abline(v=2, col="red", lty=2) A diferença é minúscula, o que indica que o valor esperado das estimativas (i.e. sua média ou esperança estatística), é muito próximo do valor do parâmetro. Agora vamos verificar o mesmo para as inclinações estimadas, que guardamos na segunda coluna da matriz de resultados: hist(result[,2]) abline(v=mean(result[,2]), lty=2, col="blue") abline(v=5, col="red", lty=2) Nos dois casos o viés das estimativas parece ser muito pequeno. Vamos calcular este viés em percentual em relação ao valor real para o intercepto: 100*(2-mean(result[,1]))/2 ## 7.8% e para a inclinação: 100*(5-mean(result[,2]))/5 ## -0.3% Agora vamos verificar as estimativas pelo intervalo de confiança. Quantas vezes os intervalos estimados incluíram o valor real? Como o intervalo foi calculado a 5%, esperaríamos que dos 2000 intervalos calculados, 1900 incluíssem o valor real. Vamos fazer esta contagem. Como guardamos os limites do intervalo de confiança do intercepto nas colunas 3 e 4 da matriz de resultados, podemos fazer a contagem como a soma de um teste lógico: sum(result[,3]<=2&result[,4]>=2) Em 1995 das 2000 simulações o intervalo incluiu a média, ou seja, a significância do intervalo é de 5/2000 = 0,0025, o que é menor do que a significância nominal de 0,05. Para a inclinação temos: sum(result[,5]<=5&result[,6]>=5) Aqui as simulações mostram que o intervalo de confiança teve uma significância um pouco maior que a nominal (0,055). Em resumo, surpreendentemente, as estimativas dos coeficientes da regressão por ponto (valores estimados) e por intervalo (intervalos de confiança) mostraram-se muito robustas a este caso de violação das premissas. No entanto, **não generalize este resultado**, pois não sabemos como variações desta simulação se comportariam. ===== Ajuste de Polinômios ===== Galileu Galilei investigou a distância percorrida por um corpo lançado de uma rampa a diferentes alturas. Para isto ele usou um plano inclinado a uma certa altura do chão, de onde soltou bolas de metal. Com este experimento ele verificou que a relação entre a altura de lançamento e a distância percorrida não é linear. Mais detalhes na ajuda do objeto [[http://finzi.psych.upenn.edu/R/library/UsingR/html/galileo.html|galileo]], no pacote [[http://finzi.psych.upenn.edu/R/library/UsingR/html/00Index.html#U|UsingR]]. Este tutorial é retirado do livro companheiro deste pacote((John Verzani. Using R for Introductory Statistics. Chapman & Hall/CRC, Boca Raton, FL, 2005.)). As alturas de lançamentos e as distâncias percorridas obtidas por Galileu são((as duas medidas estão em //puntos//, uma unidade que corresponde a 169/180 mm)): init.h = c(600, 700, 800, 950, 1100, 1300, 1500) h.d = c(253, 337, 395, 451, 495, 534, 573) Vamos fazer o gráfico de dispersão: plot(h.d~init.h) Galileu tinha razões teóricas para supor que a relação seria uma parábola, e não uma reta, como propunha o conhecimento vigente. Hoje podemos tratar isto como um problema simples de seleção de modelos, em que os modelos concorrentes são uma equação da reta: $$y \ = \ a+bx$$ e uma equação da parábola, que é um polinômio de segundo grau: $$y \ = \ a+bx+cx^2$$ O primeiro modelo é um caso particular do segundo, no qual o coeficiente $$c$$ é zero. Portanto, podemos ajustar os dois modelos e usar o comando ''anova'' para testar se o acréscimo do termo quadrático implicou em melhora. O ajuste da reta é feito do modo usual: mod1 <- lm(h.d~init.h) Mas para ajustar a parábola corretamente usamos a função ''I()'' (de //insulate//), necessária para isolar a operação matemática de elevar ao quadrado. Sem isto, o R entenderá o termo ''init.h^2'' como uma interação, detalhes [[bie5782:03_apostila:06-modelos#a_função_lm|aqui]], na seção "Uma Palavra sobre o Argumento Fórmula". O comando para ajustar a equação da parábola é: mod2 <- update(mod1,.~. +I(init.h^2)) E agora podemos comparar os dois modelos: anova(mod1,mod2) O modelo de parábola é claramente superior, o que fica evidente se acrescentamos as linhas dos previstos pelos dois modelos ao gráfico: abline(mod1) cf.m2 <- coef(mod2) curve(cf.m2[1]+cf.m2[2]*x+cf.m2[3]*x^2, add=T, lty=2) Por fim, verifique o resumo do modelo selecionado: summary(mod2) Este procedimento pode ser usado sempre que se deseje avaliar se seus dados seguem uma equação de reta. Se você suspeita de curvatura, a comparação com o modelo quadrático é um teste simples e rápido.