Interpretación de los resultados del spline

Question

Estoy tratando de ajustar un spline para un GLM usando R. Una vez que ajuste el spline, quiero ser capaz de tomar mi modelo resultante y crear un archivo de modelado en un libro de trabajo de Excel.

Por ejemplo, digamos que tengo un conjunto de datos donde y es una función aleatoria de x y la pendiente cambia abruptamente en un punto específico (en este caso @ x=500).

set.seed(1066)
x<- 1:1000
y<- rep(0,1000)

y[1:500]<- pmax(x[1:500]+(runif(500)-.5)*67*500/pmax(x[1:500],100),0.01)
y[501:1000]<-500+x[501:1000]^1.05*(runif(500)-.5)/7.5

df<-as.data.frame(cbind(x,y))

plot(df)

Ahora encajo esto usando

library(splines)
spline1 <- glm(y~ns(x,knots=c(500)),data=df,family=Gamma(link="log"))

y mis resultados muestran

summary(spline1)

Call:
glm(formula = y ~ ns(x, knots = c(500)), family = Gamma(link = "log"), 
    data = df)

Deviance Residuals: 
     Min       1Q   Median       3Q      Max  
-4.0849  -0.1124  -0.0111   0.0988   1.1346  

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)             4.17460    0.02994  139.43   <2e-16 ***
ns(x, knots = c(500))1  3.83042    0.06700   57.17   <2e-16 ***
ns(x, knots = c(500))2  0.71388    0.03644   19.59   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Gamma family taken to be 0.1108924)

    Null deviance: 916.12  on 999  degrees of freedom
Residual deviance: 621.29  on 997  degrees of freedom
AIC: 13423

Number of Fisher Scoring iterations: 9

En este punto, puedo usar la función de predicción dentro de r y obtener respuestas perfectamente aceptables. El problema es que quiero usar los resultados del modelo para construir un libro de trabajo en Excel.

Mi comprensión de la función de predicción es que dado un nuevo valor "x", r enchufa esa nueva x en la función de spline apropiada (ya sea la función para valores superiores a 500 o la de valores inferiores a 500), luego toma ese resultado y lo multiplica por el coeficiente apropiado y a partir de ese punto lo trata como cualquier otro término modelo. ¿Cómo obtengo estas funciones spline?

(Nota: Me doy cuenta de que un GLM gamma ligado al registro puede no ser apropiado para el conjunto de datos proporcionado. No estoy preguntando cómo o cuándo encajar los GLM. Estoy proporcionando ese conjunto como un ejemplo para propósitos de reproducibilidad).

Answer 1

Ya has hecho lo siguiente:

> rm(list=ls())
> set.seed(1066)
> x<- 1:1000
> y<- rep(0,1000)
> y[1:500]<- pmax(x[1:500]+(runif(500)-.5)*67*500/pmax(x[1:500],100),0.01)
> y[501:1000]<-500+x[501:1000]^1.05*(runif(500)-.5)/7.5
> df<-as.data.frame(cbind(x,y))
> library(splines)
> spline1 <- glm(y~ns(x,knots=c(500)),data=df,family=Gamma(link="log"))
> 

Ahora les mostraré cómo predecir (la respuesta) para x=12 de dos maneras diferentes: Primero usando la función de predicción (¡la manera fácil!)

> new.dat=data.frame(x=12)
> predict(spline1,new.dat,type="response")
       1 
68.78721 

La 2ª vía se basa directamente en la matriz modelo. Nota: He utilizado exp ya que la función de enlace utilizada es la de registro.

> m=model.matrix( ~ ns(df$x,knots=c(500))) 
> prd=exp(coefficients(spline1) %*% t(m)) 
> prd[12]
[1] 68.78721

Obsérvese que en la parte superior extraje el 12º elemento, ya que corresponde a x=12. Si quieres predecir para una x fuera del conjunto de entrenamiento, entonces simplemente puedes usar de nuevo la función de predicción. Digamos que queremos encontrar el valor de respuesta pronosticado para x=1100 entonces

> predict(spline1, newdata=data.frame(x=1100),type="response")
       1 
366.3483

Answer 2

Es posible que le resulte más fácil utilizar la base de potencia truncada para las ranuras de regresión cúbica, utilizando la R rms paquete. Una vez que se ajusta el modelo, se puede recuperar la representación algebraica de la función de spline ajustada usando el Function o latex funciona en rms .

Answer 3

Se puede hacer ingeniería inversa de las fórmulas de la lengüeta sin tener que entrar en el R código. Basta con saber que

• Un spline es una función polinomial a destajo.

• Polinomios de grado $d$ están determinados por sus valores en $d+1$ puntos.

• Los coeficientes de un polinomio se pueden obtener mediante una regresión lineal.

Por lo tanto, sólo tienes que crear $d+1$ puntos espaciados entre cada par de nudos sucesivos (incluyendo los puntos finales implícitos del rango de datos), predecir los valores de spline y retroceder la predicción contra las potencias de $x$ hasta $x^d$ . Habrá una fórmula separada para cada elemento de la base de la ranura dentro de cada uno de esos nudos "bin". Por ejemplo, en el ejemplo siguiente hay tres nudos internos (para cubos de cuatro nudos) y splines cúbicos ( $d=3$ ), lo que dio lugar a $4 \times 4=16$ polinomios cúbicos, cada uno con $d+1=4$ coeficientes. Debido a que los poderes relativamente altos de $x$ es imperativo preservar toda la precisión de los coeficientes. Como se puede imaginar, la fórmula completa de cualquier elemento de base de spline puede ser bastante larga!

Como mencioné hace bastante tiempo El hecho de poder utilizar la salida de un programa como la entrada de otro (sin intervención manual, que puede introducir errores irreproducibles) es una útil habilidad de comunicación estadística. Esta pregunta proporciona un buen ejemplo de cómo se aplica ese principio: en lugar de copiar los $64$ coeficientes de dieciséis dígitos manualmente, podemos hackear juntos una forma de convertir los splines computados por R en fórmulas que Excel puede entender. Todo lo que necesitamos hacer es extraer los coeficientes de spline de R como se describe arriba, que las reformatee en fórmulas similares a las de Excel, y que las copie y pegue en Excel.

Este método funcionará con cualquier software estadístico, incluso con software propietario no documentado cuyo código fuente no esté disponible.

He aquí un ejemplo tomado de la pregunta, pero modificado para tener nudos en tres puntos internos ( $200, 500, 800$ ) así como en los puntos finales $(1, 1000)$ . Las tramas muestran R seguido de la versión de Excel. Se realizó muy poca personalización en ambos entornos (aparte de especificar los colores en R para que coincidan aproximadamente con los colores por defecto de Excel).

(Las cuadrículas grises verticales en el R muestran dónde están los nudos internos).

---

Aquí está el completo R código. Es un hacking poco sofisticado, que se basa completamente en el paste para llevar a cabo la manipulación de las cuerdas. (Una mejor manera sería crear una plantilla de fórmula y rellenarla usando los comandos de coincidencia y sustitución de cadenas).

#
# Create and display a spline basis.
#
x <- 1:1000
n <- ns(x, knots=c(200, 500, 800))

colors <- c("Orange", "Gray", "tomato2", "deepskyblue3")
plot(range(x), range(n), type="n", main="R Version",
     xlab="x", ylab="Spline value")
for (k in attr(n, "knots")) abline(v=k, col="Gray", lty=2)
for (j in 1:ncol(n)) {
  lines(x, n[,j], col=colors[j], lwd=2)
}
#
# Export this basis in Excel-readable format.
#
ns.formula <- function(n, ref="A1") {
  ref.p <- paste("I(", ref, sep="")
  knots <- sort(c(attr(n, "Boundary.knots"), attr(n, "knots")))
  d <- attr(n, "degree")
  f <- sapply(2:length(knots), function(i) {
    s.pre <- paste("IF(AND(", knots[i-1], "<=", ref, ", ", ref, "<", knots[i], "), ", 
                   sep="")
    x <- seq(knots[i-1], knots[i], length.out=d+1)
    y <- predict(n, x)
    apply(y, 2, function(z) {
      s.f <- paste("z ~ x+", paste("I(x", 2:d, sep="^", collapse=")+"), ")", sep="")
      f <- as.formula(s.f)
      b.hat <- coef(lm(f))
      s <- paste(c(b.hat[1], 
            sapply(1:d, function(j) paste(b.hat[j+1], "*", ref, "^", j, sep=""))), 
            collapse=" + ")
      paste(s.pre, s, ", 0)", sep="")
    })
  })
  apply(f, 1, function(s) paste(s, collapse=" + "))
}
ns.formula(n) # Each line of this output is one basis formula: paste into Excel

La primera fórmula de salida del spline (de las cuatro que se producen aquí) es

"IF(AND(1<=A1, A1<200), -1.26037447288906e-08 + 3.78112341937071e-08*A1^1 + -3.78112341940948e-08*A1^2 + 1.26037447313669e-08*A1^3, 0) + IF(AND(200<=A1, A1<500), 0.278894459758071 + -0.00418337927419299*A1^1 + 2.08792741929417e-05*A1^2 + -2.22580643138594e-08*A1^3, 0) + IF(AND(500<=A1, A1<800), -5.28222778473101 + 0.0291833541927414*A1^1 + -4.58541927409268e-05*A1^2 + 2.22309136420529e-08*A1^3, 0) + IF(AND(800<=A1, A1<1000), 12.500000000002 + -0.0375000000000067*A1^1 + 3.75000000000076e-05*A1^2 + -1.25000000000028e-08*A1^3, 0)"

Para que esto funcione en Excel, sólo hay que quitar las comillas que lo rodean y ponerle el prefijo "=". (Con un poco más de esfuerzo podrías tener R escribir un archivo que, al ser importado por Excel, contiene copias de estas fórmulas en todos los lugares correctos). Pégalo en una caja de fórmula y luego arrastra esa celda hasta que "A1" haga referencia a la primera $x$ donde el spline debe ser computado. Copia y pega (o arrastra y suelta) esa celda para calcular los valores de otras celdas. Llené las celdas B2:E:102 con estas fórmulas, haciendo referencia $x$ en las celdas A2:A102.