Recuperación de coeficientes y varianzas brutos a partir de una regresión polinómica ortogonal

Question

Parece que si tengo un modelo de regresión como $y_i \sim \beta_0 + \beta_1 x_i+\beta_2 x_i^2 +\beta_3 x_i^3$ Puedo ajustar un polinomio bruto y obtener resultados poco fiables o ajustar un polinomio ortogonal y obtener coeficientes que no tienen una interpretación física directa (por ejemplo, no puedo utilizarlos para encontrar las ubicaciones de los extremos en la escala original). Parece que debería ser capaz de tener lo mejor de ambos mundos y ser capaz de transformar los coeficientes ortogonales ajustados y sus varianzas de nuevo a la escala en bruto. He tomado un curso de postgrado en regresión lineal aplicada (utilizando Kutner, 5ed) y miré a través del capítulo de regresión polinomial en Draper (3ed, referido por Kutner), pero no encontré ninguna discusión de cómo hacer esto. El texto de ayuda de poly() en R no lo hace. Tampoco he encontrado nada en mi búsqueda en Internet, incluido aquí. ¿Reconstruir coeficientes brutos (y obtener sus varianzas) a partir de coeficientes ajustados a un polinomio ortogonal...

1. imposible de hacer y estoy perdiendo el tiempo.
2. tal vez sea posible, pero no se sabe cómo en el caso general.
3. posible pero no se discute porque "¿quién querría hacerlo?".
4. posible pero no se discute porque "es obvio".

Si la respuesta es 3 o 4, estaría muy agradecido si alguien tuviera la paciencia de explicar cómo hacerlo o señalar una fuente que lo haga. Si es 1 ó 2, seguiría teniendo curiosidad por saber cuál es el obstáculo. Muchas gracias por leer esto, y pido disculpas de antemano si estoy pasando por alto algo obvio.

Answer 1

Sólo una adición potencialmente útil a respuesta de whuber . Mirando el código de poly puedes deducir el mapa lineal tú mismo. Sea $\vec h_{m:n} =(h_m, h_{m + 1}, \dots, h_n)^\top$ , que los índices negativos sean cero por definición, y los indefinidos $\Gamma$ sean cero. Entonces podemos encontrar que si despreciamos la escala aquí entonces el mapa al polinomio ortogonal viene dado por

$$\begin{align*} z_0 &= 1 = \gamma_{0,0:0}\cdot 1\\ z_1 &= x - \alpha_1 \\ &= \underbrace{(\vec\gamma_{0,-1:0} - \alpha_1 \vec\gamma_{0,0:1})^\top}_{\vec\gamma_{1,0:1}^\top}(1 , x) \\ z_2 &= (x - \alpha_2)z_1 - \frac{\sigma_2}{\sigma_1} z_0 \\ &= x^2 + (\gamma_{10} -\alpha_2)x - \alpha_2\gamma_{10} - \frac{\sigma_2}{\sigma_1} \\ &= (1, x, x^2)\underbrace{ (\vec\gamma_{1,-1:1}-\alpha_2\vec\gamma_{1,0:2} -\frac{\sigma_2}{\sigma_1}\vec\gamma_{0,0:2})}_{ \vec\gamma_{2,0:2}^\top}\\ z_3 &= (x - \alpha_3)z_2 - \frac{\sigma_3}{\sigma_2} z_1 \\ &= x^3 + (\gamma_{21}-\alpha_3)x^2 + (\gamma_{20}-\alpha_3\gamma_{21})x -\alpha_3\gamma_{20} -\frac{\sigma_3}{\sigma_2} x - \frac{\sigma_3}{\sigma_2}\gamma_{01} \\ &= (1, x, x^2, x^3)\underbrace{(\vec\gamma_{2,-1:2}-\alpha_3\vec\gamma_{2,0:3} -\frac{\sigma_3}{\sigma_2}\vec\gamma_{1,0:3})}_{ \vec\gamma_{3,0:3}}\\ \vdots\, &= \,\vdots \end{align*}$$

Así, podemos calcular el $\Gamma$ matriz con este código

get_poly_orth_map <- function(object){
  stopifnot(inherits(object, "poly"))
  sigs  <- attr(object, "coefs")$norm2
  alpha <- attr(object, "coefs")$alpha

  nc <- length(alpha) + 1L
  Gamma <- matrix(0., nc, nc)
  Gamma[1, 1] <- 1
  if(nc > 1){
    Gamma[ , 2] <- -alpha[1] * Gamma[, 1]
    Gamma[2, 2] <- 1
  }
  if(nc > 2)
    for(i in 3:nc){
      i_m1 <- i - 1L
      Gamma[, i] <- c(0, Gamma[-nc, i_m1]) - alpha[i_m1] * Gamma[, i_m1] - 
        sigs[i] / sigs[i_m1] * Gamma[, i - 2L]
    }

  tmp <- sigs[-1]
  tmp[1] <- 1
  Gamma / rep(sqrt(tmp), each = nc)
}

y confirmar que esto da la matriz correcta

# from whuber's answer
set.seed(1)
lm_method <- function(d, n = d * 4){
  x <- rnorm(n, mean = 2)
  x_p <- outer(x, 1:d, `^`)
  colnames(x_p) <- paste0("x", 1:d)
  poly_obj <- poly(x, d)

  list(poly_obj = poly_obj, gamma = coef(lm(cbind(1, poly_obj) ~ x_p)))
}

# check that we get the same with different degrees
for(d in 1:10){
  dat <- lm_method(d)
  stopifnot(all.equal(
    dat$gamma, get_poly_orth_map(dat$poly_obj), check.attributes = FALSE))
}

¿Reconstruir coeficientes brutos (y obtener sus varianzas) a partir de coeficientes ajustados a un polinomio ortogonal...

1. imposible de hacer y estoy perdiendo el tiempo.

Como muestra Whuber, no lo es. Como otra adición, aquí hay un ejemplo para obtener los errores estándar de las estimaciones como se menciona en los comentarios

# from `help(cars)`
fm <- lm(dist ~ speed + I(speed^2) + I(speed^3), data = cars)
summary(fm)
#R> Call:
#R> lm(formula = dist ~ speed + I(speed^2) + I(speed^3), data = cars)
#R> 
#R> Residuals:
#R>    Min     1Q Median     3Q    Max 
#R> -26.67  -9.60  -2.23   7.08  44.69 
#R> 
#R> Coefficients:
#R>             Estimate Std. Error t value Pr(>|t|)
#R> (Intercept) -19.5050    28.4053   -0.69     0.50
#R> speed         6.8011     6.8011    1.00     0.32
#R> I(speed^2)   -0.3497     0.4999   -0.70     0.49
#R> I(speed^3)    0.0103     0.0113    0.91     0.37
#R> 
#R> Residual standard error: 15.2 on 46 degrees of freedom
#R> Multiple R-squared:  0.673, Adjusted R-squared:  0.652 
#R> F-statistic: 31.6 on 3 and 46 DF,  p-value: 3.07e-11

fp <- lm(dist ~ poly(speed, 3), data = cars)
gamma <- get_poly_orth_map(poly(cars$speed, 3))
drop(gamma %*% coef(fp))
#R> [1] -19.5050   6.8011  -0.3497   0.0103
sqrt(diag(tcrossprod(gamma %*% vcov(fp), gamma)))
#R> [1] 28.4053  6.8011  0.4999  0.0113