Comment prédire de ns paramètres spline sans objet modèle

Question

j'ai les coefficients d'un GLM monté dans R, et je veux prédire les valeurs attendues pour une nouvelle série de données. Si j'avais l'objet model, ce serait simple, en utilisant predict(). Cependant, je suis maintenant hors site et pour des raisons de confidentialité des données je n'ai plus l'objet modèle. Je n'ai que l'objet de résumé, généré en utilisant résumé (modèle), qui contient les coefficients du modèle.

Il est assez facile d'utiliser les coefficients pour prédire les valeurs attendues pour un modèle simple. Cependant, je voudrais savoir comment faire lorsque le modèle comprend une ns splines cubiques(). Tous les raccourcis pour quand le modèle comprend également des variables catégoriques seraient également appréciés.

Voici un exemple simple.

library(splines) 
dat <- data.frame(x=1:500, z=runif(500), k=as.factor(sample(c("a","b"), size=500, replace=TRUE))) 
kvals <- data.frame(kn=c("a","b"),kv=c(20,30)) 
dat$y = dat$x + (40*dat$z)^2 + kvals$kv[match(dat$k,kvals$kn)] + rnorm(500,0,30) 
# Fit model 
library(splines) 
mod <- glm(y ~ x + ns(z,df=2) + k,data=dat) 
# Create new dataset 
dat.new <- expand.grid(x=1:3,z=seq(0.2,0.4,0.1),k="b") 
# Predict expected values in the usual way 
predict(mod,newdata=dat.new) 
summ <- summary(mod) 
rm(mod) 
# Now, how do I predict using just the summary object and dat.new? 

Answer 1

Il y a probablement une méthode plus efficace pour lutter contre cela, mais voici un point de départ pour vous aider à mettre en place pour mettre en œuvre la stratégie brièvement Roland suggéré. L'objet summ dispose des informations nécessaires pour définir la fonction spline, mais il est un peu enterré:

names(summ) 
    [1] "call"   "terms"   "family"   "deviance"  "aic"   
    [6] "contrasts"  "df.residual" "null.deviance" "df.null"  "iter"   
    [11] "deviance.resid" "coefficients" "aliased"  "dispersion"  "df"    
    [16] "cov.unscaled" "cov.scaled"  

Et regardant la structure de la feuille terms, nous voyons que le détail spline est enterré plus profondément encore dans la predvars subleaf :

str(summ$terms) 
Classes 'terms', 'formula' language y ~ x + ns(z, df = 2) + k 
    ..- attr(*, "variables")= language list(y, x, ns(z, df = 2), k) 
    ..- attr(*, "factors")= int [1:4, 1:3] 0 1 0 0 0 0 1 0 0 0 ... 
    .. ..- attr(*, "dimnames")=List of 2 
    .. .. ..$ : chr [1:4] "y" "x" "ns(z, df = 2)" "k" 
    .. .. ..$ : chr [1:3] "x" "ns(z, df = 2)" "k" 
    ..- attr(*, "term.labels")= chr [1:3] "x" "ns(z, df = 2)" "k" 
    ..- attr(*, "order")= int [1:3] 1 1 1 
    ..- attr(*, "intercept")= int 1 
    ..- attr(*, "response")= int 1 
    ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
    ..- attr(*, "predvars")= language list(y, x, ns(z, knots = structure(0.514993450604379, .Names = "50%"), Boundary.knots = c(0.00118412892334163, 0.99828373757191), intercept = FALSE), k) 
    ..- attr(*, "dataClasses")= Named chr [1:4] "numeric" "numeric" "nmatrix.2" "factor" 
    .. ..- attr(*, "names")= chr [1:4] "y" "x" "ns(z, df = 2)" "k" 

tirer donc l'attribut out:

str(attributes(summ$terms)$predvars) 
language list(y, x, ns(z, knots = structure(0.514993450604379, .Names = "50%"), 
       Boundary.knots = c(0.00118412892334163, 0.99828373757191), intercept = FALSE), k) 

Vous pouvez voir qu'il est po ssibles pour récupérer la spline si vous fournissez les x, y, z et k valeurs nécessaires:

with(dat, ns(z, knots = 0.514993450604379, Boundary.knots = c(0.00118412892334163, 
0.99828373757191), intercept = FALSE)) 
#--- 
       1    2 
    [1,] 5.760419e-01 -1.752762e-01 
    [2,] 2.467001e-01 -1.598936e-01 
    [3,] 4.392684e-01 4.799757e-01 
snipping .... 
[498,] 4.965628e-01 -2.576437e-01 
[499,] 5.627389e-01 1.738909e-02 
[500,] 2.393920e-02 -1.611872e-02 
attr(,"degree") 
[1] 3 
attr(,"knots") 
[1] 0.5149935 
attr(,"Boundary.knots") 
[1] 0.001184129 0.998283738 
attr(,"intercept") 
[1] FALSE 
attr(,"class") 
[1] "ns"  "basis" "matrix" 

Vous pouvez construire un dat substitut, si vous connaissez les extrêmes de vos données. Voir ?ns et les autres pages d'aide vers lesquels il pointe.