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
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
