[R] pdIdent in smoothing regression model

[Lei Liu]

Hi there,

I am reading the 2004 paper "Smoothing with mixed model software" in 
Journal of Statistical Software, by Ngo and Wand. I tried to run 
their first example in Section 2.1 using R but I had some problems. 
Here is the code:

library(nlme)
fossil <- read.table("fossil.dat",header=T)
x <- fossil$age
y <- 100000*fossil$strontium.ratio
knots <- seq(94,121,length=25)
n <- length(x)
X <- cbind(rep(1,n),x)
Z <- outer(x,knots,"-")
Z <- Z*(Z>0)
fit <- lme(y~-1+X,random=pdIdent(~-1+Z))

When I ran the code

fit <- lme(y~-1+X,random=pdIdent(~-1+Z))

I got an error message:

Error in getGroups.data.frame(dataMix, groups) :   Invalid formula for groups

I was really puzzled. I asked Dr. Ngo and he said they did it in 
S-plus but not R. Does anyone knows how to do it in R? Thanks!

Lei Liu
Associate Professor
Division of Biostatistics
Department of Public Health Sciences
University of Virginia School of Medicine

http://people.virginia.edu/~ll9f/