Hi all, Im hoping someone might be able to help me out. Forgive me if my mistake is something simple. I am new to mixed models, new to R, and new to lme4 and am struggling to figure everything out. I have two questions that I am hoping someone can answer. 1) Am I using the correct random structure for my model? 2) Can someone help me figure out what is wrong with my syntax to code for random effect variance by treatment group?
These questions somewhat go together but let me tackle number one first. I was told by our statistician (who unfortunately doesnt know R well) that my model should include random effects for ID, ID*state, and ID*season. If my understanding of lme4 code is correct, my random structure would appear as it does in this model M1<-lmer (y~trt*season*state*site+(1|ID)+(1|state)+(1|season), data=d) Is this correct? (I am starting with the fully crossed fixed effects and I will use iterative model selection to find the optimal model after I make sure that I have the correct random terms) The second question is a bit more complicated and changes my random structure as well. I had originally built my model in nlme and used the multiple variance (varIdent) function to allow different variance for two of my terms (trt and state) in my nlme model because different levels in each term had different variances. There are two levels in each of these factors. Treatment is G or S (basically treated or control groups) and state is either Pre-treatment or Post-treatment. What is happening in my model is that the variance for the treated (G) is much smaller post-treatment than pre-treatment. Thus, overall, the variance of G is less than S and the variance of Post is less than Pre. My other fixed factors are a site factor (two locations) and a season factor (breeding or non-breeding season). These ones have no issues. I found a very helpful post ( https://stat.ethz.ch/pipermail/r-sig-mixed-models/2007q3/000248.html) on how to also do this in lmer but it doesnt seem to be working correctly for me and I have not been able to figure out why. I am also trying to do this with 2 factors instead of one. The overall structure of my data is Structure of my dataset (d) ..@ frame :'data.frame': 173 obs. of 10 variables: .. ..$ y : int [1:173] 209 382 448 353 224 112 198 273 495 622 ... .. ..$ trt : Factor w/ 2 levels "G","S": 1 1 1 1 1 1 1 1 1 1 ... .. ..$ season: Factor w/ 2 levels "Breeding","NonBreeding": 2 2 2 2 2 2 2 1 1 1 ... .. ..$ state : Factor w/ 2 levels "Post","Pre": 1 1 1 1 1 1 1 1 1 1 ... .. ..$ site : Factor w/ 2 levels "Mrak","Orchard": 1 1 1 1 1 2 2 1 1 1 ... .. ..$ G : num [1:173] 1 1 1 1 1 1 1 1 1 1 ... .. ..$ ID : Factor w/ 59 levels "100","103","105",..: 11 19 22 58 6 36 50 11 20 24 ... .. ..$ S : num [1:173] 0 0 0 0 0 0 0 0 0 0 ... .. ..$ Pre : num [1:173] 0 0 0 0 0 0 0 0 0 0 ... .. ..$ Post : num [1:173] 1 1 1 1 1 1 1 1 1 1 ... Following the recommendations from the post above, here is what I did to try to obtain the different variances per treatment group. d$G <- as.numeric(d$trt == "G") d$S <- as.numeric(d$trt == "S") d$Pre <- as.numeric(d$state == "Pre") d$Post <- as.numeric(d$state == "Post") m1<-lmer(y~trt*season*state*site+(0+G|ID)+(0+S|ID) +(0+Pre|ID) + (0+Post|ID)+(1|season), data=d) I then went to check what the variances for the random effects were and it is unfortunately showing all of the variance for trt in the S group and all of the variance for state in the Pre group with zero variance allocated to the other groups. Linear mixed model fit by REML Formula: y ~ trt * season * state * site + (0 + G | ID) + (0 + S | ID) + (0 + Pre | ID) + (0 + Post | ID) + (1 | season) Data: d AIC BIC logLik deviance REMLdev 2257 2327 -1107 2387 2213 Random effects: Groups Name Variance Std.Dev. ID Post 0.0 0.000 ID Pre 30688.9 175.183 ID S 30494.6 174.627 ID G 0.0 0.000 season (Intercept) 1436.8 37.905 Residual 46594.9 215.859 Number of obs: 173, groups: ID, 59; season, 2 Im not sure what I am doing wrong or even if I am remotely on the right page here but I am out of ideas. Can anyone help? Sara [[alternative HTML version deleted]]
______________________________________________ [email protected] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
