[R-sig-ME] Simple: don't understand lmer formula and slopes

Sat Feb 1 10:39:27 CET 2014

You are fitting a model that allows the slope of the relationship between ylong and xlong to differ among species. Let's plot the data in a way that more clearly illustrates what is going on:

library(nlme)
dat.group <- groupedData(ylong ~ xlong | species, data=trial.dat)
plot(dat.group, pch=19, cex=.5)

So, it looks like that slope is pretty much the same for the different species. Not a big surprise that the slope variance is estimated to be zero.

Best,
Wolfgang

--
Wolfgang Viechtbauer, Ph.D., Statistician
Department of Psychiatry and Psychology
School for Mental Health and Neuroscience
Faculty of Health, Medicine, and Life Sciences
Maastricht University, P.O. Box 616 (VIJV1)
6200 MD Maastricht, The Netherlands
+31 (43) 388-4170 | http://www.wvbauer.com
________________________________________
From: r-sig-mixed-models-bounces at r-project.org [r-sig-mixed-models-bounces at r-project.org] On Behalf Of Robert Chatfield [chatfield at alumni.rice.edu]
Sent: Saturday, February 01, 2014 1:51 AM
To: r-sig-mixed-models at r-project.org
Cc: Robert Chatfield
Subject: [R-sig-ME] Simple: don't understand lmer formula and slopes

Here's a very simple mis-understanding of formula specification in lme4 - lmer
(Thanks to Ben and this group for help with a more fundamental question),

How does one specify slopes in lmer?
(This illustrates the essentials of a more complex problem.  I know that
I could simply using  a rescaling  of the y variables with estimated
standard deviations)

(Again) I have a single control tracer variable, and a whole set
of dependent variables with varying slopes.  The way I am specifying
the problem, I get no random effect for a, the slope, as it should
depend on species.

The lmer formal

 test.1.lmer = lmer( ylong ~ xlong + (xlong - 1 | species) + 1 , data=trial.dat)

I provide an example below suggesting that I have built the
data frame correctly, I think.

# How this was run …
# Macintosh 10.8.5
#  R 2.15.2 GUI 1.53 Leopard build 64-bit (6335)
# lme4  version 0.999999-0

Bob Chatfield
Robert.B.Chatfield at nasa.gov
Earth Science Division, NASA (Ames), Mountain View, CA 94035 USA

 #  Test an lmer with a variable slope ... specification of random effect
> #
> set.seed(101) # for reprducibility
> ## Define a perturbation of the key tracer -- log normal
> x = 10 +  rlnorm( 100, log(15), 0.3)
> #
> ## Define some "emission factors" linearly relating several dependent variables
> ## as they are related to the key tracer (above background)
> ## 10 different dependent variables or "species"
> a = 0.1*rexp(10,rate=0.6)
> ## Define the pure response variable (y ... "y unscaled")
> y = array(dim=c(100,10))
> for ( i in 1:100 ) {
+   for ( j in 1:10) y[i,j] = a[j] * xd[i]
+ }
> ## Add some noise in both the response, both w.r.t. the dependent species
> ##    and also the measurement instance
> for (i in 1:10 ) y[,i] = y[,i]*(1+rnorm(100,0.15,0.1))
> for (j in 1:100 ) y[j,] = y[j,]*(1+rnorm(10,0.15,0.1))
> #
> matplot(x,y,xlim=c(0,45),ylim=c(-1.5,max(as.vector(y))),cex=0.5)
> # points(x,y[,9],cex=2)
> points(x0,rep(0.025,100),cex=0.8)
> #
> ## Make a "long" list of the tracer-species, repeating once for every response
> ## Note:  normally I use reshape() to do this.  See: ?reshape
> xlong = rep(x,times=10)
> ylong = as.vector(y) # collapse the matrix into a similar "long" vector
> ## Make a factor variable describing each measurement
> ## instance and each dependent variable
> species = as.factor(rep(paste("species",as.character(1:10),sep=""),times=100))
> #
> xlong = rep(x,times=10)
> ylong = as.vector(y) # collapse the matrix into a similar "long" vector
> ## Make a factor variable describing each measurement
> ## instance and each dependent variable
> species = as.factor(rep(paste("species",as.character(1:10),sep=""),times=100))
> #
> trial.dat = data.frame(xlong=xlong,ylong=ylong,species=species)
> #  See example of structure below
> #
> test.1.lmer = lmer( ylong ~ xlong + (xlong - 1 | species) + 1 , data=trial.dat)
> summary(test.1.lmer)
Linear mixed model fit by REML
Formula: ylong ~ xlong + (xlong - 1 | species) + 1
   Data: trial.dat
  AIC  BIC logLik deviance REMLdev
 3896 3916  -1944     3877    3888
Random effects:
 Groups   Name  Variance Std.Dev.
 species  xlong 0.000    0.0000
 Residual       2.833    1.6832
Number of obs: 1000, groups: species, 10

Fixed effects:
            Estimate Std. Error t value
(Intercept) -0.96101    0.32930  -2.918
xlong        0.09945    0.01279   7.775

Correlation of Fixed Effects:
      (Intr)
xlong -0.987
> #
> cat( "standard deviation of original idealized data and of random effects")
standard deviation of original idealized data and of random effects
> (sd(a))
[1] 0.0821322
> sd( (ranef(test.1.lmer)$species)[["xlong"]])
[1] 0
> #
> # Add estimated slope to graph
> abline(fixef(test.1.lmer),col="darkred")
>

#  Example of the dataset

> trial.dat[1:3,]
     xlong     ylong  species
1 23.60230 0.7164348 species1
2 27.70397 0.9228682 species2
3 22.25050 0.5650994 species3
      xlong     ylong   species
9  29.75010 0.8829359  species9
10 24.02824 0.6977277 species10

> trial.dat[11:13,]
      xlong     ylong  species
11 27.56634 0.6923857 species1
12 21.81768 0.5428702 species2
13 33.02031 0.9520825 species3
> trial.dat[101:103,]
       xlong     ylong  species
101 23.60230 0.3621467 species1
102 27.70397 0.4328762 species2
103 22.25050 0.3848661 species3
> trial.dat[111:113,]
       xlong     ylong  species
111 27.56634 0.5934288 species1
112 21.81768 0.3077124 species2
113 33.02031 0.8165563 species3
_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models