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

[Robert Chatfield]

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