Arquivo

Posts Tagged ‘reshape’

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.

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

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.

Converter dados em formato amplo e formato longo

Respostas observadas em função dos fatores em experimento com a cultura da soja.

A coisa mais comum para quem analisa dados e ter que editar o formato dos dados. Dificilmente você recebe os dados no formato correto. Muitas pessoas fazem de suas planilhas de dados uma arte, com cédulas coloridas, linhas e colunas com textos em diversos formatos.

Dados de medidas repetidas tem a característica de cada nova medida repetida ser uma nova coluna na planilha (formato amplo). A maioria dos métodos de análise requer que os dados estejam no formato longo, ou seja, cada linha é um registro, cada coluna é uma variável. Portanto, temos que converter o formato dos dados. Isso vai se tornar claro/simples com os comandos que apresento nessa ridícula.

Os pacotes reshape e plyr contém funções úteis para a manipulação/restruturação dos dados. Os comandos abaixo importam dados no formato amplo, passam para o formato longo e retornam para o formato amplo.

> #-----------------------------------------------------------------------------
> # dados de espessura de fibra de Schinus
> 
> require(reshape)
> 
> # dados estão no formato amplo, colunas são secções no tronco da árvore
> fibra <- read.table("http://www.leg.ufpr.br/~walmes/ridiculas/fibra.txt",
+                     header=TRUE, sep="\t")
> str(fibra)
'data.frame':	150 obs. of  7 variables:
 $ arvore : int  1 1 1 1 1 1 1 1 1 1 ...
 $ amostra: int  1 2 3 4 5 6 7 8 9 10 ...
 $ A      : num  2.61 2.43 2.09 2.5 2.06 ...
 $ B      : num  3.17 2.68 2.43 2 1.85 ...
 $ C      : num  2.17 4.19 4.82 5.54 3.26 ...
 $ D      : num  3.7 2.2 3.5 2.61 3.43 ...
 $ E      : num  2.4 3.5 3.11 2.4 3.26 ...
> 
> # passamos os dados para o formato longo, criou uma coluna para secção
> fibra <- melt(fibra, id=c("arvore", "amostra"))
> str(fibra)
'data.frame':	750 obs. of  4 variables:
 $ arvore  : int  1 1 1 1 1 1 1 1 1 1 ...
 $ amostra : int  1 2 3 4 5 6 7 8 9 10 ...
 $ variable: Factor w/ 5 levels "A","B","C","D",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ value   : num  2.61 2.43 2.09 2.5 2.06 ...
> 
> # retornamos os dados para o formato amplo
> fibra <- adply(cast(fibra, amostra~arvore~variable,
+                     value="value", fun=NULL), c(1,2))
> str(fibra)
'data.frame':	150 obs. of  7 variables:
 $ amostra: chr  "1" "2" "3" "4" ...
 $ arvore : chr  "1" "1" "1" "1" ...
 $ A      : num  2.61 2.43 2.09 2.5 2.06 ...
 $ B      : num  3.17 2.68 2.43 2 1.85 ...
 $ C      : num  2.17 4.19 4.82 5.54 3.26 ...
 $ D      : num  3.7 2.2 3.5 2.61 3.43 ...
 $ E      : num  2.4 3.5 3.11 2.4 3.26 ...
> 
> #-----------------------------------------------------------------------------

Os dados do experimento em soja já estão no formato longo, exigido para fazer ajuste de modelos lineares. Nesse experimento foram observadas quatro variáveis resposta. Para fazer uma análise exploratória conjunta, eu passei os dados para o formato longo. Assim, tenho uma coluna que identifica a resposta observada e outra o valor. Com os recursos do pacote lattice eu fiz um gráfico de dispersão para tirar informações iniciais dos resultados do experimento.

#-----------------------------------------------------------------------------
# dados

soja <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/soja.txt",
                   header=TRUE, sep="\t", dec=",")
names(soja)[4:7] <- c("Rendimento de grãos (g/vaso)",
                      "Massa de 100 grãos (g)",
                      "K no grão (g/kg)", "P no grão (g/kg)")
str(soja)

#-----------------------------------------------------------------------------
# criando uma coluna que identifica as variáveis resposta

soja2 <- melt(soja, id=1:3)
str(soja2)

#-----------------------------------------------------------------------------
# gráfico para análise exploratória

require(lattice)

#png("f006.png", w=500, h=400)
mycols <- colors()[c(47,81,614)]
trellis.par.set(superpose.symbol=list(col=mycols, pch=c(15,17,19)),
                superpose.line=list(col=mycols),
                strip.background=list(col="gray90"))
xyplot(value~potassio|variable, groups=agua, data=soja2,
       xlab="Potássio no solo (mg/dm³)", ylab="Valor observado",
       key=list(title="Níveis de umidade do solo (%)", cex.title=1,
         columns=3, type="o", divide=1,
         lines=list(pch=trellis.par.get("superpose.symbol")$pch,
           col=trellis.par.get("superpose.symbol")$col),
         text=list(c(format(unique(soja$agua))))),
       type=c("p","a","g"), jitter.x=TRUE,
       scales=list(y=list(relation="free")), layout=c(2,2))
#dev.off()

#-----------------------------------------------------------------------------
# para ver as configurações gráficas da lattice

library(TeachingDemos)
TkListView(trellis.par.get())

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

Fatorial com desdobramento da interação, gráfico e tabela de médias

Rendimento de grãos em função dos níveis de potássio e umidade do solo. Letras minúsculas comparam níveis de potassio em um mesmo nível de umidade, as letras maiúsculas fazem o inverso.

O R não é um aplicativo específico de análise de experimentos. Mesmo assim, ele oferece funções que permitem a análise dos mais diversos tipos de experimentos, ficando a cargo do usuário saber usar corretamente as funções/pacotes que o R oferece.

Pode-se dizer que os experimentos fatoriais são os mais comuns na ciência. Nas ciências agrárias são amplamente usados. Têm-se uma tradição de representar a interação entre os fatores (qualitativos) empregando testes de médias para níveis de um fator fixando níveis dos demais.

Nessa ridícula, apresento procedimentos para obter o desdobramento de interação em um fatorial duplo. Os dados são do rendimento de grãos de soja em função dos fatores umidade do solo e dose de fertilizante potássico. O experimento foi instalado em casa de vegetação com as unidades experimentais (vasos) agrupadas em blocos que controlaram variações de luminosidade/umidade/ventilação.

Os procedimentos incluem o ajuste do modelo, o teste sobre os efeitos dos fatores pela análise de variância, a checagem das pressuposições do modelo, o fatiamento das somas de quadrados, o teste de médias. No final são apresentados procedimentos para representar os resultados em gráfico e tabela de dupla entrada. É importante dizer que os fatores deveriam ser analisados como quantitativos, onde espera-se uma interpretação mais interessante sobre o fenômeno. Entretanto, o método abaixo é ilustrativo do procedimento para fatores qualitativos. Em uma futura oportunidade analisarei dos dados considerando os fatores como quantitativos.

#-----------------------------------------------------------------------------
# importa e prepara os dados
 
rend <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/rendimento.txt",
                   header=TRUE, sep="\t")
rend <- transform(rend, K=factor(K), A=factor(A), bloc=factor(bloc))
str(rend)
 
#-----------------------------------------------------------------------------
# análise gráfica
 
require(lattice)
xyplot(rg~K|A, groups=bloc, data=rend, type="b", auto.key=TRUE)
 
#-----------------------------------------------------------------------------
# ajuste do modelo e faz checagem
 
m0 <- aov(rg~bloc+A*K, data=rend)
summary(m0)
par(mfrow=c(2,2)); plot(m0); layout(1)
 
#-----------------------------------------------------------------------------
# desdobrando somas de quadrados para a variação de K dentro de A e vice-versa

m1 <- aov(rg~bloc+K/A, data=rend)
desAinK <- sapply(paste("K", levels(rend$K), sep=""), simplify=FALSE,
                  grep, x=names(coef(m1)[m1$assign==3]))
summary(m1, split=list("K:A"=desAinK))

m2 <- aov(rg~bloc+A/K, data=rend)
desKinA <- sapply(paste("A", levels(rend$A), sep=""), simplify=FALSE,
                  grep, x=names(coef(m2)[m2$assign==3]))
summary(m2, split=list("A:K"=desKinA))
 
#-----------------------------------------------------------------------------
# desdobrando a interação em testes de Tukey
 
require(agricolae)
require(plyr)

KinA <- sapply(levels(rend$A), simplify=FALSE,
               function(a){
                 with(subset(rend, A==a),
                      HSD.test(rg, K,
                               DFerror=df.residual(m0),
                               MSerror=deviance(m0)/df.residual(m0)))
               })

KinA <- ldply(KinA, NULL)
KinA$M <- gsub(" ", "", KinA$M, fixed=TRUE)
KinA$trt <- as.factor(as.numeric(as.character(KinA$trt)))
str(KinA)

AinK <- sapply(levels(rend$K), simplify=FALSE,
               function(k){
                 with(subset(rend, K==k),
                      HSD.test(rg, A,
                               DFerror=df.residual(m0),
                               MSerror=deviance(m0)/df.residual(m0)))
               })

AinK <- ldply(AinK, NULL)
AinK$M <- toupper(gsub(" ", "", AinK$M, fixed=TRUE))
AinK$trt <- as.factor(as.numeric(as.character(AinK$trt)))
str(AinK)

#-----------------------------------------------------------------------------
# combinando os data.frames para obter o gráfico e a tabela

KinA$comb <- paste(KinA$trt, KinA$.id)
AinK$comb <- paste(AinK$.id, AinK$trt)
desd <- merge(KinA, AinK, by=c("comb","comb"))
desd$let <- with(desd, paste(M.y, M.x, sep=""))
desd <- desd[,c("trt.x","trt.y","means.x","let")]
names(desd) <- c("K","A","means","let")
str(desd)

#-----------------------------------------------------------------------------
# gráfico das médias com as letras sobre as barras

#png("f005.png", w=500, h=250)
barchart(means~K|A, data=desd, horiz=FALSE, layout=c(3,1),
         ylab="Rendimento de grãos (g/vaso)", ylim=c(NA,45),
         xlab="Dose de potássio (mg/dm³)",
         key=list(text=list("Níveis de umidade do solo (%)")),
         sub=list("*letras minúsculas comparam médias em um mesmo nível de umidade.", cex=0.8, font=3),
         strip=strip.custom(bg="gray90"), col=colors()[51],
         panel=function(x, y, subscripts, ...){
           panel.barchart(x, y, subscripts=subscripts, ...)
           panel.text(x, y, label=desd[subscripts,"let"], pos=3)
         })
#dev.off()

#-----------------------------------------------------------------------------
# tabela de dupla entrada com as médias das combinações A e K

require(reshape)
desd$means.let <- with(desd, paste(means, let))
cast(desd, K~A, value="means.let")

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