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/
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: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:> 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 %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
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 %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%