Interpretación de resultados de spline

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

Por ejemplo, supongamos que tengo un conjunto de datos donde y es una función aleatoria de xy 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.

Entiendo que la función de predicción es que, dado un nuevo valor de "x", r conecta 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 se multiplica por el coeficiente apropiado y desde ese punto lo trata como cualquier otro término modelo. ¿Cómo obtengo estas funciones de spline?

(Nota: me doy cuenta de que un GLM gamma vinculado a un registro puede no ser apropiado para el conjunto de datos proporcionado. No estoy preguntando cómo o cuándo ajustar los GLM. Estoy proporcionando ese conjunto como un ejemplo para fines de reproducibilidad).

Puede aplicar ingeniería inversa a las fórmulas de spline sin tener que ir al Rcódigo. Es suficiente saber que

• Una spline es una función polinómica por partes.

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

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

$d+1$$x$$x^d$$d=3$$4\times 4=16$$d+1=4$$x$

$64$RR

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

$200, 500, 800$$(1, 1000)$RR

(Las líneas grises verticales en la Rversión muestran dónde están los nudos internos).

Aquí está el Rcódigo completo . Es un truco poco sofisticado, que se basa completamente en la pastefunción para lograr la manipulación de la cadena. (Una mejor manera sería crear una plantilla de fórmula y completarla utilizando 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 de spline (de las cuatro producidas 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)"

R$x$$x$

Ya hiciste 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 le mostraré cómo predecir (la respuesta) para x = 12 de dos maneras diferentes: Primero, use 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 segunda forma se basa directamente en la matriz del modelo. Nota que utilicé expya que la función de enlace utilizada es log.

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

Tenga en cuenta que en el anterior extraje el elemento 12, ya que corresponde a x = 12. Si desea predecir una x fuera del conjunto de entrenamiento, simplemente puede volver a utilizar la función de predicción. Digamos que queremos encontrar el valor de respuesta pronosticado para x = 1100 y luego

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

Puede resultarle más fácil utilizar la base de potencia truncada para splines de regresión cúbica, utilizando el rmspaquete R. Una vez que ajuste el modelo, puede recuperar la representación algebraica de la función de spline ajustada utilizando las funciones Functiono latexen rms.