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@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.
