Como fazer regressão polinomial com diferentes graus

July 18, 2011
By

This post was kindly contributed by Ridículas - go there to comment and to read the full post.

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()

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

Tags: , , , , , , , , , ,

Comments are closed.