Re: [R] lme, spline

From: Gavin Simpson <gavin.simpson_at_ucl.ac.uk>
Date: Tue, 15 Jun 2010 23:26:46 +0100

On Tue, 2010-06-15 at 06:28 -0700, Farhad Shokoohi wrote:
> Dear All,
> I revise my question about the problem I have.
> Take a look at the article
> http://www.jstatsoft.org/v09/i01
> and download the attached code.
> try to run one of the codes for example section 2.1 in R
> here is the code

That is an older (can't believe I'm saying that for something published in 2004!) paper and contains code for SAS and S-PLUS. It doesn't surprise me that there are problems running it in R.

IIRC these ideas (or some of them) are in the SemiPar package, including the fossil data. Try there.

Alternative, I've been fiddling with modelling these data with various functions in R and found gam in package mgcv very useful. If your interests is in fitting similar style models you could try

require(SemiPar)
require(mgcv)
mod <- gam(strontium.ratio ~ s(age), data = fossil,

           method = "REML")
plot(strontium.ratio ~ age, data = fossil) pdat <- with(fossil,

             data.frame(age = seq(min(age), max(age),
                        length = 100)))

p1 <- predict(mod, newdata = pdat)
lines(p1 ~ age, data = pdat, col = "red")

Although I'm not 100% clear on this, I think you can get the same (or similar) to your lme model (though using a different smoother type) by doing

mod2 <- gamm(strontium.ratio ~ s(age), data = fossil) mod2$lme
mod2$gam

as gamm uses lme to fit the additive model in much the same manner as the code you show.

If you specifically want to fit exactly this code or know more about pdIdent, I can't help there I'm afraid. Try the maintainer of the nlme package or better still, try the R-SIG-Mixed mailing list where there are lots of knowledgeable people who use lme and lmer and who may be able to help.

HTH G

>
> 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))
>
> there is an error which prevents running the lme function. The error is
> something like this
> Error in getGroups.data.frame(dataMix, groups) :
> Invalid formula for groups
>
> Let me know what's wrong?
> Best

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

______________________________________________
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.
Received on Tue 15 Jun 2010 - 22:31:47 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Tue 15 Jun 2010 - 22:50:30 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive