Simple: don't understand lmer formula and slopes
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
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