Página Inicial > delineamento, gráficos, regressão > Predição em modelos de regressão com blocos

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.

#-----------------------------------------------------------------------------
About these ads

Deixe um comentário

Preencha os seus dados abaixo ou clique em um ícone para log in:

Logotipo do WordPress.com

Você está comentando utilizando sua conta WordPress.com. Sair / Alterar )

Imagem do Twitter

Você está comentando utilizando sua conta Twitter. Sair / Alterar )

Foto do Facebook

Você está comentando utilizando sua conta Facebook. Sair / Alterar )

Foto do Google+

Você está comentando utilizando sua conta Google+. Sair / Alterar )

Conectando a %s

Seguir

Obtenha todo post novo entregue na sua caixa de entrada.

Junte-se a 52 outros seguidores

%d blogueiros gostam disto: