Skip to content
Prev 20241 / 20628 Next
Original Raw

Random effects in R vs SAS

John Maindonald 
Will, an immediate comment is that what may be accommodated are
variance ?components'.  They are not the variances of any identifiable
quantity. It is surely misleading to call them variances.  For residuals,
the component is at the same time a variance.  With a suitably balanced
experimental design, one can work with the function aov() to obtain a
breakdown from which the negative component can be extracted,

The following inflates within block variances ? it mirrors a context
where blocks have been designed to be more heterogenous than
plots generally.
library(lme4)
   df0 <- data.frame(block=rep(1:100, rep(10,100)), trt=rep(1:5,200),
                  y=rnorm(1000,4,1))
  ## Plot variances are all set to 1.0
  df0$y <- df0$y + df0$trt*0.1  ## Differences between treatments
  ## Add systematic within block variation
  ## 10 different systematic effects for the 10 different plots in
  ## a block are randomly assigned  
  systeff <- 0.4*(1:10)
  df <- df0
  for(i in 1:100)df$y[df$block==i] <- df$y[df$block==i]+sample(systeff)
  ## sample does the random assignment
  df$block = factor(df$block); df$trt =factor(df$trt)
  y.aov <- aov(y ~ trt +Error(block), data=df)
Error: block
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals 99  127.7    1.29               

Error: Within
           Df Sum Sq Mean Sq F value Pr(>F)
trt         4   23.5   5.865   2.302 0.0569
Residuals 896 2282.8   2.548
boundary (singular) fit: see help('isSingular?)
Linear mixed model fit by REML ['lmerMod']
Formula: y ~ trt + (1 | block)
   Data: df

REML criterion at convergence: 3776.3

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.78060 -0.72441  0.04553  0.70139  2.79109 

Random effects:
 Groups   Name        Variance Std.Dev.
 block    (Intercept) 0.000    0.000   
 Residual             2.537    1.593   
Number of obs: 1000, groups:  block, 100

. . .
Linear mixed-effects model fit by REML
  Data: df 
       AIC      BIC    logLik
  3790.336 3824.655 -1888.168

Random effects:
 Formula: ~1 | block
         (Intercept) Residual
StdDev: 6.727467e-05 1.592659

Fixed effects:  y ~ trt 
               Value Std.Error  DF  t-value p-value
(Intercept) 6.096161 0.1126180 896 54.13131  0.0000
trt2        0.305667 0.1592659 896  1.91922  0.0553
trt3        0.418814 0.1592659 896  2.62965  0.0087
trt4        0.515445 0.1592659 896  3.23638  0.0013
trt5        0.561525 0.1592659 896  3.52571  0.0004

. . .

In principle, one should be able to set up a within blocks correlation
structure that accommodates the negative component of variance.
John Maindonald
Statistics Reseach Associates
Wellington NZ.

  
  


  


  


  
  

    

  

    
      
    
    
Thread (4 messages)

  

  

    
            
        Will Hopkins
        
        Random effects in R vs SAS
      
      Aug 17
    
    
      
            
        John Maindonald
        
        Random effects in R vs SAS
      
      Aug 18
    
    
      
            
        Mollie Brooks
        
        Random effects in R vs SAS
      
      Aug 24
    
    
      
            
        Mollie Brooks
        
        Random effects in R vs SAS
      
      Aug 24