José Rafael Ferrer Paris
2007-Mar-06 03:12 UTC
[R] different random effects for each level of a factor in lme
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