Arquivo

Posts Tagged ‘polinômio’

Reparametrização do modelo quadrático para ponto crítico

Perfil de log-verossimilhaça representada pelo módulo da raíz da função deviance para cada um dos parâmetros do modelo quadrático reparametrizado. Linhas tracejadas indicam os limites dos intervalos de confiança de 99, 95, 90, 80 e 50%. Padrões de perfil diferentes de um V indicam intervalos assimétricos.

O emprego de funções polinomiais é muito frequente em ciências agrárias. O polinômio de segundo grau ou modelo quadrático é a escolha mais natural e imediata de modelo para quando a variável dependente apresentam uma tendência/relação não linear com a variável independente. Não raramente, o emprego desse modelo motiva também a estimação do ponto crítico (mínimo ou máximo). Para isso, o que se faz é estimar os parâmetros do modelo e via regras do cálculo obter a estimativa pontual do ponto crítico. O problema é que normalmente se necessita e não se faz, por desconhecimento de métodos, associação de incerteza à essa estimativa. É sobre isso que vou tratar agora: como associar incerteza ao ponto crítico de um modelo quadrático.

O modelo quadrático tem essa equação

m(x) = a+b\cdot x + c\cdot x^2,

em que a é o intercepto, ou seja, E(Y|x=0) = ab é valor da derivada do modelo na origem, ou seja, d(m_0(x))/dx = b, e c mede a curvatura. Derivando o modelo, igualando a zero e resolvendo, obtemos que o ponto crítico é

x_m = -b/(2\cdot c),

portanto, x_m é uma função das estimativas dos parâmetros, e como tal, deve apresentar incerteza associada.

Uma forma de associar incerteza é empregar o método delta para encontrar erros padrões. No R, o mesmo está disponível em msm::deltamethod(). Outra forma é usar um simulação boostratp e aplicar nas amostras bootstrap a relação acima. E o terceiro, e que eu mais gosto pois tem uma série de vantagens, é usar o perfil de log-verossimilhança.

Para usar o perfil de verossimilhança eu tive que escrever o modelo de forma que o ponto crítico fosse um parâmetro dele. Reparametrizar um modelo nem sempre é uma tarefa fácil. Nesse caso o procedimento escrito é curto e simples de entender. Pórem eu demorei cerca de uma hora com tentativa-erro até usar a seguinte abordagem: imagine uma função que tenha ponto de máximo. Podemos aproximar essa função ao redor de um valor x_0 por série de Taylor. Para uma aproximação de segundo grau temos

f(x) \approx f(x_0)+f'(x_0)(x-x_0)+0.5\cdot f''(x_0)(x-x_0)^2.

Agora imagine que x_0 é um ponto crítico da função. O termo linear da aproximação é 0 no ponto crítico por definição. Se a função que estamos aproximando for perfeitamente quadrática, esse procedimento nos leva a reparametrização que desejamos. Assim, o modelo quadrático reparametrizado é

n(x) = y_m+c\cdot (x-x_m)^2.

em que x_m é o nosso parâmetro de interesse, o ponto crítico, y_m é o valor da função no ponto crítico, ou seja, E(Y|x=x_m) = y_m, c é a métade da segunda derivada da função no ponto crítico, ou seja, c=0.5\cdot f''(x_m), é uma curvatura. Então nos reparametrizamos o model quadrático para ter o ponto crítico como parâmetro. Nessa reparametrização obtvemos 3 parâmetros com interpretação simples. Em termos de interpretação, se c=0 não temos curvatura, portanto, não temos ponto crítico e devemos ajustar um modelo linear da reta. Se temos curvatura podemos estimar o ponto crítico e o valor da função no ponto crítico. Então imagine que você estuda doses de nutrientes para uma cultura, não é fundamental saber qual a dose que confere o máximo e qual a produção esperada nessa dose? Melhor que isso é poder associar incerteza à essas quantidades.

O modelo reparametrizado, diferente do original, é um modelo não linear. Por isso, para estimar os parâmetros no R usaremos a função nls() e não a lm(). (Visite os posts passados para ver mais coisas relacionadas à modelos não lineares.) Para obter os intervalos de perfil de log-verossimilhança usaremos os métodos confint() e profile(). Consulte a documentação dessas funções para saber o que elas de fato fazem.

No CMR abaixo eu ilustro a aplicação desses métodos. É recomendando ao leitor que compare os intervalos de confiança, que repita os procedimento simulando os dados com outros parâmetros, tamanhos e espaços, e se for o caso, que use dados reais. Futuros artigos do Rídiculas vão incluir comparação de modelos quadráticos (teste de hipótese) e o modelo quadrático para duas regressoras (bidimensional). Até a próxima rídicula.

#-----------------------------------------------------------------------------
# dados simulados

set.seed(23456)
x <- 1:12
y <- 3+0.17*x-0.01*x^2+rnorm(x,0,0.05)
plot(x, y)                         # dados observados
curve(3+0.17*x-0.01*x^2, add=TRUE) # curva que gerou os dados

#-----------------------------------------------------------------------------
# ajuste do y = a+b*x+c*x², é um modelo linear

m0 <- lm(y~x+I(x^2))
summary(m0)

#-----------------------------------------------------------------------------
# método delta

require(msm)
coef <- coef(m0)[2:3]; vcov <- vcov(m0)[2:3,2:3]
xm <- -0.5*coef[1]/coef[2]; xm                         # estimativa
sd.xm <- deltamethod(~(-0.5*x1/x2), coef, vcov); sd.xm # erro padrão
xm+c(-1,1)*qt(0.975, df=df.residual(m0))*sd.xm         # IC 95%

#-----------------------------------------------------------------------------
# simulação (bootstrap paramétrico)

require(MASS)
aa <- mvrnorm(1000, mu=coef, Sigma=vcov)
xm.boot <- apply(aa, 1, function(i) -0.5*i[1]/i[2])
mean(xm.boot)                      # média bootstrap
sd(xm.boot)                        # desvio-padrão bootstrap
quantile(xm.boot, c(0.025, 0.975)) # IC boostrap

#-----------------------------------------------------------------------------
# ajuste do y = ym+c*(x-xm)², é um modelo não linear

plot(x, y)
abline(h=3.7, v=8.5, lty=3) # aproximação visual do chute
start <- list(ym=3.7, xm=8.5, c=coef(m0)[3])
max(fitted(m0))             # boa aproximação de chute
x[which.max(fitted(m0))]    # boa aproximação de chute
n0 <- nls(y~ym+c*(x-xm)^2, start=start)
summary(n0) # estimativa do xm que é parâmetro do modelo, estimado diretamente

#-----------------------------------------------------------------------------
# gráfico dos observados com a curva ajustada

plot(x, y)
with(as.list(coef(n0)),{
  curve(ym+c*(x-xm)^2, add=TRUE)
  abline(h=ym, v=xm, lty=3)
})

#-----------------------------------------------------------------------------
# veja a igualdade dos modelos

cbind(coef(m0),
      with(as.list(coef(n0)), c(ym+c*xm^2, -2*c*xm, c)))
deviance(m0)
deviance(n0)
plot(fitted(m0), fitted(n0)); abline(0,1)

#-----------------------------------------------------------------------------
# intervalos de confiança

confint.default(n0) # IC assintóticos, são simétricos por construção
confint(n0)         # IC baseado em perfil de log-verossimilhança

#-----------------------------------------------------------------------------
# gráficos de perfil de verossimilhança

#png(file="f030.png", width=500, height=175)
par(mfrow=c(1,3))
plot(profile(n0))
layout(1)
#dev.off()

#-----------------------------------------------------------------------------
# representação do ajuste

pred <- data.frame(x=seq(0,13,l=30))
pred <- cbind(as.data.frame(predict(m0, newdata=pred, interval="confidence")),
              pred)
str(pred)
ic <- confint(n0); ic

plot(x, y, type="n")              # janela vazia
with(pred,                        # bandas de confiança
     polygon(c(x,rev(x)), c(lwr,rev(upr)), col="gray90", border=NA))
points(x, y)                      # valores observados
with(pred, lines(x, fit))         # valores preditor
abline(v=ic[2,], h=ic[1,], lty=2) # IC para xm e ym

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

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.

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