# [R] Question about Natural Splines (ns function)

David Winsemius dwinsemius at comcast.net
Wed Sep 7 05:25:13 CEST 2011

```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) * mm[,1] +
>                   coefficients(glm1) * mm[,2] +
>                   coefficients(glm1) * mm[,3] +
>                   coefficients(glm1) * mm[,4] +
>                   coefficients(glm1) * 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

```