[R] Question about Natural Splines (ns function)

[David Winsemius]

On Sep 6, 2011, at 5:08 PM, Axel Urbiz wrote:

> Hi - How can I 'manually' reproduce the results in 'pred1' below?

Well, not by constructing the prediction function from a different  
data basis that the glm function gets.

> My attempt
> is pred_manual, but is not correct. Any help is much appreciated.
>
> library(splines)
> set.seed(12345)
> y <- rgamma(1000, shape =0.5)
> age <- rnorm(1000, 45, 10)
> glm1 <- glm(y ~ ns(age, 4), family=Gamma(link=log))
> dd <- data.frame(age = 16:80)
> mm <- model.matrix( ~ ns(dd$age, 4))
> pred1 <- predict(glm1, dd, type = "response")
> pred_manual <- exp(coefficients(glm1)[1] * mm[,1] +
>                   coefficients(glm1)[2] * mm[,2] +
>                   coefficients(glm1)[3] * mm[,3] +
>                   coefficients(glm1)[4] * mm[,4] +
>                   coefficients(glm1)[5] * mm[,5])

 > attr(glm1$terms, "predvars")
list(y, ns(age, knots = c(38.7407342480734, 44.9093960482465,
51.5913373894399), Boundary.knots = c(14.7723335249845,  
76.5536098692015), intercept = FALSE))

 > pred_man <- exp(coefficients(glm1) %*% t(mm))
 > str(pred_man)
  num [1, 1:65] 0.327 0.335 0.343 0.351 0.36 ...
  - attr(*, "dimnames")=List of 2
   ..$ : NULL
   ..$ : chr [1:65] "1" "2" "3" "4" ...

 > str(pred1)
  Named num [1:65] 0.336 0.343 0.351 0.359 0.366 ...
  - attr(*, "names")= chr [1:65] "1" "2" "3" "4" ...
 > plot(pred_man, pred1)

So a 4 knot natural spline derived from a random rnorm  against a  
random gamma series is somewhat similar to one derived from from a  
series of integers over the same range.

 > dd2 <- data.frame(age = age)
 > mm2 <- model.matrix( ~ ns(dd2$age, 4))
 > pred_man2 <- exp(coefficients(glm1) %*% t(mm2))
 > plot(age,pred_man2)
 > points(dd$age, pred_man)

Different bases, different predictions.

-- 

David Winsemius, MD
Heritage Laboratories
West Hartford, CT