Interface gráfica para ajuste de modelos de regressão não linear
Modelos de regressão não linear são considerados quando alguma informação a priori existe sobre o fenômeno. Essa informação pode ser, por exemplo, que a curva seja sempre crescente, típico para curvas de crescimento/acúmulo. Em geral, esses modelos têm parâmetros com interpretação física/química/biológica e alguns parâmetros não têm, mas estão presentes para conferir flexibilidade.
Talvez um dos grandes gargalos no ajuste de modelos não lineares seja atribuição de valores iniciais para o método numérico de estimação. Para parâmetros com interpretação é possível obter valores ao ver gráficos de dispersão das variáveis ou pelo menos se tem idéia do intervalo em que o parâmetro se encontra. Para parâmetros sem interpretação a obtenção de valores iniciais é realmente trabalhosa.
Um procedimento muito útil para obter valores iniciais é fazer um grid de valores iniciais e plotar a curva corresponde para cada conjunto sobre o diagrama de dispersão dos dados. Você usa como valores iniciais os correspondentes à curva que melhor se aproxima do comportamento dos dados. Por outro lado, criar um grid de valores, sobrepor cada curva ao gráfico, escolher à curva que melhor se aproxima aos dados, fornecer os valores iniciais e estimar os parâmetros é algo de demanda tempo se o número de curvas à ajustar for grande. Seria de grande ajuda se esse processo pudesse ser automatizado ou facilitado, por exemplo, por meio de uma interface gráfica com botões de ação e controle de parâmetros por meio do mouse.
O pacote rpanel
possui um conjunto de funções básicas que permitem construir interfaces gráficas para o usuário interagir com o R. O conjunto de funções compreende desde caixas de texto, caixas de seleção, deslizadores, até botões para disparo de alguma função. Usando essas ferramentas foi possível construir uma interface que permite selecionar o conjunto de dados, pré-ajustar a curva aos dados por meio de deslizadores e ajustar o modelo de regressão não linear. O código abaixo fornece as ferramentas para obter o que se vê nos gifs animados dessa matéria. Modificações devem ser feitas para outros dados/modelos. Minha intenção no futuro é aprimorar essas funções de forma que o usuário forneça o modelo e os limites dos intervalos para os parâmetros diretamente pela interface. Em caso de sucesso, isso pode virar um pacote ou complementar algum dos pacotes existentes dedicados à regressão não linear.
O primeiro exemplo usa gráficos do pacote graphics
, básico do R. Os dados são de 10 curvas de retenção de água no solo onde ajustou-se o modelo van Genuchten (1980). O segundo conjunto de dados é de um experimento para avaliar o impacto da desfolha na produtividade do algodão. Para esse caso adotou-se o modelo potência e gráficos do pacote lattice
. Clique na caixa abaixo para ver o código. Até à próxima ridícula.
#----------------------------------------------------------------------------- # definições da sessão require(rpanel) require(lattice) require(latticeExtra) require(grid) #----------------------------------------------------------------------------- # dados de curva de retenção de água no solo cra <- read.table("http://www.leg.ufpr.br/~walmes/data/cra_manejo.txt", header=TRUE, sep="\t") str(cra) cra$tens[cra$tens==0] <- 0.1 #----------------------------------------------------------------------------- # ver cras <- subset(cra, condi=="LVA3,5") cras <- with(cras, aggregate(cbind(umid), list(posi=posi, tens=tens, prof=prof), mean)) cras$caso <- with(cras, interaction(posi, prof)) xyplot(umid~log10(tens)|posi, groups=prof, data=cras, type=c("p","a")) #----------------------------------------------------------------------------- # ajuste do modelo van Genuchten de modo iterativo com rpanel # nlsajust é excutada quando clicar no botão "Ajustar" nlsajust <- function(panel){ ## ajuste do modelo não linear n0 <- try(nls(umid~tr+(ts-tr)/(1+(a*tens)^n)^(1-1/n), # modelo data=da, # dados start=start)) # valores iniciais ## em caso de não convergência imprime mensagem de erro if(class(n0)=="try-error"){ par(usr=c(0, 1, 0, 1)) text(0.5, 0.5, "Não convergiu!\nAproxime mais.", col="red", cex=2) } else { ## coloca à curva ajusta sobre os pontos with(as.list(coef(n0)), curve(tr+(ts-tr)/(1+(a*10^x)^n)^(1-1/n), add=TRUE, col=2)) ## salva o ajuste numa lista aju[[i]] <<- n0 } panel } vg <- function(panel){ ## seleciona o subconjunto dos dados i <<- panel$caso da <<- subset(cras, caso==i) ## lista com valores iniciais vindos dos deslizadores start <<- panel[c("tr","ts","a","n")] ## diagrama de dispersão plot(umid~log10(tens), data=da, ylab=expression(Conteúdo~de~água~(g~g^{-1})), xlab=expression(log[10]~Tensão~matricial~(kPa)), main=expression(f(x)==tr+(ts-tr)/(1+(a*x)^n)^{1-1/n})) ## sobrepõe a curva controlada pelos deslizadores with(start, curve(tr+(ts-tr)/(1+(a*10^x)^n)^(1-1/n), add=TRUE, col=2, lty=2)) panel } par(mar=c(4.1,4.2,3.1,1)) # cria objetos vazios que serão preenchidos durante processo da <- c(); start <- list(); aju <- list(); i <- c() # abre interface gráfica panel <- rp.control() # controla parâmetros do modelo rp.slider(panel, ts, 0.4, 0.8, initval=0.6, showvalue=TRUE, action=vg) rp.slider(panel, tr, 0, 0.5, initval=0.3, showvalue=TRUE, action=vg) rp.slider(panel, a, 0.1, 5, initval=1.3, showvalue=TRUE, action=vg) rp.slider(panel, n, 1, 5, initval=1.6, showvalue=TRUE, action=vg) # seleciona o conjunto de dados rp.listbox(panel, caso, vals=levels(cras$caso), title="Subconjunto", action=vg) # cria botão "Ajustar" rp.button(panel, action=nlsajust, title="Ajustar") sapply(aju, coef) # estimativas dos coeficientes lapply(aju, summary) # summary dos modelos ajustados #----------------------------------------------------------------------------- # dados de peso de capulhos em função da desfolha do algodão cap <- read.table("http://www.leg.ufpr.br/~walmes/data/algodão.txt", header=TRUE, sep="\t", encoding="latin1") str(cap) cap$desf <- cap$desf/100 cap <- subset(cap, select=c(estag, desf, pcapu)) cap$estag <- factor(cap$estag, labels=c("vegetativo","botão floral", "florescimento","maçã","capulho")) str(cap) xyplot(pcapu~desf|estag, data=cap, layout=c(5,1), xlab="Nível de desfolha artificial", ylab="Peso de capulhos") #----------------------------------------------------------------------------- # o exemplo a seguir mostra como usar com gráficos do pacote lattice # função adapatada para outro modelo de regressão não linear nlsajust <- function(panel){ n0 <- try(nls(pcapu~f0-f1*desf^exp(C), data=da, start=start)) if(class(n0)=="try-error"){ trellis.focus("panel", nivel, 1, highlight=FALSE) grid.text(x=0.5, y=0.5, label="Não convergiu!\nAproxime mais.", gp=gpar(col="red")) trellis.unfocus() } trellis.focus("panel", nivel, 1, highlight=FALSE) with(as.list(coef(n0)), panel.curve(f0-f1*x^exp(C), add=TRUE, col=2)) trellis.unfocus() print(coef(n0)) aju[[i]] <<- n0 panel } # função adapatada para outro modelo de regressão não linear ptn <- function(panel){ nivel <<- as.numeric(panel$nivel) i <<- levels(cap$estag)[nivel] da <<- subset(cap, estag==i) start <<- panel[c("f0","f1","C")] print({ xyplot(pcapu~desf|estag, data=cap, layout=c(5,1), col=1, xlab="Nível de desfolha artificial", ylab="Peso de capulhos", main=expression(f(x)==f[0]-f[1]*x^exp(C)), strip=strip.custom(bg="gray90")) }) trellis.focus("panel", nivel, 1, highlight=FALSE) with(start, panel.curve(f0-f1*x^exp(C), add=TRUE, col=2, lty=2)) trellis.unfocus() panel } da <- c(); start <- list(); aju <- list(); i <- c(); nivel <- c() panel <- rp.control() rp.slider(panel, f0, 10, 50, initval=30, showvalue=TRUE, action=ptn) rp.slider(panel, f1, 0, 30, initval=10, showvalue=TRUE, action=ptn) rp.slider(panel, C, -5, 5, initval=0, showvalue=TRUE, action=ptn) rp.listbox(panel, nivel, vals=1:nlevels(cap$estag), title="Painel", action=ptn) rp.button(panel, action=nlsajust, title="Ajustar") sapply(aju, coef) lapply(aju, summary) #-----------------------------------------------------------------------------
Fatorial duplo com dois tratamentos adicionais
Experimentos fatoriais são aqueles em que dois ou mais fatores estão presentes. Nos fatoriais incompletos nem todas as combinações possíveis estão presentes. Algumas das combinações podem não ser de interesse e por isso não são realizadas, alguma combinação é acidentalmente perdida durante o experimento ou certas combinações não identificam novos níveis, como por exemplo aplicação do produto A e B na dose 0 que, apesar dos rótulos diferentes, são o mesmo tratamento.
Um caso comum de fatorial incompleto são os fatoriais (completos) com a presença de tratamentos adicionais. Experimentos fatorais são realizados para verificar se os fatores interagem e então escolher combinação de níveis mais promissora para o processo estudado. Caso não exista interação, entende-se que os fatores têm efeito aditivos. Em qualquer caso, adicionam-se tratamentos adicionais ao experimento quando se deseja comparar o promissor aos de referência ou extremos.
Esse é um tipo de ensaio típico na avaliação do efeito de insumos agrícolas (fungicidas, insecticidas e herbicidas). Os ensaios testam produtos e suas doses e precisam ter as chamadas testemunhas para avaliar a eficiência das combinações geradas. Devido não ser um delineamento de níveis completamente cruzados a análise do experimento fatorial com tratamentos adicionais têm alguns detalhes que exigem mais atenção.
Nessa matéria faço análise de um fatorial com dois tratamentos adicionais. Os dados analisados foram retirados do Zimmermann (2004) de um experimento que estudou três fungicidas (benlate, captam, derosal) em três formas de aplicação (puro, antes de polímero, misturado com polímero). Adicionou-se dois tratamentos que foram a testemunha e aplicação do polímero puro. A resposta
avaliada foi o número de sementes infectadas com Fusarium em um total de 40 sementes.
O código está comentado para informar o leitor. São realizados a importação e reordenação dos níveis, a análise de variância com decomposição das somas de quadrados, contrastes entre efeitos, e testes sobre as médias ajustadas. Foi dada prioridade para os procedimentos matriciais de análise que podem ser adaptados para situações de desbalanceamento. Até a próxima ridícula.
#----------------------------------------------------------------------------- # definições da sessão require(lattice) require(latticeExtra) require(grid) require(MASS) require(doBy) require(gmodels) require(multcomp) #----------------------------------------------------------------------------- # leitura fu <- read.table( "http://www.leg.ufpr.br/~walmes/data/zimmermann_fusarium.txt", header=TRUE, sep="\t") str(fu) #----------------------------------------------------------------------------- # modificar ordem dos níveis, primeiro níveis fatoriais, depois adicionais levels(fu$fungicida) # ordem correta levels(fu$aplicacao) # precisa alterar levels(fu$aplicacao)[c(4,1,2,3,5)] # para ficar assim fu$aplicacao <- factor(fu$aplicacao, levels=levels(fu$aplicacao)[c(4,1,2,3,5)]) levels(fu$aplicacao) # agora ordem correta #----------------------------------------------------------------------------- # análise exploratória xtabs(~fungicida+aplicacao, fu) # os níveis não se cruzam completamente, portanto é um fatorial incompleto # no caso 3x3+2, fatorial com 2 tratamentos adicionais xyplot(infec40~tratamentos, fu) # resposta do tipo contagem de sucessos xyplot(infec40~tratamentos, fu, jitter.x=TRUE) #----------------------------------------------------------------------------- # Vamos primeito analisar fora da estrutura fatorial, ou seja, considerando # 3x3+2 = 11 níveis de um único fator. # Vamos considerar distribuição para reproduzir os resultados do livro # do Zimmermann (2004). Todavia, não é a análise mais correta devido os dados # serem discretos. A distribuição binomial parece ser mais indicada. m0 <- lm(infec40~tratamentos, data=fu) par(mfrow=c(2,2)); plot(m0); layout(1) # não tão ruim, tentar boxcox? m0 <- lm(infec40+1~tratamentos, data=fu) # soma 1 para eliminar os 0 boxcox(m0) # aponta log m1 <- lm(log(infec40+1)~tratamentos, data=fu) par(mfrow=c(2,2)); plot(m1); layout(1) # alguma melhoria? # Zimmermann usa asin(sqrt(x/40)). m2 <- lm(asin(sqrt(infec40/40))~tratamentos, data=fu) par(mfrow=c(2,2)); plot(m2); layout(1) # alguma melhoria? #----------------------------------------------------------------------------- # Visualmente todos os resíduos são semelhantes, apresentam fugas. lapply(list(m0, m1, m2), function(i) shapiro.test(residuals(i))) lapply(list(m0, m1, m2), function(i) ks.test(scale(residuals(i)), "pnorm")) # Pelos p-valores, a tranformação Box-Cox foi melhor, embora os resíduos ainda # não apresentem um comportamento satísfatório. # O foco não é testar transformações, então vamos seguir o que Zimmermann fez # e trabalhar com asin(sqrt(x/40)). #----------------------------------------------------------------------------- # análise exploratória com a variável transformada fu$y <- asin(sqrt(fu$infec40/40)) # cria a transformada xyplot(y~tratamentos, fu, jitter.x=TRUE) xyplot(y~fungicida, groups=aplicacao, fu, jitter.x=TRUE, type=c("p","a"), auto.key=TRUE) # Apresenta indicativos de efeito de interação e que os adicionais # são piores que o os fatoriais. #----------------------------------------------------------------------------- # análise na estrutura fatorial 3x3+2 m3 <- lm(y~fungicida*aplicacao, data=fu) anova(m3) # Os graus de liberdade são o que são pois: # fungicida 4: são 5 niveis (3 fatoriais, 2 adicionais), 5-1=4; # aplicacao 4: são 5, mas 2 estão confundidos (os adicionais), 5-2-1=2; # fun:apli 4: (3-1)*(3-1)=4; # O que queremos é a anova com somas de quadrados particionadas da seguinte # forma # fonte de variação gl # fungicidas 2 # aplicacao 2 # fun:apli 4 # fatorial vs adicionais 1 # polimero vs testemunha 1 #----------------------------------------------------------------------------- # Na minha opinião, obter esse quadro de análise de variância é desperdício # de tempo pois: # * não tem sentido prático comparar o efeito médio dos fatoriais contra o # efeito médio dos adicionais, principalmente sabendo que a interação # fungicida:aplicacao foi significativa; # * a interação fungicida:aplicacao significativa torna não interpretável # o teste para os efeitos principais; # * polímero vs testemunha tem só um grau de liberdade, posso comprar os # efeitos por um testes t; coef(m3) # tem NA para os efeitos não estimados (não presentes) efstm <- colnames(vcov(m3)) # nomes dos efeitos que foram estimados b <- coef(m3)[efstm] # só os estimados b # matriz de coeficientes para estimação das médias ajustadas X <- popMatrix(m3, effect=c("fungicida", "aplicacao")) str(X) Xn <- attr(X, "grid") # 25 linhas, 5x5=25 (faz pensando que é completo) Xo <- unique(fu[,c("fungicida","aplicacao")]) # que ocorrem (3x3+2) Xo #----------------------------------------------------------------------------- # calculando as médias ajustadas matricialmente rownames(X) <- apply(Xn, 1, paste, collapse=":") on <- apply(Xo, 1, paste, collapse=":") X <- X[on,efstm] # só linhas e colunas que correspondem à efeitos que ocorrem X%*%b # são as médias ajustadas with(fu, tapply(y, list(fungicida, aplicacao), mean)) # são médias amostrais # elas coincidem por ser um experimento regular (balanceado/ortogonal) #----------------------------------------------------------------------------- # fazendo contrastes # não funcionam com presenta de NA # contrast::contrast() # multcomp::glht() # para isso usar a gmodels::estimable() m4 <- aov(formula(m3), data=fu) # aov não NA em coef # polímero contra a testemunha p.vs.t <- matrix(X[10,]-X[11,], nrow=1) rownames(p.vs.t) <- "p.vs.t" estimable(m4, cm=p.vs.t) # fatorial contra adicional f.vs.a <- matrix(apply(X[1:9,], 2, mean)-apply(X[10:11,], 2, mean), 1) rownames(f.vs.a) <- "f.vs.a" estimable(m4, cm=f.vs.a) estimable(m4, cm=rbind(p.vs.t, f.vs.a)) # os dois de uma só vez # E aí, complicado fazer tudo isso, né? Qual a sensação em saber de que # pouco disso vai ser útil? Se deu interação teremos que estudar níveis # de um fator fixando os níveis do outro. #----------------------------------------------------------------------------- # fazendo a decomposição das somas de quadrados dentro da anova cf <- contrasts(C(fu$fungicida, helmert)) cf # o segredo está em montar os contrastes corretamente cf[,3] <- c(-1,-1,-1,1.5,1.5) # fatorial contra adicional cf[,4] <- c(0,0,0,-1,1) # polímero vs testemunha cf contrasts(fu$fungicida) <- cf m5 <- aov(y~fungicida*aplicacao, data=fu) summary(m5, expand=FALSE, split=list("fungicida"=list( "fungicida"=1:2, "fat.vs.adi"=3, "pol.vs.tes"=4))) # Esse procedimento depende de balanceamento pois as somas de quadrados são # sequênciais. #----------------------------------------------------------------------------- # Uma vez que deu interação é de prática fazer o desdobramento que consiste # em estudar o efeito de aplicação dentro dos níveis de fungicida e # vice-versa. Além do mais, por ter tratamentos adicionais, pode-se comparar # com os níveis dos fatoriais. Vamos listas o número de hipóteses que teremos # que fazer: # * 3 níveis de aplicação em cada fungicida: 3*choose(3,2)=9 # * 3 níveis de fungicida em cada aplicação: 3*choose(3,2)=9 # * testemunha contra cada um dos 9 níveis do 3x3: 9 # * polimero contra cada um dos 9 níveis do 3x3: 9 # Com isso são 4*9=36 comparações. Já que é assim, vou fazer todos # contra todos, choose(11, 2)=55, que o usuário pode olhar a que lhe # interessa, afinal de contas, de 36 para 55 não têm diferença mm <- model.matrix(m4) dim(mm) mm <- mm[,efstm] # modelo de classe lm equivalente à m5, sem NA m6 <- lm(y~-1+mm, data=fu) anova(m6) cbind(coef(m4), coef(m6)) # são o mesmo modelo #----------------------------------------------------------------------------- # médias ajustadas (ma) cima <- confint(glht(m6, linfct=X)) # intervalos de confiança para ma cima # tem cobertura *global* de 95% (e não indivídual) str(cima) cima$confint result <- cbind(Xo, cima$confint) str(result) segplot(fungicida:aplicacao~lwr+upr, data=result, xlab="Plantas infectadas (asin(sqrt(x/40))", ylab="Fungicida:aplicação", centers=Estimate, draw.bands=FALSE, segments.fun=panel.arrows, ends="both", angle=90, length=1, unit="mm") #----------------------------------------------------------------------------- # carrega algumas funções que serão úteis para comparar médias e fazer # gráficos com intervalos de confiança objs <- ls() objs <- c(objs, "apc", "prepanel.cbH", "panel.cbH") source("http://www.leg.ufpr.br/~walmes/ridiculas/ridiculas_functionsi.R", encoding="latin1") rm(list=ls()[!(ls()%in%objs)]) ls() #----------------------------------------------------------------------------- # comparando as médias ajustadas, choose(11, 2), contrastes de Tukey str(X) Xc <- apc(X) str(Xc) #----------------------------------------------------------------------------- # fazendo as comparações com a glht, método para correção de p-valor # fdr, o single-step demora muito c0 <- summary(glht(m6, linfct=Xc), test=adjusted(type="fdr")) c0 c0$focus <- "comparacoes" # para poder usar a cld() cld(c0) result$cld <- cld(c0)$mcletters$Letters #----------------------------------------------------------------------------- segplot(fungicida:aplicacao~lwr+upr, data=result, xlab="Plantas infectadas (asin(sqrt(x/40))", ylab="Fungicida:aplicação", centers=Estimate, draw.bands=FALSE, segments.fun=panel.arrows, ends="both", angle=90, length=1, unit="mm", panel=function(x, y, z, centers, subscripts, ...){ panel.segplot(x, y, z, centers=centers, subscripts=subscripts, ...) panel.text(centers, z, paste(format(centers, digits=2), result$cld[subscripts], sep=" "), pos=3) }) #----------------------------------------------------------------------------- # gráfico final result$desloc <- 0.05*c(0,-1,1,0,0)[match(result$aplicacao, levels(fu$aplicacao))] #png("f048.png", 500, 400) xyplot(y~fungicida, groups=aplicacao, data=fu, jitter.x=TRUE, amount=0.05, xlab="Fungicida", ylab=expression(asin(sqrt(x/40))), prepanel=prepanel.cbH, ly=result$lwr, uy=result$upr, auto.key=list(title="Aplicação", cex.title=1.2, columns=2, type="o", divide=1, lines=TRUE, points=FALSE))+ as.layer(xyplot(Estimate~fungicida, groups=aplicacao, data=result, type="l", ly=result$lwr, uy=result$upr, desloc=result$desloc, cty="bars", panel.groups=panel.cbH, panel=panel.superpose)) trellis.focus("panel", 1, 1, h=FALSE) x <- as.numeric(result$fungicida)+result$desloc y <- result$Estimate lab <- paste(format(y, digits=2), result$cld, sep=" ") a <- paste(rep("0", max(nchar(lab))), collapse="") grid.rect(x=unit(x, "native"), y=unit(y, "native"), width=unit(1, data=a, "strwidth"), height=unit(1.2, data=a, "strheight"), hjust=-0.1, gp=gpar(fill="gray90", col=NA, fontsize=10)) grid.text(lab, x=unit(x, "native"), y=unit(y, "native"), hjust=-0.2, gp=gpar(col="black", fontsize=10)) trellis.unfocus() #dev.off() #-----------------------------------------------------------------------------
Como fazer e interpretar o gráfico quantil-quantil
O gráfico quantil-quantil (q-q) é uma ferramenta muito útil para checar adequação de distribuição de frequência dos dados à uma distribuição de probabilidades. Situações como essa ocorrem principalmente na análise de resíduos de modelos de regressão onde o gráfico q-q é usado para verificar se os resíduos apresentam distribuição normal. O gráfico q-q é melhor que o histograma e o gráfico de distribuição acumulada empírica porque nós temos mais habilidade para verificar se uma reta se ajusta aos pontos do que se uma curva de densidade se ajusta a um histograma ou uma curva de probabilidade acumulada se ajusta à acumulada empírica. Compare às três visualizações com o código a seguir.
#----------------------------------------------------------------------------- # q-q vs histograma vs ecdf y <- rnorm(50) par(mfrow=c(1,3)) qqnorm(y); qqline(y) plot(density(y)) curve(dnorm(x, mean(y), sd(y)), add=TRUE, col=2) plot(ecdf(y)) curve(pnorm(x, mean(y), sd(y)), add=TRUE, col=2) # qual você sente mais segurança para verificar adequação? # nossos olhos têm mais habilidade para comparar alinhamentos #-----------------------------------------------------------------------------
Apesar de muito usado, poucos usuários conhecem o procedimento para fazê-lo e interpretá-lo. O procedimento é simples e pode ser estendido para outras distribuições de probabilidade, não apenas para a distribuição normal como muitos podem pensar. Além do mais, alguns padrões do gráfico q-q obdecem à certas características dados dados, como assimetria, curtose e discreticidade. Saber identificar essas características é fundamental para indicar uma transformação aos dados. Abaixo o código para o gráfico q-q para distribuição normal, q-q para qualquer distribuição e o q-q com envelope de confiança obtido por simulação. A execução e estudo do código esclarece o procedimento. O envelope de confiança por simulação torna-se proibitivo para grandes amostras pelo tempo gasto na simulação.
#----------------------------------------------------------------------------- # função para fazer o gráfico quantil-quantil da normal qqn <- function(x, ref.line=TRUE){ x <- na.omit(x) # remove NA xo <- sort(x) # ordena a amostra n <- length(x) # número de elementos i <- seq_along(x) # índices posicionais pteo <- (i-0.5)/n # probabilidades teóricas qteo <- qnorm(pteo) # quantis teóricos sob a normal padrão plot(xo~qteo) # quantis observados ~ quantis teóricos if(ref.line){ qrto <- quantile(x, c(1,3)/4) # 1º e 3º quartis observados qrtt <- qnorm(c(1,3)/4) # 1º e 3º quartis teóricos points(qrtt, qrto, pch=3) # quartis, passa uma reta de referência b <- diff(qrto)/diff(qrtt) # coeficiente de inclinação da reta a <- b*(0-qrtt[1])+qrto[1] # intercepto da reta abline(a=a, b=b) # reta de referência } } x <- rnorm(20) par(mfrow=c(1,2)) qqn(x) qqnorm(x); qqline(x) layout(1) #----------------------------------------------------------------------------- # função para fazer o gráfico quantil-quantil de qualquer distribuição qqq <- function(x, ref.line=TRUE, distr=qnorm, param=list(mean=0, sd=1)){ x <- na.omit(x) # remove NA xo <- sort(x) # ordena a amostra n <- length(x) # número de elementos i <- seq_along(x) # índices posicionais pteo <- (i-0.5)/n # probabilidades teóricas qteo <- do.call(distr, # quantis teóricos sob a distribuição c(list(p=pteo), param)) plot(xo~qteo) # quantis observados ~ quantis teóricos if(ref.line){ qrto <- quantile(x, c(1,3)/4) # 1º e 3º quartis observados qrtt <- do.call(distr, # 1º e 3º quartis teóricos c(list(p=c(1,3)/4), param)) points(qrtt, qrto, pch=3) # quartis, por eles passa uma reta de referência b <- diff(qrto)/diff(qrtt) # coeficiente de inclinação da reta a <- b*(0-qrtt[1])+qrto[1] # intercepto da reta abline(a=a, b=b) # reta de referência } } x <- rnorm(20) par(mfrow=c(1,2)) qqq(x) qqnorm(x); qqline(x) layout(1) x <- runif(20) qqq(x, ref.line=TRUE, distr=qunif, param=list(min=0, max=1)) x <- rgamma(20, shape=4, rate=1/2) qqq(x, ref.line=TRUE, distr=qgamma, param=list(shape=4, rate=1/2)) #----------------------------------------------------------------------------- # envelope para o gráfico de quantis (simulated bands) qqqsb <- function(x, ref.line=TRUE, distr=qnorm, param=list(mean=0, sd=1), sb=TRUE, nsim=500, alpha=0.95, ...){ x <- na.omit(x) # remove NA xo <- sort(x) # ordena a amostra n <- length(x) # número de elementos i <- seq_along(x) # índices posicionais pteo <- (i-0.5)/n # probabilidades teóricas qteo <- do.call(distr, # quantis teóricos sob a distribuição c(list(p=pteo), param)) plot(xo~qteo, ...) # quantis observados ~ quantis teóricos if(ref.line){ qrto <- quantile(x, c(1,3)/4) # 1º e 3º quartis observados qrtt <- do.call(distr, # 1º e 3º quartis teóricos c(list(p=c(1,3)/4), param)) points(qrtt, qrto, pch=3) # quartis, passa uma reta de referência b <- diff(qrto)/diff(qrtt) # coeficiente de inclinação da reta a <- b*(0-qrtt[1])+qrto[1] # intercepto da reta abline(a=a, b=b) # reta de referência } if(sb){ rdistr <- sub("q", "r", # função que gera números aleatórios deparse(substitute(distr))) aa <- replicate(nsim, # amostra da distribuição de referência sort(do.call(rdistr, c(list(n=n), param)))) lim <- apply(aa, 1, # limites das bandas 100*alpha% quantile, probs=c((1-alpha)/2,(alpha+1)/2)) matlines(qteo, t(lim), # coloca as bandas do envelope simulado lty=2, col=1) } } x <- rnorm(20) #png("f047.png", 400, 300) qqqsb(x, xlab="Quantis teóricos da distribuição normal padrão", ylab="Quantis observados da amostra", main="Gráfico quantil-quantil da normal") #dev.off() x <- rpois(50, lambda=20) qqqsb(x, distr=qpois, param=list(lambda=20)) #-----------------------------------------------------------------------------
A normalidade dos resíduos é um dos pressupostos da análise de regressão. São os resíduos e não os dados que devem apresentar normalidade. Se a distribuição dos dados, ou melhor, da sua variável resposta (Y) condicional ao efeito das suas variáveis explicativas for normal, os resíduos terão distribuição normal. Porém, se você aplica um teste de normalidade aos dados (Y) você não está considerando os efeitos das variáveis explicativas, ou seja, você está aplicando um teste na distribuição marginal de Y que não tem porque atender a normalidade. Todo teste de normalidade supõe que os dados têm uma média e uma variância e só os resíduos atendem essa premissa porque os dados (Y) têm média dependente do efeito das variáveis explicativas.
#----------------------------------------------------------------------------- # distribuição condicional vs distribuição marginal layout(1) x <- rep(1:4, e=50) y <- rnorm(x, mean=0+0.75*x, sd=0.1) da <- data.frame(x, y) plot(y~x, da) db <- split(da, f=da$x) #png("f046.png", 400, 300); par(mar=c(1,1,1,1)) plot(y~x, da, xlim=c(1,6), xaxt="n", yaxt="n", xlab="", ylab="") abline(a=0, b=0.75) lapply(db, function(d){ dnst <- density(d$y) lines(d$x[1]+dnst$y*0.1, dnst$x) abline(v=d$x[1], lty=2) }) points(rep(5, length(y)), da$y, col=2) abline(v=5, lty=2, col=2) dnst <- density(da$y) lines(5+dnst$y*2, dnst$x, col=2) text(3, 1, "distribuição\ncondicional") text(5.5, 1, "distribuição\nmarginal", col=2) #dev.off() par(mfrow=c(1,2)) qqnorm(da$y) qqnorm(residuals(lm(y~x, da))) layout(1) #-----------------------------------------------------------------------------
Muitos usuários preferem aplicar um teste de normalidade do que olhar para o gráfico q-q. Isso têm duas razões: (1) costume, o usuário sempre usou aplicativos para análise de dados que não dispõem de recursos gráficos, eles conduzem toda análise sem sequer ver os dados em um gráfico, (2) consideram subjetiva a análise gráfica. Do meu ponto de vista, a subjetividade está presente ao aplicar um teste também pois o usuário é quem escolhe o teste e o nível de significância. Outro ponto é que os testes de normalidade assumem independência entre as observações e os resíduos não são independentes e foram obtidos após a estimação de p parâmetros, o que não é considerado por nenhum teste. Ou seja, qualquer teste aplicado fornece um resultado aproximado. Mas o que de fato eu defendo é que a análise gráfica é indiscutivelmente mais informativa. Veja, se o teste rejeita a normalidade é porque os dados não apresentam distribuição normal por algum motivo. Quando você visualiza o q-q é possível explicar a fuga de normalidade que pode ser sistemática: (1) devido à desvio de assimetria, (2) de curtose, (3) à mistura de distribuições, (4) à presença de um outlier e (5) ao dados serem discretos. Essas são as principais causas de afastamento. Cada uma delas sugere uma alternativa para corrigir a fuga: tranformação, modelagem da variância, deleção de outlier, etc. Nesse sentido, como identificar esses padrões de fuga? É o que os gráficos animados vão mostrar. Até a próxima ridícula.
#----------------------------------------------------------------------------- # assimetria dir.create("frames") setwd("frames") png(file="assimet%03d.png", width=300, height=150) par(mar=c(1,1,1,1)) par(mfrow=c(1,2)) for(i in 10*sin(seq(0.01, pi-0.01, l=100))){ curve(dbeta(x, i, 10-i), 0, 1, xaxt="n", yaxt="n", xlab="", ylab="") y <- rbeta(100, i, 10-i) qqnorm(y, xaxt="n", yaxt="n", main=NULL, xlab="", ylab=""); qqline(y) } dev.off() # converte os pngs para um gif usando ImageMagick system("convert -delay 10 assimet*.png assimet.gif") # remove os arquivos png file.remove(list.files(pattern=".png")) #-----------------------------------------------------------------------------
#----------------------------------------------------------------------------- # mistura de distribuições png(file="mistura%03d.png", width=300, height=150) par(mar=c(1,1,1,1)) par(mfrow=c(1,2)) for(i in sin(seq(0, pi, l=100))){ curve(i*dnorm(x,0,1)+(1-i)*dnorm(x,0,6), -20, 20, xaxt="n", yaxt="n", xlab="", ylab="") y <- c(rnorm(ceiling(500*i),0,1), rnorm(floor(500*(1-i)),0,6)) qqnorm(y, xaxt="n", yaxt="n", main=NULL, xlab="", ylab=""); qqline(y) } dev.off() # converte os pngs para um gif usando ImageMagick system("convert -delay 10 mistura*.png mistura.gif") # remove os arquivos png file.remove(list.files(pattern=".png")) #-----------------------------------------------------------------------------
#----------------------------------------------------------------------------- # discreticidade png(file="discret%03d.png", width=300, height=150) par(mar=c(1,1,1,1)) par(mfrow=c(1,2)) for(i in c(1:100, 99:1)){ x <- 0:i px <- dbinom(x, i, 0.5) plot(px~x, type="h", xaxt="n", yaxt="n", xlab="", ylab="") y <- rbinom(100, i, 0.5) qqnorm(y, xaxt="n", yaxt="n", main=NULL, xlab="", ylab=""); qqline(y) } dev.off() # converte os pngs para um gif usando ImageMagick system("convert -delay 10 discret*.png discret.gif") # remove os arquivos png file.remove(list.files(pattern=".png")) #-----------------------------------------------------------------------------
Gerando código LATEX a partir do R para tabelas e matrizes
Oi Pessoal, ai vai mais uma dica!
Sabe aquela tabela (de dados ou resultados) ou matriz imensa que sai do R e que vc precisa montar ela no latex para gerar um trabalho, artigo, poster, etc? Afff q trabalhão, né?!
Mas não esquenta tem um comandinho q faz milagre no R, mas você precisa instalar o pacote xtable e pronto!
Olha um exemplinho ai para você se deliciar!!! É show!
###—###—###—###—###
x=rnorm(10,2,2)
y=rnorm(10,4,2)
A=data.frame(x,y)
require(xtable)
xtable(A)
## outro modelo
Modelo=summary(lm(x~y))
xtable(Modelo)
###—###—###—###—###
divirtam-se!
Si
Interpretação da matriz de covariância das estimativas
Quero falar de um erro comum de interpretação de resultados em análise de regressão através de um exemplo. Considere que você tem pessoas com diferentes pesos e diferentes alturas. Facilmente você aceita que uma pessoa mais alta tem maior peso, certo? Ou seja, existe uma correlação positiva entre peso e altura. Pois bem, vamos simular observações desse experimento e em seguida ajustar uma regressão linear simples para relação peso-altura.
require(MASS) set.seed(12345) da <- as.data.frame(mvrnorm(80, c(altura=30, peso=60), matrix(c(2,1.9,1.9,3), 2,2))) plot(peso~altura, da) m0 <- lm(peso~altura, data=da) # ajusta a reta summary(m0) # estimativa dos parâmetros abline(m0) # adiciona uma reta ao gráfico vcov(m0) # covariância das estimativas
Aqui está o ponto que eu quero comentar: a interpretação da matriz de covariância. Perceba que a covariância foi negativa. Já vi gente interpretando isso da seguinte forma: espera-se que pessoas mais altas tenham menor peso. Duplamente errado. Primeiro nos sabemos por experiência que a correlação de peso e altura é positiva. Segundo, a matriz de covariância se refere as estimativas dos parâmetros e não as variáveis envolvidas. Nunca a interprete dessa forma.
Então como interpretar? Bem, a matriz de covariância das estimativas, é um reflexo da função objetivo ao redor da solução. A função objetivo nesse caso é minimizar a soma de quadrados (SQ). Então se eu aumento o valor do parâmetro b0, o parâmetro b1 tem que diminuir para compensar, ou seja, para diminuir a SQ. Ocorre um efeito compensatório aqui. Além do mais, esse efeito pode ser eliminado facilmente aplicando uma reparametrização, como por exemplo, centrar a covariável na média.
da$alturac <- with(da, altura-mean(altura)) # centrando a covariável par(mfrow=c(1,2)) plot(peso~altura, da) abline(m0) plot(peso~alturac, da) m1 <- lm(peso~alturac, data=da) # ajusta a reta summary(m1) # estimativa dos parâmetros abline(m1) # adiciona uma reta ao gráfico vcov(m1) # covariância das estimativas require(ellipse) par(mfrow=c(1,2)) # regiões de confiança conjunta plot(ellipse(vcov(m0)), type="l") plot(ellipse(vcov(m1)), type="l") cov2cor(vcov(m0)) cov2cor(vcov(m1))
Perceba que os resultados são os mesmos em termos de estatísticas de teste e medidas de ajuste. Isso é esperado por eu só fiz uma traslação da altura. Mas o importante é que agora a matriz de covariância tem covariâncias praticamente nulas, que é resultado da translação. O intercepto é então o valor esperado para alturac igual a zero (centro dos dados). Veja porque essa covariância é zero: se eu aumentar b0, não adiante eu alterar b1 que não vai diminuir a SQ, e vice-versa, porque agora eles são ortogonais. Ortogonalidade entre parâmetros é uma característica desejável pois permite você inferir sobre um deles sem considerar o outro. Além disso, tem vantagens do ponto de vista de estimação por métodos numéricos de otimização. Em outras palavras, pegando conceitos de probabilidade, se a distribuição amostral de duas variáveis aleatórias é normal bivariada com covariância nula, a distribuição condicional de A|B é igual a distribuição marginal de A pela independência. Reforçando, eu não preciso conhecer valores de B para descrever A. Considerando tudo que foi argumentado, é sempre preferível que você adote o modelo que apresente menor covariância entre os parâmetros. Até a próxima ridícula.