[R] pdIdent in smoothing regression model

[Gavin Simpson]

On Sun, 2011-10-09 at 10:25 -0400, Lei Liu wrote:
> 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:

[whoops - fingers slipped and sent the message prematurely]

I'm not going to try to fix this nor work out how to do it via lme(),
but just wanted to note that you can use mgcv:::gam() for example to fit
something very similar.

require(SemiPar)
require(mgcv)

data(fossil)
mod <- gam(strontium.ratio ~ s(age), data = fossil,
           method = "REML")

That is treating the smoothness terms as random effects.

Also note that Matt Wand wrote the SemiPar package as support software
for the book Semiparametric Regression, where these example data are
taken from:

http://cran.r-project.org/web/packages/SemiPar/index.html

You could use the function contained therein.

The above presumes you weren't wanting to make the example work for some
pedagogic reason.

HTH

G

> 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/
> 
> ______________________________________________
> R-help at 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.

-- 
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
 Dr. Gavin Simpson             [t] +44 (0)20 7679 0522
 ECRC, UCL Geography,          [f] +44 (0)20 7679 0565
 Pearson Building,             [e] gavin.simpsonATNOSPAMucl.ac.uk
 Gower Street, London          [w] http://www.ucl.ac.uk/~ucfagls/
 UK. WC1E 6BT.                 [w] http://www.freshwaters.org.uk
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%