R help - Sep 2011 - Question about Natural Splines (ns function)

Hi - How can I 'manually' reproduce the results in 'pred1'
below? 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])

Regards,
Axel.

	[[alternative HTML version deleted]]

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