Archive

Posts Tagged ‘xyplot’

Equação junto ao gráfico na lattice

Reposta à salinidade para quatro espécies de craca. A regulação osmótica está ligada à capacidade adaptativa do organismo ao ambiente. O gráfico mostra os valores observados e a curva de regressão estimada acompanhada de bandas de confiança (95%). As equações em cada gráfico representam o modelo ajustado que descreve a curva.

Os polinômios são os modelos de regressão mais empregados em algumas áreas do conhecimento. Isso não quer dizer que eles sejam modelos bons. Na verdade, isso depende do cada um define como um bom modelo. Os polinômios são uma aproximação local. Isso significa que eles não devem ser usados para predições fora da amplitude observada das covariáveis (região experimental). Em outras palavras, não se deve fazer extrapolações, mas pode-se, com certos cuidados, fazer interpolações. Os polinômios são desprovidos de interpretação a partir do segundo grau. A partir do segundo grau, apenas o intercepto tem interpretação pois é o valor da função quando a covariável (x) é zero. Por isso, interpretações em geral são perdas de tempo. O que pode-se fazer são julgamentos quando ao sinal do coeficiente de maior grau que por exemplo, num modelo quadrático, indica se a curvatura é para cima (ou não). Outro ponto é que os polinômios apresentam problema de colinearidade. Para superar essa deficiência pode-se centrar a covariável ou usar polinômios ortogonais. Apesar do polinômio já ser sem interpretação, o uso de polinômios ortogonais é menos interpretável ainda pois os coeficientes não são taxas em relação à x, x², x³…, mas para uma transformação linear de x, essa que induz a ortogonalidade. Resumindo, os polinômios ortogonais são úteis pois minimizam o problema de colineridade. Então, do ponto de vista de reportar uma equação de predição, aquela que descreve o valor esperado de y (resposta) em função de x, é mais útil usar o polinômio cru. Isso porque, mesmo sem ter interpretação, toda pessoa sabe traçar a curva de um polinômio, é algo familiar, visto desde o ensino médio, e novamente, não que isso seja um ponto positivo para um modelo. De qualquer forma, o post de hoje ensina a ajustar polinômios e representar no gráfico os valores observados, os valores ajustados com intervalo de confiança e a equação dentro do gráfico. Os procedimentos abaixo eu desenvolvi para atender demandas de consultoria. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# definições da sessão

require(lattice)
require(latticeExtra)
require(grid)

#-----------------------------------------------------------------------------
# leitura

osmo <- read.table("http://www.leg.ufpr.br/~walmes/ridiculas/osmo.txt",
                   header=TRUE, sep="\t")
str(osmo)

xyplot(Osm~Salinidade|Espécie, data=osmo, type=c("p","smooth"))

#-----------------------------------------------------------------------------
# vou desconsiderar o fato dos dados serem heterocedasticos e ajustar um lm
# em outra matéria trataremos os desvios de pressupostos
# ajuste usando polinômios ortogonais

m0 <- aov(Osm~Espécie*poly(Salinidade,2), data=osmo)
par(mfrow=c(2,2)); plot(m0); layout(1)

qqmath(~residuals(m0)|osmo$Espécie)
xyplot(sqrt(abs(residuals(m0)))~fitted(m0)|osmo$Espécie, type=c("p","smooth"))
# desvio de pressupostos, trataremos em outra matéria do blog

summary(m0) # inferências sem validade pois existe desvio de pressupostos

#-----------------------------------------------------------------------------
# partição das somas de quadrados, declara Salinidade dentro de Espécie

m1 <- aov(Osm~Espécie/poly(Salinidade,2), data=osmo)

names(coef(m1))
attr(terms(m1), "term.labels")
ncoef <- 2 # número de efeitos por expécie=2 (linear e quadrático)
ngrup <- 4 # número de níveis de espécie=4
split.list <- split(matrix(1:(ncoef*ngrup), ngrup, ncoef), f=1:ngrup)
names(split.list) <- levels(osmo$Espécie)
split.list <- lapply(split.list,
                     function(x){ names(x) <- c("lin","qua"); x})
split.list <- list(a=split.list)
names(split.list) <- attr(terms(m1), "term.labels")[2]

summary(m1, split=split.list)

split.list <- as.list(1:(ncoef*ngrup))
names(split.list) <- c(outer(levels(osmo$Espécie), c("lin","qua"),
                             paste, sep=": "))
split.list <- list(a=split.list)
names(split.list) <- attr(terms(m1), "term.labels")[2]

summary(m1, split=split.list)

#-----------------------------------------------------------------------------
# prepara as equações para colocar no gráfico, função para polinômiais

coefs <- coef(lm(Osm~-1+Espécie/(Salinidade+I(Salinidade^2)), data=osmo))
coefs <- as.numeric(formatC(coefs, format="f", digits=3, flag="-"))

get.expression <- function(coef, digits){
  coef <- as.numeric(formatC(coef, format="f", digits=digits, flag="-"))
  postext <- c("","*x", paste("*x^", 2:(length(coef)-1), sep=""))
  sig <- c(ifelse(coef[1]<0, "-", ""), ifelse(coef[-1]<0, "-", "+"))
  txt <- paste(paste(sig, abs(coef), postext, sep=""), collapse=" ")
  eval(parse(text=paste("expression(y==", txt, ")", sep="")))
}

coefs <- split(matrix(coefs, ncol=3), f=1:4)
eqn <- lapply(coefs, get.expression, digits=3)
eqn

#-----------------------------------------------------------------------------
# nomes científicos em itálico nas tarjas

fl <- as.expression(lapply(levels(osmo$Espécie),
    function(x){ bquote(italic(.(x))) }))
fl

#-----------------------------------------------------------------------------
# predição

pred <- expand.grid(Espécie=levels(osmo$Espécie), Salinidade=seq(0,40,by=2))
pred <- cbind(pred, predict(m1, newdata=pred, interval="confidence"))

#-----------------------------------------------------------------------------
# carrega funções que serão necessárias para fazer as bandas de confiança

source("http://www.leg.ufpr.br/~walmes/ridiculas/ridiculas_functions.R")
grep("ciH", ls(), value=TRUE) # nome das funções (painéis) que vamos usar

#-----------------------------------------------------------------------------
# gráfico final

#png("f035.png", width=600, height=500)
xyplot(Osm~Salinidade|Espécie, data=osmo, col=1,
       xlab=expression(Salinidade~(ppm)),
       ylab=expression(Osmolalidade~(mOsm~kg^{-1}~de~H[2]*O)),
       strip=strip.custom(bg="gray90", factor.levels=fl))+
  as.layer(xyplot(fit~Salinidade|Espécie, data=pred, type="l", col=1,
                  ly=pred$lwr, uy=pred$upr,
                  prepanel=prepanel.ciH, panel=panel.ciH))
dimLay <- trellis.currentLayout()
for(col in 1:ncol(dimLay)){
  for(row in 1:nrow(dimLay)){
    trellis.focus("panel", col, row, highlight=FALSE)
    grid.text(label=eqn[[dimLay[row,col]]],
              x=0.03, y=1-0.03, just=c("left","top"))
  }
}
trellis.unfocus()
#dev.off()

#-----------------------------------------------------------------------------
Categorias:gráficos Tags:, ,

Uma alternativa para apresentar bandas de confiança em regressão

Valores preditos para índice agronômico de milho em função da dose de nitrogênio e cultivar. Regiões em cinza são os intervalos de confiança (95%) para os valores preditos.

Neste post novamente vou apresentar a análise dos dados de índice agronômico de milho que iniciei em outros posts (Predição em modelos de regressão com blocos, Como fazer regressão polinomial com diferentes graus). A única novidade será a representação gráfica, que ao invés de separar as curvas por nível de cultivar, vou apresentar no mesmo gráfico usando uma região escura para representar os intervalos de confiança para os valores preditos. Esse gráfico facilita perceber se os intervalos têm sobreposição ou não. No entanto, ter (ou não) sobreposição não aponta aceitação (ou rejetição) da hipótese de igualdade dos valores preditos. Hipóteses devem ser definidas antes da realização do experimento e avaliadas por testes adequados. Foi necessário desenvolver uma função complemento para os gráficos da lattice para produzir as regiões de confiança. No script abaixo eu reproduzo as análises e faço o gráfico que mencionei. Até a próxima ridícula!

#-----------------------------------------------------------------------------
# dados de índice agronômico de milho

da <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/anovareg.txt",
                 header=TRUE, sep="\t")
str(da)

#-----------------------------------------------------------------------------
# ajuste do modelo. veja os posts anteriores a esse mencionados no texto.

levels(da$cultivar)
da$Pi.C <- with(da, ifelse(da$cultivar=="Pioneer-B815", dose^3, 0))
m3 <- lm(indice~bloco+cultivar*(dose+I(dose^2))+Pi.C, data=da,
         contrast=list(bloco=contr.sum))
summary(m3)

#-----------------------------------------------------------------------------
# faz o gráfico final de predição pelo método matricial

# data.frame da predição
pred <- expand.grid(cultivar=levels(da$cultivar), dose=seq(0,300,10))
pred$Pi.C <- with(pred, ifelse(cultivar=="Pioneer-B815", dose^3, 0))

# matriz de predição sem colunas para os efeitos de blocos
X <- model.matrix(~cultivar*(dose+I(dose^2))+Pi.C, data=pred)

# vetor de estimativas sem os efeitos de blocos
b <- coef(m3)[-grep("bloco", names(coef(m3)))]

# margem de erro, produto do erro padrão e quantil
me <- predict(m3, newdata=cbind(bloco=levels(da$bloco)[1], pred),
              se.fit=TRUE)$se*qt(0.975, df=df.residual(m3))

# data.frame com valores preditos e IC
pred$fit <- X%*%b
pred$lwr <- pred$fit-me; pred$upr <- pred$fit+me
str(pred)

#-----------------------------------------------------------------------------
# paineis feitos para lattice para apresentar os intervalos de confiança como
# regiões de outra cor, úteis para sobrepor ICs

prepanel.ciH <- function(x, y, ly, uy, subscripts, ...){
  x <- as.numeric(x)
  ly <- as.numeric(ly[subscripts])
  uy <- as.numeric(uy[subscripts])
  list(ylim=range(uy, ly, finite=TRUE), xlim=range(x))
}

panel.ciH <- function(x, y, ly, uy, subscripts, ...){
  y <- as.numeric(y)
  x <- as.numeric(x)
  or <- order(x)
  ly <- as.numeric(ly[subscripts])
  uy <- as.numeric(uy[subscripts])
  panel.polygon(c(x[or], x[rev(or)]),
                c(ly[or], uy[rev(or)]),
                col=1, alpha=0.03, border=NA)
  panel.lines(x[or], ly[or], lty=3, lwd=0.5, col=1)
  panel.lines(x[or], uy[or], lty=3, lwd=0.5, col=1)
  panel.xyplot(x, y, subscripts=subscripts, ...)
}

#-----------------------------------------------------------------------------
# gráfico com observados, preditos e IC pelo método matricial

require(lattice)
require(latticeExtra) # as.layer()
require(grid)

#png("f020.png", w=500, h=350)
with(pred,
     xyplot(fit~dose, groups=cultivar, type="l", lwd=2,
            strip=strip.custom(bg="gray90"),
            xlab="Dose de nitrogênio (kg/ha)", ylab="Índice agronômico",
            sub=list("predição livre do efeito de bloco", cex=0.8, font=3),
            ly=lwr, uy=upr,
            prepanel=prepanel.ciH,
            panel=panel.superpose,
            panel.groups=panel.ciH,
            key=list(space="top", columns=1,
              points=list(pch=1), text=list("valores observados"),
              lines=list(col=1), text=list("valores preditos"),
              rectangle=list(col="gray90", border=NA), text=list("IC 95%"))
            ))+
  as.layer(xyplot(indice~dose, groups=cultivar, data=da, jitter.x=TRUE))
draw.key(list(columns=1,
              lines=list(lty=1, lwd=2,
                col=trellis.par.get("superpose.symbol")$col[1:3]),
              text=list(levels(da$cultivar))),
         draw=TRUE, vp=viewport(x=unit(0.75, "npc"), y=unit(0.25, "npc")))
#dev.off()

#-----------------------------------------------------------------------------

Como fazer regressão polinomial com diferentes graus

Ajuste de modelos de regressão polinomial com diferentes graus para o polinômio. Para as cultivares da esquerda o polinômio é quadrático, para a da direita o polinômio é cúbico.

No post Predição em modelos de regressão com blocos foram apresentados procedimentos para estimação e predição de parâmetros em modelos de regressão com blocos. Nesse post aqui eu vou estender o exemplo permitindo o ajuste com diferentes graus para o polinômio que descreve a relação da resposta com a preditora contínua. Nos códigos abaixo vou apresentar aspectos dos polinômios ortogonais, do cálculo do R² e de como fazer o ajuste com diferentes graus de polinômio em um mesmo modelo, ou em outras palavras, em uma mesma lm(). No final, como de costume, confecciono o gráfico resultado do modelo ajustado com as bandas de confiança. Mais uma vez deixei o código comentado passo a passo para orientar o leitor. Até a próxima ridícula!

#-----------------------------------------------------------------------------
# dados de índice agronômico de milho de um experimento

da <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/anovareg.txt",
                 header=TRUE, sep="\t")
str(da)

#-----------------------------------------------------------------------------
# objetivo: apresentar procedimentos para ajuste de polinômios e como ajustar
# polinômios de diferentes graus em uma mesma análise

#-----------------------------------------------------------------------------
# vamos apresentar essas funções para o caso do ajuste de uma função e fazer
# de conta que não temos blocos apenas por facilidade de exposição

levels(da$cultivar)
br <- subset(da, cultivar=="BR-300")
str(br)

#-----------------------------------------------------------------------------
# ajustar diminuindo o nível do polinômio até que não tenhamos falta de ajuste
# fazer isso de diferentes formas:
# 1) fórmula explicita e anova seqüencial
# 2) polinômios ortogonais e teste t

# 1) fórmula explicita e anova seqüencial
m1.5 <- lm(indice~dose+I(dose^2)+I(dose^3)+I(dose^4)+I(dose^5), data=br)
anova(m1.5) # testes F apontam que os termos > 2 grau não contribuem
m1.2 <- lm(indice~dose+I(dose^2), data=br)
anova(m1.2, m1.5) # o abandono desses termos não resulta em falta de ajuste
anova(m1.2)

# 2) polinômios ortogonais e teste t
m2.5 <- lm(indice~poly(dose, degree=5), data=br)
anova(m2.5)   # anova não separa os termos
summary(m2.5) # quando usamos polinômios orotoganis, o valor t^2 = F
cbind("t^2"=summary(m2.5)$coeff[-1,3]^2, "F"=anova(m1.5)[-6,4]) # viu?

m2.2 <- lm(indice~poly(dose, degree=2), data=br)
summary(m2.2)

#-----------------------------------------------------------------------------
# antigamente experimentos agronômicos com fatores quantitativos eram
# planejados para polinômios ortogonais. assim usava espaçamento regular entre
# níveis. muitas pessoas acham que isso é regra mas era apenas algo imposto
# pelas limitações de computação da época. com polinômios ortogonais a matriz
# X'X é ortogonal, então facilmente invertível e obtém-se as estimativas dos
# parâmetros. com essas estimativas e erros padrões aplicava-se o teste t (ao
# invés do teste F) para testar a significância dos termos. os temos não
# significativos eram abandonados mas a estimação não precisava ser refeita
# por causa da ortogonalidade. observe os procedimentos abaixo.

X <- model.matrix(m2.5) # matriz de delineamento ortogonal
round(colSums(X), 4) # soma nas colunas dá zero
round(t(X)%*%X, 4)   # ortogonal, portanto a inversa é a inversa da diagonal
t(X)%*%br$indice     # as estimativas dos parâmetros, intercepto dividir por n

cbind(X[,2:3], model.matrix(m2.2)[,2:3])[1:6,] # iguais colunas
coef(m2.5); coef(m2.2) # as estimativas são iguais, não requer outro ajuste

#-----------------------------------------------------------------------------
# no entanto, esses coeficientes não tem interpretação pois a matriz não está
# na escala das doses. aliás, polinômio não tem interpretação em nenhuma
# escala. eles representam apenas uma aproximação local. interpretar a
# magnitude das estimativas é perda de tempo. compare apenas o rank e a
# direção (sinal). faça o teste de hipótese e interprete os valores preditos.
# caso queira apresentar a equação ajustada, terá que obter as estimativas
# sem polinômios ortogonais.

#-----------------------------------------------------------------------------
# obter as estimativas sem polinômios ortogonais com a função poly()

m2.2 <- lm(indice~poly(dose, degree=2, raw=TRUE), data=br)
summary(m2.2) # estimativas na escala do fator, são as mesmas do modelo m1.2

#-----------------------------------------------------------------------------
# muitos tutoriais/aplicativos obtém o R² do ajuste fazendo ajuste às médias.
# *eu* não considero essa prática porque esse R² não reflete a razão
# ruído/sinal dos dados para com o modelo, mas da média dos dados para com o
# modelo. esse tipo de procedimento só é possível quando temos repetição de
# níveis. veja os três exemplos simulados a seguir

desv <- rnorm(12,0,1)      # desvios aleatórios
a <- 0; b <- 1             # valores dos parâmetros
x1 <- seq(1,12,1)          # sem repetição de nível
x2 <- rep(seq(3,12,3),e=3) # com repetição de nível
y1 <- a+x1*b+desv          # valores observados de y1
y2 <- a+x2*b+desv          # valores observados de y2

m1 <- lm(y1~x1); summary(m1) # R² correto
m2 <- lm(y2~x2); summary(m2) # R² correto

# R² calculado com relação ao desvio dos ajustados para as médias
md <- tapply(y2, x2, mean)            # tiramos a média das repetições
aj <- unique(round(fitted(m2),5))     # valores ajustados pelo modelo
r2 <- 1-sum((md-aj)^2)/sum((md-mean(md))^2); r2 # R² como desvio das médias

#-----------------------------------------------------------------------------
# não sei dizer qual a justificativa para o uso desse procedimento para obter
# o R². é um procedimento atraente porque o R² será *sempre* maior e talvez
# essa seja a razão da sua adoção (que na *minha opinião* é errada). esse R²
# não associa desvio das observações. então o valor desse R² não é a
# porcentagem de explicação que o modelo dá para as observações, e sim para a
# média delas! como assim se estou modelando as observações e não a média
# amostral delas? veja esse exemplo mais abrupto

#-----------------------------------------------------------------------------
# criando um novo par x, y

x3 <- x2+c(-0.001,0,0.001) # deslocamento mínimo faz perder as repetições
y3 <- a+x3*b+desv

cbind(y2,y3,x2,x3)

# ambos pares x, y são praticamente a mesma coisa e o R² será muito próximo
# para ambos. porém, para x2 pode-se usar o R² relativo a desvio da média.
# se eu calcular esse R² ele vai diferir muito pois representa outra razão
# ruído/sinal. *eu* não gosto disso.

m3 <- lm(y3~x3)
c(summary(m1)$r.sq, r2, summary(m2)$r.sq, summary(m3)$r.sq) # compara R²

#-----------------------------------------------------------------------------
# para mais críticas sobre R² procure sobre a análise do quarteto de Anscombe.

#-----------------------------------------------------------------------------
# ajustar modelos de regressão para os 3 níveis de cultivar e decompor a SQ.
# aqui usaremos o fator de níveis ordenados que tem como padrão o contraste
# de polinômios ortogonais. na poly() conseguimos especificar o grau e no
# fator ordenado não conseguimos, mas as somas de quadrados podem ser
# facilmente desdobradas se usarmos aov().

m0 <- aov(indice~bloco+cultivar*ordered(dose), data=da)
summary(m0)    # pode-se usar anova(m0) também
summary.lm(m0) # quadro de estimativas dos parâmetros

#-----------------------------------------------------------------------------
# para obter o quadro de análise de variância com desdobramentos dos termos
# do polinômio dentro dos níveis de cultivar precisamos mudar a fórmula do
# modelo, é apenas uma mudança a nível de colunas mas o espaço do modelo ainda
# é o mesmo. compare as colunas das matrizes obtidas com os comando abaixo.

aux <- expand.grid(cultivar=gl(2,1), dose=1:3) # estrutura experimental
model.matrix(~cultivar*ordered(dose), aux)     # cultivar cruzado com dose
model.matrix(~cultivar/ordered(dose), aux)     # dose aninhado em cultivar

#-----------------------------------------------------------------------------
# ajustando modelo com fórmula corresponde ao desdobramento desejado
# o procedimento abaixo é para criar a lista que especifica a partição nas
# somas de quadrados, envolve algumas funções que vale a pena consultar o help

da$D <- ordered(da$dose) # cria uma nova coluna apenas por conforto
m1 <- aov(indice~bloco+cultivar/D, data=da) # modelo com a fórmula conveniente

m1$assign # valores iguais a 3 representam os estimativas do 3 termo
tm <- names(coef(m1))[m1$assign==3]                  # nomes dos termos
tm <- gsub("^cultivar(.*):.*(.{2})$", "\\1-\\2", tm) # edita os nomes
ord <- c(sapply(levels(da$cultivar), grep, x=tm, simplify=TRUE)) # ordena
split <- as.list(ord)   # cria a lista com a posição dos termos
names(split) <- tm[ord] # atribui nome aos elementos da lista

#-----------------------------------------------------------------------------
# quadro de análise de variância com somas que quadrados particionadas.
# com ela escolhemos o grau do polinômio adequado para cada cultivar.

summary(m1, split=list("cultivar:D"=split))

#-----------------------------------------------------------------------------
# uma das cultivares requer 3 grau, as outras ficam com o 2 grau. juntar os
# termos não significativos do polinômio para criar o termo falta de ajuste.

do.call(c, split) # ver os índices e fazer o agrupamento gerando nova split
split <- list(Ag.L=1, Ag.Q=4,         Ag.lof=c(7,10,13),
              BR.L=2, BR.Q=5,         BR.lof=c(8,11,14),
              Pi.L=3, Pi.Q=6, Pi.C=9, Pi.lof=c(12,15))
summary(m1, split=list("cultivar:D"=split)) # lof não significativos, ok!

#-----------------------------------------------------------------------------
# não iremos apresentar as estimativas do modelo m1 nem seus valores ajustados
# porque esses se referem ao modelo completo. temos que ajustar o modelo
# reduzido removendo os temos apontado pelo teste F. este modelo é o modelo
# que apresentaremos as estimativas e faremos predição. note que o grau do
# polinômio agora muda com o nível de cultivar. teremos que arrumar um
# artifício ajustar um modelo considere isso. então vou criar uma indicadora
# para multiplicar o termo cúbico e sumir com ele quando não for necessário.

levels(da$cultivar)
da$ind <- 0; da[da$cultivar=="Pioneer-B815",]$ind <- 1 # a indicadora do nível
da$ind # o valor 1 é associado ao nível de cultivar que tem o termo cúbico

#-----------------------------------------------------------------------------
# o termo cúbico só vai existir para o nível Pionner por causa da indicadora.
# vai aparecer um NA para os outros níveis

m2 <- lm(indice~bloco+cultivar*(dose+I(dose^2)+I(ind*dose^3)), data=da)
summary(m2)

#-----------------------------------------------------------------------------
# os NA implicam em algumas conseqüência como no funcionamento da predict().
# para evitá-los podemos criar uma coluna para representar esse termo cúbico e
# e declara uma fórmula de modelo que não retorne NA.

da$Pi.C <- with(da, ind*dose^3) # é a coluna do cubo da dose para Pioneer
m3 <- lm(indice~bloco+cultivar*(dose+I(dose^2))+Pi.C, data=da,
         contrast=list(bloco=contr.sum))
summary(m3)

#-----------------------------------------------------------------------------
# faz o gráfico final de predição pelo método matricial

# data.frame da predição
pred <- expand.grid(cultivar=levels(da$cultivar), dose=seq(0,300,10))
pred$Pi.C <- with(pred, ifelse(cultivar=="Pioneer-B815", dose^3, 0))

# matriz de predição sem colunas para os efeitos de blocos
X <- model.matrix(~cultivar*(dose+I(dose^2))+Pi.C, data=pred)
head(X)

# vetor de estimativas sem os efeitos de blocos
b <- coef(m3)[-grep("bloco", names(coef(m3)))]

# vetor de erros padrões das estimativas
se <- predict(m3, newdata=cbind(bloco=levels(da$bloco)[1], pred),
              se.fit=TRUE)$se

# quantil superior da distribuição t para um IC de 95%
tc <- qt(0.975, df=df.residual(m3))

# margem de erro
me <- se*tc

# data.frame com valores preditos e IC
pred$fit <- X%*%b
pred$lwr <- pred$fit-me; pred$upr <- pred$fit+me

require(reshape)

pred <- melt(pred[,-3], id=c("cultivar","dose"))
names(pred)[ncol(pred)] <- "indice"
da$variable <- "obs"
pred <- rbind(da[,names(pred)], pred)

#-----------------------------------------------------------------------------
# gráfico com observados, preditos e IC pelo método matricial

require(lattice)

#png("f019.png", w=500, h=250)
xyplot(indice~dose|cultivar, groups=variable, data=pred,
       distribute.type=TRUE, layout=c(3,1),
       type=c("l","l","p","l"), col=c(2,3,1,3), lty=c(1,2,0,2),
       strip=strip.custom(bg="gray90"),
       xlab="Dose de nitrogênio (kg/ha)", ylab="Índice agronômico",
       sub=list("predição livre do efeito de bloco", cex=0.8, font=3),
       key=list(space="top", columns=3, type="o", divide=1,
         lines=list(pch=c(1,NA,NA), lty=c(0,1,2), col=c(1,2,3)),
         text=list(c("valores observados","valores preditos","IC 95%"))))
#dev.off()

#-----------------------------------------------------------------------------

Predição em modelos de regressão com blocos

Padrão de reposta de cultivares de milho à adubação com nitrogênio.

Em experimentos onde não se dispõe de unidades experimentais uniformes (UI), usualmente se agrupa UI semelhantes em grupos que usualmente recebe o nome de bloco. Em cada bloco as unidades experimentais são uniformes, e entre classes são diferentes. Uniformidade aqui significa que entre essas UI não é possível atribuir uma causa/origem as pequenas diferenças existentes. São exemplos de blocos em experimentos agrícolas: nível topográfico, idade inicial, posição na casa de vegetação, operador do sistema, local de origem, horário de avaliação, dentre muitos outros que podem ser citados.

O uso de um fator de blocagem permite que parâmetros e contrastes sejam estimados com maior precisão porque parte da diferença entre UI foi identificada pelo arranjo em blocos. No entanto, o uso de blocos sem a necessidade pode elevar a taxa de erro tipo I e o uso de blocos de forma incorreta pode reduzir o poder dos testes de hipótese. Ao iniciar uma pesquisa, consulte um especialista em planejamento de experimentos.

Os blocos representam uma causa de variação presente mas que não é o foco da pesquisa. O pesquisador deve declarar o efeito de bloco no modelo estatístico da análise mas as inferências só serão aplicadas aos demais fatores do estudo. Em algum ponto, além de estimar e comparar parâmetros, há o interesse de fazer a predição da resposta. Nessa predição, normalmente, deseja-se que estejam presentes apenas os termos de interesse da pesquisa e remova-se o efeito dos blocos.

Neste post vou apresentar procedimentos para obter valores preditos livre dos efeitos dos blocos. Farei uso da função predict() e mostrarei como fazer por operação de matrizes por meio de funções lineares. Deve-se fazer algumas adaptações nos procedimentos quando o experimento possui observações não disponíveis (parcelas perdidas). Algumas adaptações serão necessárias também para uso modelos lineares generalizados. Nesses procedimentos o efeito de bloco é fixo e não aleatório.

Os dados são de experimento conduzido com parcelas dispostas em blocos. Os fatores de interesse são a cultivar de milho e a dose de nitrogênio. O objetivo da análise é obter o padrão de resposta de cada cultivar à adubação com nitrogênio. Os resultados serão apresentados por meio de gráfico. O código está comentado em casa sessão para orientar o leitor. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# dados de índice agronômico de milho de um experimento

da <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/anovareg.txt",
                 header=TRUE, sep="\t")
str(da)

#-----------------------------------------------------------------------------
# estrutura dos dados, experimento balanceado ortogonal
# 4 níveis de bloco (identifica diferenças/variações de ambiente e manejo),
# 3 níveis de cultivar de milho (fator de interesse)
# 6 doses de nitrogênio igualmente espaçadas (fator de interesse)
# 72 observações correspondentes ao produto cartesiano dos níveis dos fatores

#-----------------------------------------------------------------------------
# contraste polinomial, cria matriz de delineamento com colunas ortogonais.
# para dois conjuntos de valores diferentes, o contr.poly é o mesmo
# pois é estabelecido pelo rank e não pelos valores em si.
# só é interpretável para seqüências regulares de valores.
# isso é só demonstração, não iremos usar ortogonais nos modelos reduzidos.

x <- ordered(seq(1,10,1))                     # seqüência regular
y <- ordered(c(1,1.5,2.5,3,3.5,4,4.5,6,8,10)) # seqüência irregular
X <- model.matrix(~x)    # matriz do modelo para níveis de x
Y <- model.matrix(~y)    # matriz do modelo para níveis de y
cbind(X[,2:3], Y[,2:3])  # observe que as matrizes são iguais
round(t(X)%*%X)          # colunas ortogonais implica em matriz de covariância
                         # diagonal para as estimativas dos parâmetros

#-----------------------------------------------------------------------------
# objetivo: ajustar modelos de regressão para cada cultivar em função da dose,
# modelos de regressão polinomial, no máximo de grau 3, fazer a predição livre
# do efeito de bloco pois não há interesse em diferenças entre níveis.

#-----------------------------------------------------------------------------
# cria coluna de fator com níveis ordenados para dose

unique(da$dose)           # doses únicas usadas
da$ds <- ordered(da$dose) # fator ordenado
str(da$ds)                # mostra a ordem dos níveis

#-----------------------------------------------------------------------------
# ajuste do modelo saturado, ou seja, grau máximo permitido para o polinômio
# número de parâmetros estimados para dose é número de níveis - 1

m0 <- lm(indice~bloco+cultivar*ds, data=da)
par(mfrow=c(2,2)); plot(m0); layout(1) # checa as pressuposições

#-----------------------------------------------------------------------------
# ajuste do modelo de regressão linear, anova() verifica falta de ajuste (lof)

m1 <- lm(indice~bloco+cultivar*dose, data=da)
par(mfrow=c(2,2)); plot(m1); layout(1) # gráfico 1,1 aponta que requer termos
anova(m0, m1)                          # compara modelos aninhados, testa lof

#-----------------------------------------------------------------------------
# ajusta o modelo quadrático para dar curvatura, 3 alternativas para isso
# 1) fórmula explícita;
# 2) polinômios ortogonais de grau k=degree;
# 3) polinômios não ortogonais, assim, estimativas iguais em 1 e 3;
# os valores ajustados/preditos não mudam com as alternativas.

m2.1 <- lm(indice~bloco+cultivar*(dose+I(dose^2)), data=da)
m2.2 <- lm(indice~bloco+cultivar*poly(dose, degree=2), data=da)
m2.3 <- lm(indice~bloco+cultivar*poly(dose, degree=2, raw=TRUE), data=da)
sapply(list(m2.1, m2.2, m2.3), coef)   # coeficientes lado a lado
sapply(list(m2.1, m2.2, m2.3), fitted) # valores ajustados lado a lado

m2 <- m2.1
par(mfrow=c(2,2)); plot(m2); layout(1) # ok
anova(m0, m2)                          # modelo quadrático não apresenta lof
anova(m2)                              # teste de hipóteses com o modelo final

#-----------------------------------------------------------------------------
# como fazer predição
# predição com fitted(): faz a predição total, considerando todos os termos do
#                        modelo, tem dimensão dos dados, chama-se de valores
#                        ajustados.
# predição com predict(): faz a predição para quaisquer valores, mas todos os
#                         termos estão presentes.
# predição com produto de matrizes: faz o que você quiser.

#-----------------------------------------------------------------------------
# valores observados (o) e ajustados (f)

of <- cbind(o=da$indice, f=fitted(m2)); of
plot(of)

#-----------------------------------------------------------------------------
# usando a função predict() para prever usando apenas um nível de bloco.
# vou usar o primeiro nível. como bloco ortogonal aos demais efeitos, os
# contrastes entre quaisquer outros níveis dos fatores são os mesmos
# independente do nível do bloco. fazer uma sequencia mais fina para a dose.

bl <- levels(da$bloco)[1] # nível de bloco usando na predição
ct <- levels(da$cultivar) # todos os níveis de cultivar serão usados
ds <- seq(0,300,by=2)     # seqüência fina de valores para dose

pred <- expand.grid(bloco=bl, cultivar=ct, dose=ds) # data.frame da predição
aux <- predict(m2, new=pred, interval="confidence") # IC 95%
str(aux) # matriz com colunas fit, lwr, upr

#-----------------------------------------------------------------------------
# tratamento dos dados para gráfico com bandas usando lattice::xyplot()

pred1 <- cbind(pred, as.data.frame(aux))
str(pred1) # contém apenas o nível I de bloco

require(reshape)

pred1 <- melt(pred1, id=c("cultivar","dose","bloco"))
names(pred1)[ncol(pred1)] <- "indice"
str(pred1)
da$variable <- "obs"
str(da)

pred1 <- rbind(da[,c("cultivar","dose","bloco","variable","indice")], pred1)
str(pred1) # junção dos valores observados e preditos

#-----------------------------------------------------------------------------
# gráfico com observados, preditos e IC

require(lattice)

xyplot(indice~dose|cultivar, groups=variable, data=pred1,
       distribute.type=TRUE, sub="predição para o bloco I",
       type=c("l","l","p","l"), col=c(2,3,1,3), lty=c(1,2,0,2))

#-----------------------------------------------------------------------------
# o gráfico aponta desvio sistemático dos observados para os preditos porque
# nossa predição considerou só o bloco I. se mudarmos de bloco as curvas serão
# movimentadas verticalmente pela mudança do intercepto para o de outro bloco.
# para que o predito passe no meio dos dados é preciso fazer o efeito de bloco
# se anular. se usarmos o contraste soma zero para blocos, podemos saber
# quanto o bloco I difere de zero e sutrair esse valor do fit, lwr e upr

#-----------------------------------------------------------------------------
# usar contraste soma zero e diminuir dos preditos o efeito efeito do bloco I

m3 <- lm(formula(m2), data=da, contrasts=list(bloco=contr.sum))
cbind(coef(m2), coef(m3)) # os coeficientes de bloco mudaram

blI <- coef(m3)["bloco1"]; blI # efeito do bloco I sob a restrição soma zero

# remove o efeito do bloco I dos dados
pred1[pred1$variable!="obs",]$indice <- pred1[pred1$variable!="obs",]$indice-blI

#-----------------------------------------------------------------------------
# faz o gráfico sem o efeito de bloco

xyplot(indice~dose|cultivar, groups=variable, data=pred1,
       distribute.type=TRUE, sub="valores preditos sem efeito de bloco",
       type=c("l","l","p","l"), col=c(2,3,1,3), lty=c(1,2,0,2))

#-----------------------------------------------------------------------------
# predição pelo procedimento matricial. é um método mais geral. é valido para
# modelos da classe glm() com algumas modificações.

# matriz de predição sem colunas para os efeitos de blocos
X <- model.matrix(~cultivar*(dose+I(dose^2)), data=pred)

# vetor de estimativas sem os efeitos de blocos
b <- coef(m3)[-grep("bloco", names(coef(m3)))]

# vetor de erros padrões das estimativas
se <- predict(m3, newdata=pred, se.fit=TRUE)$se

# quantil superior da distribuição t para um IC de 95%
tc <- qt(0.975, df=df.residual(m3))

# margem de erro
me <- se*tc

# data.frame com valores preditos e IC
pred2 <- pred[,-1]
pred2$fit <- X%*%b
pred2$lwr <- pred2$fit-me; pred2$upr <- pred2$fit+me
pred2 <- melt(pred2, id=c("cultivar","dose"))
names(pred2)[ncol(pred2)] <- "indice"

pred2 <- rbind(da[,c("cultivar","dose","variable","indice")], pred2)
str(pred2) # data.frame dos valores preditos pelo método matricial

#-----------------------------------------------------------------------------
# gráfico com observados, preditos e IC pelo método matricial

#png("f018.png", w=500, h=250)
xyplot(indice~dose|cultivar, groups=variable, data=pred2,
       distribute.type=TRUE, layout=c(3,1),
       type=c("l","l","p","l"), col=c(2,3,1,3), lty=c(1,2,0,2),
       strip=strip.custom(bg="gray90"),
       xlab="Dose de nitrogênio (kg/ha)", ylab="Índice agronômico",
       sub=list("predição livre do efeito de bloco", cex=0.8, font=3),
       key=list(space="top", columns=3, type="o", divide=1,
         lines=list(pch=c(1,NA,NA), lty=c(0,1,2), col=c(1,2,3)),
         text=list(c("valores observados","valores preditos","IC 95%"))))
#dev.off()

#-----------------------------------------------------------------------------
# compara os métodos de predição

round(pred1$indice-pred2$indice, 3) # pred1 e pred2 são iguais

#-----------------------------------------------------------------------------
# nesse procedimento consideramos os efeitos de bloco como fixos.
# caso fossem de efeito aleatório poderiamos usar a função nlme::lme() e
# declarar o efeito de blocos como aleatório. os valores preditos seriam
# obtidos com predict(..., level=0). algumas manipulações seriam necessárias
# para obter o intervalo de confiança.

#-----------------------------------------------------------------------------

Como adicionar legendas em gráficos da lattice

Gráfico com bandas de confiança, legenda e assinatura de autor e versão.

No post Como fazer legendas em gráficos eu mostrei procedimentos para gráficos do pacote graphics. Nesse post vou apresentar procedimento para fazer legendas em gráficos da lattice, desenvolvido por Deepayan Sarkar. Iremos empregar as opções auto.key=, key=, colorkey=, e as funções draw.key() e draw.colorkey(). O pacote grid, desenvolvido por Paul Murrell será necessário para posicionarmos as legendas na página gráfica. Nos exemplos temos legendas para pontos, linhas, pontos e linhas, barras, escala de cores para gráficos de três dimensões, gráfico com duas legendas, posicionamento da legenda com o mouse e gráfico com legenda para representar o ajuste de um modelo de regressão com valores observados, preditos e intervalo de confiança.

#-----------------------------------------------------------------------------
# pacotes

require(lattice) # para fazer os gráficos
require(grid)    # para posicionar as legendas

#-----------------------------------------------------------------------------
# dados de indice agronômico em função da cultivar de milho e dose de N

da <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/anovareg.txt",
                 header=TRUE, sep="\t")
str(da)

#-----------------------------------------------------------------------------
# legenda para gráfico de dispersão

# gráfico com legenda padrão
xyplot(indice~dose, groups=cultivar, data=da, jitter.x=TRUE, auto.key=TRUE)

# gráfico com legenda em colunas e alteração dos padrões dos pontos
xyplot(indice~dose, groups=cultivar, data=da, jitter.x=TRUE,
       auto.key=list(columns=3),
       par.settings=list(superpose.symbol=list(pch=15:17, col=3:5, cex=1.2)))

# muda a posição da legenda
xyplot(indice~dose, groups=cultivar, data=da, jitter.x=TRUE,
       auto.key=list(space="right", columns=1),
       par.settings=list(superpose.symbol=list(pch=15:17, col=3:5, cex=1.2)))

# coloca a legenda dentro do gráfico
xyplot(indice~dose, groups=cultivar, data=da, jitter.x=TRUE,
       auto.key=list(corner=c(0.95,0.05), columns=1),
       par.settings=list(superpose.symbol=list(pch=15:17, col=3:5, cex=1.2)))

#-----------------------------------------------------------------------------
# legenda para gráfico de linhas

db <- with(da, aggregate(indice, list(cultivar=cultivar, dose=dose), mean))
col <- trellis.par.get("superpose.line")$col[1:nlevels(db$cultivar)]

# gráfico com legendas de linhas, nomes *antes* das linhas
xyplot(x~dose, groups=cultivar, data=db, type="l",
       key=list(corner=c(0.95,0.05), columns=1,
         text=list(levels(db$cultivar)),
         lines=list(col=col)))

# gráfico com legendas de linhas, nomes *depois* das linhas
xyplot(x~dose, groups=cultivar, data=db, type="l",
       key=list(corner=c(0.95,0.05), columns=1,
         lines=list(col=col),
         text=list(levels(db$cultivar))))

#-----------------------------------------------------------------------------
# legenda para gráfico de linhas e pontos

xyplot(x~dose, groups=cultivar, data=db, type="o", pch=19,
       key=list(corner=c(0.95,0.05), columns=1,
         type="o", divide=1,
         lines=list(col=col, pch=19),
         text=list(levels(db$cultivar))))

#-----------------------------------------------------------------------------
# legenda para pontos e reta de regressão

# ajuste de modelo e criação de objeto com valores observados e preditos
m0 <- lm(indice~cultivar*(dose+I(dose^2)), data=da)
dc <- expand.grid(cultivar=levels(da$cultivar), dose=seq(0,300,l=100))
dc$indice <- predict(m0, newdata=dc)
dc$ind <- "predito"
da$ind <- "obsevado"
dc <- rbind(da[,-3], dc)
str(dc)

# legenda fora do gráfico
xyplot(indice~dose|cultivar, groups=ind, data=dc,
       distribute.type=TRUE, type=c("p","l"), col=1, layout=c(3,1),
       key=list(space="top", columns=2, type="o", divide=1,
         lines=list(pch=c(1,NA), lty=c(0,1)),
         text=list(c("valores observados","valores preditos"))))

# legenda dentro do gráfico
xyplot(indice~dose|cultivar, groups=ind, data=dc,
       distribute.type=TRUE, type=c("p","l"), col=1, layout=c(3,1),
       key=list(corner=c(0.95,0.05), columns=1, type="o", divide=1,
         lines=list(pch=c(1,NA), lty=c(0,1)),
         text=list(c("valores observados","valores preditos"))))

#-----------------------------------------------------------------------------
# dados de mortalidade de aves em função do tipo de resfriamento do aviário

ma <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/mortes.txt",
                 header=TRUE, sep="\t")
str(ma)

#-----------------------------------------------------------------------------
# legenda para gráficos de barras

mb <- with(ma, aggregate(mortes, list(galpao=galpao, sistema=asper), sum))
str(mb)

col <- c("forestgreen", "palegreen")

# legenda fora da região gráfica
barchart(x~galpao, groups=sistema, data=mb, horiz=FALSE, col=col,
         key=list(space="top", columns=2,
           rectangles=list(col=col),
           text=list(c("com aspersão","sem aspersão"))))

# legenda fora da região gráfica com título
barchart(x~galpao, groups=sistema, data=mb, horiz=FALSE, col=col,
         key=list(title="Sistema de resfriamento", cex.title=1.1,
           space="top", columns=2,
           rectangles=list(col=col),
           text=list(c("com aspersão","sem aspersão"))))

# legenda dentro da região gráfica com título
barchart(x~galpao, groups=sistema, data=mb, horiz=FALSE, col=col,
         key=list(title="Sistema de resfriamento", cex.title=1.1,
           corner=c(0.05,0.95), columns=1, padding.text=2,
           rectangles=list(col=col, size=3, height=0.8, border=TRUE),
           text=list(c("com aspersão","sem aspersão"))))

#-----------------------------------------------------------------------------
# escolhendo a posição da legenda com o mouse

barchart(x~galpao, groups=sistema, data=mb, horiz=FALSE, col=col)
trellis.focus(name="page")
loc <- grid.locator(unit="npc") # clicar no gráfico
trellis.unfocus()

barchart(x~galpao, groups=sistema, data=mb, horiz=FALSE, col=col,
         key=list(title="Sistema de\nresfriamento", cex.title=1.1,
           x=loc$x, y=loc$y, corner=c(0.5,0.5), columns=1, padding.text=2,
           rectangles=list(col=col, size=3, height=0.8, border=TRUE),
           text=list(c("com\naspersão","sem\naspersão"))))

#-----------------------------------------------------------------------------
# nota de rodapé no gráfico com assinatura, data e versão

# ajuste de modelo e criação de objeto com valores observados, preditos e IC
m0 <- lm(indice~cultivar*(dose+I(dose^2)), data=da)
dc <- expand.grid(cultivar=levels(da$cultivar), dose=seq(0,300,l=100))
indice <- predict(m0, newdata=dc, interval="confidence")
dc <- cbind(dc, stack(data.frame(indice)))
names(dc)[3] <- "indice"
dc <- rbind(da[,-3], dc)

#png("f017.png", w=500, h=350)
xyplot(indice~dose|cultivar, groups=ind, data=dc,
       distribute.type=TRUE, type=c("l","l","p","l"),
       col=1, lty=c(1,2,0,2), layout=c(3,1),
       strip=strip.custom(bg="gray90"),
       xlab="Dose de nitrogênio (kg/ha)", ylab="Índice agronômico",
       key=list(space="top", columns=3, type="o", divide=1,
         lines=list(pch=c(1,NA,NA), lty=c(0,1,2)),
         text=list(c("valores observados","valores preditos","IC 95%"))))
pushViewport(viewport())
grid.text(label=paste("Walmes Zeviani --",
            format(Sys.time(), "%d/%m/%Y %H:%M:%S --"),
            R.version.string), rot=90,
          x=unit(1, "npc")-unit(2, "mm"),
          y=unit(0, "npc")+unit(2, "mm"),
          just=c("left", "bottom"),
          gp=gpar(cex=0.8, col="gray50"))
popViewport()
#dev.off()

#-----------------------------------------------------------------------------
# legenda para gráficos de 3 dimensões

x <- seq(-10,10,l=50)
dd <- expand.grid(x=x, y=x)
dd$z <- with(dd, 500+x+2*y-0.5*x*y-x^2-y^2)

# gráfico sem cores e sem legenda
wireframe(z~x+y, data=dd)

# gráfico com cores e legenda
wireframe(z~x+y, data=dd, drape=TRUE)

# remove a legenda
wireframe(z~x+y, data=dd, drape=TRUE, colorkey=FALSE)

# coloca a legenda embaixo
wireframe(z~x+y, data=dd, drape=TRUE, colorkey=list(space="bottom"))

# redimensionando a legenda
wireframe(z~x+y, data=dd, drape=TRUE,
          colorkey=list(width=1, height=0.7))

# gráfico de níveis de cor
levelplot(z~x+y, data=dd)

# com contornos e redimensiomento da legenda
levelplot(z~x+y, data=dd, contour=TRUE,
          colorkey=list(width=1, height=0.7))

#-----------------------------------------------------------------------------
# usando a draw.key(), usar fora da função gráfica ou dentro da panel()

xyplot(indice~dose|cultivar, groups=ind, data=dc,
       distribute.type=TRUE, type=c("l","l","p","l"),
       col=1, lty=c(1,2,0,2), layout=c(2,2),
       strip=strip.custom(bg="gray90"),
       xlab="Dose de nitrogênio (kg/ha)", ylab="Índice agronômico")
draw.key(list(space="top", columns=1, type="o", divide=1,
         lines=list(pch=c(1,NA,NA), lty=c(0,1,2)),
         text=list(c("valores observados","valores preditos","IC 95%"))),
         draw=TRUE,
         vp=viewport(x=unit(0.75, "npc"), y=unit(0.65, "npc")))

#-----------------------------------------------------------------------------
# usando a draw.colorkey(), usar fora da função gráfica ou dentro da panel()

wireframe(z~x+y, data=dd, drape=TRUE, colorkey=FALSE,
          panel=function(...){
            panel.wireframe(...)
            draw.colorkey(list(space="bottom", at=seq(250,500,25),
                               height=0.8, width=1), draw=TRUE,
                          vp=viewport(x=unit(0.5,"npc"), y=unit(0.07,"npc")))
          })

#-----------------------------------------------------------------------------
# usando duas legendas

de <- expand.grid(w=gl(3,5,la="w"), z=gl(4,3,la="z"))
de$x <- with(de, rnorm(w, as.numeric(w), 0.5))
de$y <- with(de, rnorm(z, as.numeric(z), 0.5))
pch <- c(1,2,5,7); col=1:3

xy <- xyplot(y~x, data=de, groups=z:w,
             pch=rep(pch, e=3), col=rep(col, 4))

print(xy, position=c(0,0,0.8,1), more=FALSE)
draw.key(list(text=list(levels(de$z)), points=list(pch=pch, col=1),
         space="top", columns=1, title="Níveis de Z", cex.title=1.1),
         draw=TRUE, vp=viewport(x=unit(0.85, "npc"), y=unit(0.6, "npc")))
draw.key(list(text=list(levels(de$w)), points=list(pch=15, col=col),
         space="top", columns=1, title="Níveis de W", cex.title=1.1),
         draw=TRUE, vp=viewport(x=unit(0.85, "npc"), y=unit(0.4, "npc")))

#-----------------------------------------------------------------------------

Método gráfico para valores iniciais em regressão não linear com RStudio

Obtenção de valores iniciais para o modelo duplo van Genuchten de 7 parâmetros com os recursos da função manipulate do RStudio.

No post Método gráfico interativo para valores iniciais em regressão não linear eu apresentei alguns procedimentos gráficos para obter bons valores iniciais para ajuste de modelos de regressão não linear. Ao obter valores iniciais pelo procedimento gráfico facilita-se a convergência do método de estimação e você passa a entender/interpretar melhor os parâmetros do modelo.

Nesse post vou apresentar os métodos gráficos construídos com os recursos do pacote manipulate que é parte do RStudio, um recente editor de scripts R que tem ganhado muitos adeptos por ser muito amigável. Eu já postei o uso dessas funções em RStudio e a função manipulate() onde falo sobre as opções das função.

Vou apresentar o procedimento com dados reais. Preciso encontrar chutes iniciais para o modelo duplo van Genuchten (ver artigo original). Esse modelo tem 7 parâmetros na seguinte equação

\theta(\Psi) = \theta_{res}+\displaystyle \frac{\theta_{pmp}-\theta_{res}}{(1+(\alpha_{tex}\Psi)^{n_{tex}})^{1-1/n_{tex}}}+\displaystyle \frac{\theta_{sat}-\theta_{pmp}}{(1+(\alpha_{est}\Psi)^{n_{est}})^{1-1/n_{est}}},

em que \theta e \Psi são valores observados e o restante são os parâmetros da função. Em geral, quanto maior o número de parâmetros a serem estimados num modelo de regressão mais demorara será a convergência porque o espaço da solução é de alta dimensão. Além do mais, com muitos parâmetros a serem estimados, pode haver fraca de identificabilidade do modelo. Dessa forma, bons chutes iniciais podem ajudar na convergência. Uma ligeira noção do campo de variação de cada parâmetro o usuário precisa ter (limites inferior e superior). Você pode inicialmente usar intervalos amplos e ir reduzindo o comprimento do intervalo a medida que observar quais são os valores mais próximos dos ótimos.

#-----------------------------------------------------------------------------
# dados para a Curva de Retenção de Água, umidade (theta) e tensão (psi)

cra <- data.frame(psi=c(0.01,1,2,4,6,8,10,33,60,100,500,1500,2787,4727,6840,
                    7863,9030,10000,10833,13070,17360,21960,26780,44860,
                    69873,74623,87287,104757,113817,147567,162723,245310,
                    262217,298223),
                  theta=c(0.779,0.554,0.468,0.406,0.373,0.36,0.344,0.309,
                    0.298,0.292,0.254,0.241,0.236,0.233,0.223,0.202,0.172,
                    0.187,0.138,0.098,0.07,0.058,0.052,0.036,0.029,0.0213,
                    0.0178,0.0174,0.0169,0.0137,0.0126,0.0109,0.0106,0.0053))

#-----------------------------------------------------------------------------
# gráfico dos valores observados

plot(theta~log10(psi), data=cra)

#-----------------------------------------------------------------------------
# função do modelo duplo van Genuchten, 7 parâmetros

dvg <- function(x, ts, ti, tr, ae, at, ne, nt){
  tr+(ti-tr)/((1+(at*10^x)^nt)^(1-1/nt))+(ts-ti)/((1+(ae*10^x)^ne)^(1-1/ne))
}

dvg(log10(c(1,10,100,1000,10000)),         # log 10 das tensões
    0.7, 0.2, 0.05, 1.3, 0.0001, 1.5, 3.5) # valor dos parâmetros

#-----------------------------------------------------------------------------
# procedimento gráfico para obter chutes e conhecer a função dos parâmetros

require(manipulate) # carrega pacote que permite manipulação gráfica
start <- list()     # cria uma lista vazia para receber os valores finais

manipulate({
             plot(theta~log10(psi), data=cra)
             curve(dvg(x, ts=ts, ti=ti, tr=tr, ae=ae, at=at, ne=ne, nt=nt), add=TRUE)
             start <<- list(ts=ts, ti=ti, tr=tr, ae=ae, at=at, ne=ne, nt=nt)
           },
           ts=slider(0.7, 0.9, initial=0.8),
           ti=slider(0.15, 0.25, initial=0.2),
           tr=slider(0, 0.10, initial=0.05),
           ae=slider(1.01, 3, initial=1.3),
           at=slider(0, 0.0001, initial=0.00005),
           ne=slider(1.01, 3, initial=1.65),
           nt=slider(1.8, 5, initial=4.3))

#-----------------------------------------------------------------------------

start # valores salvos do último movimento

#-----------------------------------------------------------------------------
# ajuste do modelo aos dados com os chutes do procedimento gráfico

n0 <- nls(theta~tr+(ti-tr)/((1+(at*psi)^nt)^(1-1/nt))+
          (ts-ti)/((1+(ae*psi)^ne)^(1-1/ne)),
          data=cra, start=start, trace=TRUE)
summary(n0) # quadro de estimativas

#-----------------------------------------------------------------------------
# gráfico dos dados com a curva estimada

lis <- c(list(x=NULL), as.list(coef(n0)), body(dvg))
plot(theta~log10(psi), cra,
     ylab=expression(Conteúdo~de~água~no~solo~(theta*","~g~g^{-1})),
     xlab=expression(Tensão~matricial~(Psi*","~kPa)), xaxt="n")
tmp <- as.function(lis)
curve(tmp, add=TRUE, col=2, lwd=1.5)
axis(1, at=-2:5, label=as.character(10^(-2:5)), lwd.ticks=2)
s <- log10(sort(sapply(1:9, function(x) x*10^(-3:6))))
axis(1, at=s, label=NA, lwd=0.5)
abline(v=-2:6, h=seq(0,1,0.05), lty=3, col="gray50")

#-----------------------------------------------------------------------------

O procedimento permite encontrar chutes iniciais para duas ou mais curvas. Basta escrever adequadamente as funções dentro da manipulate(). Abaixo estão os códigos para obter os valores inicias para o modelo logístico usado no contexto de epidemiologia. Em duas variedades de uma espécie vegetal foram observadas a severidade da doença ao longo do tempo. O objetivo é encontrar estimativas dos parâmetros para a função (que no caso é a logística) que associa severidade ao tempo e comparar as funções entre as variedades. A comparação das funções pode ser feita fazendo-se testes de hipóteses para funções dos parâmetros (em geral contrastes) ou testes envolvendo modelos aninhados, como você verá a seguir.

#-----------------------------------------------------------------------------
# dados de severidade em função do tempo para duas variedades

da <- read.table("http://www.leg.ufpr.br/~walmes/ridiculas/epidem.txt",
                 header=TRUE, sep="\t")
str(da)

#-----------------------------------------------------------------------------
# gráficos exploratórios

require(lattice)
xyplot(V1+V2~DAI, data=da, type=c("p","a"))

# ajustar modelos. separado ou junto?
# qual modelo? quais os chutes?

#-----------------------------------------------------------------------------
# mudando a disposição dos dados

require(reshape)
db <- melt(da, id=c("DAI","rep"))
names(db)[3:4] <- c("vari","sever")
str(db)

#-----------------------------------------------------------------------------
# carrega pacote e faz a função para encontrar via deslizadores valores ótimos

require(manipulate)
start <- list()

manipulate({
  plot(sever~DAI, db, col=ifelse(db$vari=="V1",1,2))
  A1 <- a1; S1 <- s1; X01 <- x01
  A2 <- a2; S2 <- s2; X02 <- x02
  curve(A1/(1+exp(-(x-X01)/S1)), add=TRUE, col=1)
  curve(A2/(1+exp(-(x-X02)/S2)), add=TRUE, col=2)
  start <<- list(A=c(A1,A2), S=c(S1,S2), X0=c(X01,X02))
  },
  a1=slider(90, 110, initial=95, step=1),
  a2=slider(90, 110, initial=95, step=1),
  s1=slider(2, 8, initial=6, step=0.1),
  s2=slider(2, 8, initial=6, step=0.1),
  x01=slider(10, 40, initial=20, step=1),
  x02=slider(10, 40, initial=20, step=1)
  )

start # valores do último movimento

#-----------------------------------------------------------------------------
# ajuste do modelo não restrito: cada variedade com os seus parâmetros

n0 <- nls(sever~A[vari]/(1+exp(-(DAI-X0[vari])/S[vari])),
          data=db, start=start)
summary(n0)

#-----------------------------------------------------------------------------
# ajuste do modelo restrito: as variades possuem mesmo parâmetro S

start$S <- mean(start$S)
start

n1 <- nls(sever~A[vari]/(1+exp(-(DAI-X0[vari])/S)),
          data=db, start=start)
summary(n1)

#-----------------------------------------------------------------------------
# testa a hipótese do modelo nulo n1 tendo como alternativa o modelo n0

anova(n1, n0)

#-----------------------------------------------------------------------------

Obtenção de valores iniciais para o ajuste de duas curvas.

Ainda é necessário fazer a checagem das pressuposições envolvidas no modelo de regressão clássico: homocedasticidade (a variância dos desvios de regressão é constante ao longo a curva), normalidade (os erros aleatórios têm distribuição normal) e independência entre observações (normalmente a independência é garantida por um plano de amostragem adequado). Consulte Análise de resíduos em regressão não linear para saber como obter os 4 gráficos de análise de resíduos.

Esse post foi um aprimoramento de uma resposta na R-help e estimulado por citações na lista. Até a próxima ridícula.

Gráfico com escala logarítmica em notação de potência

Gráfico (verde) em escala logarítmica com notação de potência.

O pacote latticeExtra traz funções que permitem aprimorar as funções gráficas disponíveis no pacote lattice, ambos desenvolvidos por Deepayan Sarkar. Recentemente precisei confeccionar gráficos em escala logarítmica para os quais eu tinha um enorme roteiro de funções construídas à custo de muita dedicação e tempo. Felizmente, hoje não preciso mais reaprender as implementações que fiz porque a latticeExtra fornece funções para escrever as escalas dos eixos em notação de potência.

As funções disponibilizadas são: xscale.components.logpower, xscale.components.fractions, xscale.components.log10ticks, xscale.components.log10.3 e
xscale.components.subticks, todas elas na versão para o eixo y também. Essas funções permitem, por exemplo, eixos escritos com potência e fração e os traços marcas (tick) para datas considerando o diferente número de dias de cada mês. Tudo isso você pode conferir rodando os exemplos das funções. Até a próxima ridícula.

#-----------------------------------------------------------------------------

require(lattice)

data(Earthquake, package="MEMSS")

p1 <-
  xyplot(accel~distance, data=Earthquake, pch=19, lwd=2,
         prepanel=prepanel.loess, col="blue",
         type=c("p", "g", "smooth"), scales=list(log=2),
         xlab="Distância do epicentro (km)",
         ylab="Aceleração máxima horizontal (g)")

require(latticeExtra)

p2 <-
  xyplot(accel~distance, data=Earthquake, pch=19, lwd=2,
         prepanel=prepanel.loess, col="forestgreen",
         type=c("p", "g", "smooth"), scales=list(log=2),
         xlab="Distância do epicentro (km)",
         ylab="Aceleração máxima horizontal (g)",
         xscale.components=xscale.components.logpower,
         yscale.components=yscale.components.logpower)

#png("f011.png", w=500, h=300)
print(p1, split=c(1,1,2,1), more=TRUE)
print(p2, split=c(2,1,2,1))
#dev.off()

#-----------------------------------------------------------------------------