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/
pdIdent in smoothing regression model
3 messages · Lei Liu, Gavin Simpson
7 days later
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 %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%