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

April 16, 2012
By

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

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-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 é o dobro da segunda derivada da função no ponto crítico, ou seja, c=2\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. Para estimar os parâmetros no R usaremos a função nls() e não a lm(). Para obter os intervalos de perfil de log-verossimilhança usaremos os métodos confint() e profile().

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 e se for o caso, que use dados reais. Futuros artigos do Rídiculas vão incluir comparação de modelos quadráticos 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

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


Tags: , , ,

Comments are closed.