[R-sig-ME] Specifying random effects in lme()
Gavin Simpson
gavin.simpson at ucl.ac.uk
Wed May 9 19:40:10 CEST 2007
On Wed, 2007-05-09 at 16:58 +0200, Nilsson Fredrik X wrote:
> Dear Gavin,
>
> It SHOULD work if you add
> random = ~log.so4+log.cl|site
Hi Nilsson,
Many thanks for this. I had actually tried this but when I tried coef()
on the resulting model, I got what looked like only random effects for
the log.cl parameter. Now that I have looked at the resulting object
more closely after your email, I have noticed that the magnitudes of the
random effects for log.so4 are very small indeed, and so when coef()
prints the values for log.so4, they do all look constant, but only
because they are printed to finite decimal places.
I should look more closely next time!
Many thanks, I now have some indication that the effect of log.cl does
vary between sites, but not log.so4.
All the best,
G
>
> I made a script which is able to do so, but perhaps the addition of
> this extra variability in your case over-specifies the model?
>
> Hope this helps. Otherwise I suppose that you have to supply more details.
>
> Cheers,
>
> Fredrik
> --------------------------------------------------------------------------
> Example:
>
> nsite<-10
> nwithin<-10
> prosw<-nsite*nwithin
>
> conc<-vector(length=prosw)
> temp<-rnorm(prosw)
> log.so4<-seq(-5,5, length=nwithin)
> log.cl<-rnorm(prosw)
> int<-17
> acl<-sqrt(2)
> aso4<-21
> atemp<- -14
> rho<--0.9
> for (i in 1:nsite)
> {
> sf<-0.1*rnorm(1)
> sc<-2*rnorm(1)
> ss<-rnorm(1)
> err<-0.1*rnorm(1)
> for (j in 1:nwithin)
> {
> conc[nwithin*(i-1)+j]<-sf + int + (acl+sc)*log.cl[nwithin*(i-1)+j] + (aso4+ss)*log.so4[j] +
> atemp*temp[nwithin*(i-1)+j] + err
> err<-err*rho + 0.1*rnorm(1)
> }
> }
>
> temp<-rep(temp,nsite)
> log.so4<-rep(log.so4, nsite)
> log.cl<-rep(log.cl,nsite)
> site<-rep(seq(1,nsite), each=nwithin)
>
> DF<-data.frame(conc,temp, log.so4, log.cl, site)
> DF$site<-as.factor(site)
>
> DFgd<-groupedData(conc~temp|site, data=DF)
> plot(DFgd)
>
> DF.lme<-lme(conc~log.so4+temp+log.cl, random=~log.so4+log.cl|site, data=DF)
> plot(ACF(DF.lme, resType="n", form=~1|site, maxLag=floor(nwithin/2-.5)), alpha=0.05)
> DF2.lme<-update(DF.lme, correlation=corAR1(form=~1|site))
>
>
>
> -----Ursprungligt meddelande-----
> Från: r-sig-mixed-models-bounces at r-project.org [mailto:r-sig-mixed-models-bounces at r-project.org] För Gavin Simpson
> Skickat: den 9 maj 2007 15:24
> Till: R-SIG-Mixed-Model
> Ämne: [R-sig-ME] Specifying random effects in lme()
>
> Dear List,
>
> I am modelling repeated measures data from 11 sites in the UK using
> lme() from the nlme package. I have the following call for a random
> intercept model with an AR1 error structure:
>
> mod <- lme(doc ~ temp + log.so4 + log.cl, random = ~ 1 | site,
> correlation = corAR1(), data = hydro)
>
> I would like to fit a model that has random effects for log.so4 and
> log.cl within site, so that I am fitting a random slop and intercept
> model --- we want the effect of log.so4 and log.cl to vary between sites
> and to compare the fits of the two models.
>
> If I am following lmer() in lme4 correctly, I believe I could fit this
> model as:
>
> hydro.lmer <- lmer(doc ~ temp + log.so4 + log.cl +
> (log.so4 | site) + (log.cl | site),
> data = hydro)
>
> but without the AR1 correlation.
>
> I am struggling to get the syntax correct for random in lme() to add
> random effects to log.so4 and log.cl. I have only managed to get a model
> where one of the coefficients (log.cl) also has a random effect.
>
> I would be most grateful if someone could explain how to fit the model
> in lme() by showing me the correct form for the random argument.
>
> Many thanks,
>
> Gav
>
--
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
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
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
More information about the R-sig-mixed-models
mailing list