Interprétation des résultats de spline

20

J'essaie d'adapter une spline pour un GLM à l'aide de R. Une fois que j'ai ajusté la spline, je veux pouvoir prendre mon modèle résultant et créer un fichier de modélisation dans un classeur Excel.

Par exemple, supposons que j'ai un ensemble de données où y est une fonction aléatoire de x et la pente change brusquement à un point spécifique (dans ce cas @ 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)

J'adapte maintenant ceci en utilisant

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

et mes résultats montrent

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

À ce stade, je peux utiliser la fonction de prédiction dans r et obtenir des réponses parfaitement acceptables. Le problème est que je veux utiliser les résultats du modèle pour créer un classeur dans Excel.

Ma compréhension de la fonction prédire est que, étant donné une nouvelle valeur "x", r branche ce nouveau x dans la fonction spline appropriée (soit la fonction pour les valeurs supérieures à 500 ou celle pour les valeurs inférieures à 500), alors il prend ce résultat et multiplie il par le coefficient approprié et à partir de ce point le traite comme tout autre terme de modèle. Comment obtenir ces fonctions splines?

(Remarque: je me rends compte qu'un GLM gamma lié au journal peut ne pas être approprié pour l'ensemble de données fourni. Je ne demande pas comment ni quand adapter les GLM. Je fournis cet ensemble comme exemple à des fins de reproductibilité.)

Eric
la source
7
Je suggère, si possible, d'éviter d'inclure du code qui supprime toutes les variables ( rm(list=ls())), surtout pas sans avertissement. Quelqu'un peut copier-coller votre code dans une session ouverte de R où ils ont des variables déjà (mais pas appelés x, y, dfou spline1) et manque que votre code efface leur travail. Est-ce un peu stupide pour eux de faire ça? Oui. Mais il est toujours poli de les laisser décider quand supprimer leurs propres variables.
Glen_b -Reinstate Monica

Réponses:

25

Vous pouvez effectuer une rétro-ingénierie des formules splines sans avoir à entrer dans le Rcode. Il suffit de savoir que

  • Une spline est une fonction polynomiale par morceaux.

  • dd+1

  • Les coefficients d'un polynôme peuvent être obtenus par régression linéaire.

d+1xxdd=34×4=16d+1=4x

64RR

Cette méthode fonctionnera avec tous les logiciels statistiques, même les logiciels propriétaires non documentés dont le code source n'est pas disponible.

200,500,800(1,1000)RR

Tracés R

Tracés Excel

(Les quadrillages gris verticaux dans la Rversion indiquent où se trouvent les nœuds internes.)


Voici le Rcode complet . C'est un hack non sophistiqué, reposant entièrement sur la pastefonction pour effectuer la manipulation de chaîne. (Une meilleure façon serait de créer un modèle de formule et de le remplir à l'aide de commandes de correspondance et de substitution de chaînes.)

#
# 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 première formule de sortie spline (sur les quatre produites ici) est

"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)"

Rxx

Extrait Excel

whuber
la source
2
ns.formula.. pensez- vous en R?! Sérieusement, votre méthode semble très utile, mais il semble ironique de devoir pirater un hack pour obtenir ces paramètres. Serait très utile pour sortir un tableau ..
geotheory
Cela peut être une question stupide: mais s'agit-il de 4 splines que vous tracez ou de 4 bases d'une spline?
Erosennin
@Erosennin I dépend de ce que vous entendez par "une spline". Ces quatre courbes sont à la base d'une spline qui est cubique par morceaux en quatre intervalles et continuellement seconde différenciable aux trois points où ces intervalles se rencontrent, comme décrit par les trois puces qui introduisent ma réponse.
whuber
Merci! Je ne voulais pas être tatillonne, il semble qu'il y ait quatre splines (de la réponse), et non quatre courbes qui sont une base. Encore une fois, j'essaie juste de comprendre ...
Erosennin
1
@Erosennin Pas de problème. Cela peut peut-être aider: la "spline" est la combinaison linéaire de ces quatre courbes qui est déterminée par le processus d'ajustement de régression. Une autre façon de le dire: la spline est constituée d'un espace vectoriel de courbes qui peut être créé en prenant des combinaisons linéaires de ces quatre courbes.
whuber
4

Vous avez déjà fait ce qui suit:

> 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"))
> 

Maintenant, je vais vous montrer comment prédire (la réponse) pour x = 12 de deux manières différentes: D'abord en utilisant la fonction de prédiction (le plus simple!)

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

La 2ème voie est basée directement sur la matrice du modèle. Remarque J'ai utilisé expcar la fonction de liaison utilisée est log.

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

Notez que ci-dessus j'ai extrait le 12ème élément, car cela correspond à x = 12. Si vous souhaitez prédire un x en dehors de l'ensemble d'entraînement, vous pouvez simplement utiliser à nouveau la fonction de prédiction. Disons que nous voulons trouver la valeur de réponse prédite pour x = 1100 puis

> predict(spline1, newdata=data.frame(x=1100),type="response")
       1 
366.3483 
Stat
la source
Merci pour votre réponse! Mais, je suis toujours confus: /. Je ne suis pas sûr de savoir quoi faire avec cette matrice. Par exemple, si j'avais x = 12, alors prédire dit y = 68,78721, mais en levant 12 de cette matrice, j'obtiens 0,016816392. L'ordonnée à l'origine et le coefficient pour x <500 sont 4,174603 et 3,830416, respectivement. exp (4,174603 + 3,8304116 * 0,016816392) <> 68,78721. De plus, comment pourrais-je obtenir des valeurs pour x si x n'était pas dans l'ensemble d'apprentissage?
Eric
J'ai changé ma réponse.
Stat
J'ai ajouté un code pour le cas où x n'était pas dans l'ensemble de formation.
Stat
2
Existe-t-il un moyen d'obtenir 366,3483 pour x = 1100 sans utiliser la fonction de prédiction?
Eric
4

Vous pouvez trouver plus facile d'utiliser la base de puissance tronquée pour les splines de régression cubique, en utilisant le rmspackage R. Une fois le modèle ajusté, vous pouvez récupérer la représentation algébrique de la fonction spline ajustée à l'aide des fonctions Functionou .latexrms

Frank Harrell
la source
Je vous remercie. J'ai lu votre réponse ici stats.stackexchange.com/questions/67607/… avant de poster. Je suppose que j'ai juste besoin de mieux comprendre ce que je peux faire avec la RMS.
Eric
La documentation de Function()ne dit pas vraiment ce qu'elle fait. Dans mon cas (voir les détails sur Rpubs rpubs.com/EmilOWK/rms_splines ), j'obtiens function(x = NA) {-2863.7787+245.72672* x-0.1391794*pmax(x-10.9,0)^3+0.27835881*pmax(x-50.5,0)^3-0.1391794*pmax(x-90.1,0)^3 } <environment: 0x556156e80db8>La -2863.7787valeur est le premier coef du modèle, 245.72672le second et le dernier coef -873.0223n'est pas vu dans l'équation nulle part. Il en va de même pour la sortie de latex().
Deleet
Functionfonctionne avec Glm()lorsque vous utilisez rcsla fonction spline. La sortie reformule la spline sous sa forme la plus simple en écrivant comme si les restrictions de queue linéaires n'étaient pas là (mais elles le sont) comme détaillé dans mes notes de cours RMS .
Frank Harrell