[R] different random effects for each level of a factor in lme

[José Rafael Ferrer Paris]

I have an interesting lme - problem. The data is part of the Master
Thesis of my friend, and she had some problems analysing this data,
until one of her Jurors proposed to use linear mixed-effect models. I'm
trying to help her since she has no experience with R. I'm very used to
R but have very few experience with lme.

The group calls of one species of parrot were recorded at many
localities in mainland and islands. Within the localities the parrots
move in groups, several calls were recorded for each group but the calls
can not be separated by individuals.

We use the variable s1 to measure one property of the calls (the length
of the first part of the call). We are interested in explaining the
variability of the calls, one hypothesis is that variability of calls
tends to be greater in islands compared with mainland.

So we have 
s1   : as a measure of interest
loc  : as a grouping variable (locality)
grp  : as a grouping variable (group) nested in loc
isla : is an factor that identify which localities are in island 
       and which are in mainland (it is outer to loc)

I began with a simple model with fixed effects in isla (since there are
some differences in the length of s1) and nested random effects: 

f00 <- lme(s1~isla,data=s1.ap,
           random=~1|loc/grp)

My final model should have fixed effects in isla, different nested
random for both levels of isla, and different error per stratum,
something like:

f11 <-
  lme(s1~isla,data=s1.ap,
      random=~isla|loc/grp,
      weights=varIdent(form=~1|isla))

or perhaps:

f11b <-
  lme(s1~isla,s1.ap,
      random=~isla-1|loc/grp,
      weights=varIdent(form=~1|isla))

Is this a valid formulation? I have seen that ~x|g1/g2 is usually used
for modelling random effects in (the intercept and slope of) covariates
and that ~1|g1/f1 is used to model interactions between grouping factors
and treatment factors. 

I fitted the above models (and a few other variants between the simpler
and the complex ones) and found f11 to be the best model, f11 and f11b
are both identical in terms of AIC or LR-Tests since they are the same
model with different parametrization (I guess...).

Now, I suppose I did everything right, and I want to compare the
variance decomposition in islands and mainland, I use 

> VarCorr(f11)

            Variance        StdDev   Corr  
loc =       pdLogChol(isla)                
(Intercept) 1643.5904       40.54122 (Intr)
islaT        962.2991       31.02095 -0.969
grp =       pdLogChol(isla)                
(Intercept)  501.7315       22.39936 (Intr)
islaT        622.5393       24.95074 -0.818
Residual     547.0888       23.38993       

> VarCorr(f11b)

         Variance            StdDev   Corr 
loc =    pdLogChol(isla - 1)               
islaI    1643.4821           40.53988 islaI
islaT     168.6514           12.98659 0    
grp =    pdLogChol(isla - 1)               
islaI     501.7357           22.39946 islaI
islaT     209.8698           14.48688 0    
Residual  547.0871           23.38989      

The variance for islands (islaI) is always greater  than the ones for
mainland (islaT) as expected, and the estimates of Intercept in f11 are
nearly equal to the estimates for islaI in f11b. However, the estimates
for islaT are completely different. It seems to me that the estimates in
f11b are the "correct ones" and the ones in f11 are obtained by
reparametrization of the variance-covariance matrix. Am I right?

I want to say what percentage of the variance is explained by each
level, something like this:

> tmp <-
data.frame(I=c(1643.48,501.74,(547.09)),T=c(168.65,209.86,(0.7823097*547.09)))
> rownames(tmp) <- c("loc","grp","res")
> t(t(tmp)/(colSums(tmp)))
          I       T
loc  0.6104349 0.2091125
grp  0.1863604 0.2602096
res  0.2032047 0.5306780

(0.7823097 was the result of varIdent for islaT)

If I compare the sum of variances in f11b for each level of isla with
the variances of the data frame I get similar results:

> colSums(tmp)
        I         T 
2692.3100  806.5038 
> tapply(s1.ap$s1,s1.ap$isla,var)
        I         T 
2417.1361  731.8165 

So I guess this is the right way to interpret the variances in the
fitted model. Or, is it not?

thanks,

JR

-- 
Dipl.-Biol. JR Ferrer Paris
~~~~~~~~~~~~~~~~~~~~~~~~~~~
Laboratorio de Biología de Organismos --- Centro de Ecología
Instituto Venezolano de Investigaciones Científicas (IVIC) 
Apdo. 21827, Caracas 1020-A 
República Bolivariana de Venezuela

Tel: (+58-212) 504-1452
Fax: (+58-212) 504-1088

email: jferrer at ivic.ve
clave-gpg: 2C260A95