Calcular la moda de una regresión cuadrática con intervalos de confianza

Question

Tengo una regresión cuadrática y contra x y me interesa el valor x donde y es el máximo (ymax->x). Puedo calcular x(ymax) pero también me interesa el error estándar o el intervalo de confianza de x(ymax) para obtener x(ymax) +/- error estándar o, respectivamente, intervalo de confianza inferior y superior de xmax.

Una solución podría ser utilizar métodos de bootstrap y calcular nuevas estimaciones de la beta, pero los datos reales que tengo son grandes y se necesita una hora para calcular el modelo (modelo mixto con estructura de varianza-covarianza).

Otra posibilidad sería utilizar los errores estándar de las betas para calcular una banda de confianza y encontrar los puntos de corte de la horizontal a través del máximo de la curva y la "curva de confianza" superior (véase más abajo). Pero no estoy seguro de que este método sea una forma válida de obtener el intervalo de confianza.

¿Hay alguien que pueda darme algún consejo sobre cómo afrontar esto?

Gracias de antemano.

---

# simulate data set
b0 <- 100
b1 <- 100
b2 <- -2
(x <- rep(1:50,2))
(y <- b0 + b1*x + b2*x^2 + rnorm(length(x),0,200))

# fit model
fit <- lm(y~x + I(x^2))
summary(fit)

# visualize
plot(y~x)
lines(fitted(fit), col="red")

# retrieve coefficients
beta <- fit$coef

# build function for optimization
f1 <- function(x) sum(beta * c(1,x,x^2))

# compute xmax
(xymax <- optimize(f1,interval=c(0,50), maximum=TRUE))

# add line at xmax
# add line at xmax
abline(v=xymax$max, col="blue")
abline(h=xymax$obj, col="blue")

# compute confidence band
cilbeta <- beta - qnorm(0.975,0,1)*summary(fit)$coef[,"Std. Error"]
ciubeta <- beta + qnorm(0.975,0,1)*summary(fit)$coef[,"Std. Error"]
f2 <- function(x) { sum(c(1,x,x^2)*cilbeta) }
f3 <- function(x) { sum(c(1,x,x^2)*ciubeta) }
# curve(f2,from=1,to=50)  does not work !!
ycil = c()
yciu = c()
for (i in 1:50) {
    ycil <- c(ycil,f2(i))
    yciu <- c(yciu,f3(i))
}
lines(1:50,ycil, lty=2)
lines(1:50,yciu, lty=2)

# cutting points between upper "confidence curve" horizontal line
# through y.max

# x.max = -b2/(2b3) = -0.5*b2/b3
# var(x.max) = 0.5^2*var(b2/b3)
# var.b23 = b2^2/b3^2*[var(b2)/b2^2 -2*cov(b2,b3)/(b2*b3) + var(b3)/b3^2]
(var.b23 <- coef0[2]^2/coef0[3]^2*(vcov0[2,2]/coef0[2]^2 
            -2*vcov0[2,3]/(coef0[2]*coef0[3]) + vcov0[3,3]/coef0[3]^2
            )
        )
(var.xmax <- (-0.5)^2*var.b23)
(se.xmax <- sqrt(var.xmax))
# confidence interval
paste(formatC(x.max,digits=2,format="f")," ("
      , formatC(x.max-2*se.xmax,digits=2, format="f"),"-"
      , formatC(x.max+2*se.xmax,digits=2,format="f"),")"
      , sep=""
      )
# [1] "24.77 (24.12-25.42)"

Answer 1

El modo de la función cuadrática $$ b_1 + b_2 x + b_3 x^2 $$ está en $M = -b_2/(2 b_3)$ siempre y cuando $b_3$ es negativo. (Compare con la solución en su código).

La desviación estándar de M puede obtenerse por bootstrap o analíticamente:

En http://www.stat.cmu.edu/~hseltman/files/ratio.pdf se puede encontrar una aproximación de la varianza de un cociente de dos variables aleatorias correlacionadas. (Obsérvese que la varianza exacta parece ser intrínsecamente complicada, incluso en el caso normal bivariante). Así que basta con sustituir las expectativas, varianzas y covarianzas en la última fórmula del artículo por las estimaciones de la regresión.