[R-sig-ME] simulation of factor-specific random effect variance estimation, the order of assigned variance matters?

[Paul Debes]

Hi Chen,

I think you may only have to change the random effect specification from  
variance contrast to variance means:

random = list(groupID = pdDiag(~prog))
to
random = list(groupID = pdDiag(~0+prog))

It may also be good to change "groupID" from a continuous covariate to a  
factor.
This should give the expected standard deviations for the random-term  
groups in both simulations.
Just as a note, please be aware that you do not specify variance but  
standard deviation in your simulations (model output is also in standard  
deviation).

Best,
Paul


On Wed, 27 Jul 2016 23:48:01 +0300, Chen Chun <talischen at hotmail.com>  
wrote:

> Dear all,
>
> I have a nested data structure and I would like to estimate  
> factor-specific random effect variance. Inspired by the article:
>
> http://rstudio-pubs-static.s3.amazonaws.com/6B98_c0ae95A0AAA44A3Aa8a9b6ba9703fcf5.html<http://rstudio-pubs-static.s3.amazonaws.com/6298_c0ae951011144131a8a9b6ba9703fcf5.html>
>
> I conducted some simulations:
>
> ###############
> ##--generate data-----
> ngroup <- 25
> nrep<- 5
> ##----random effect of group---
> meanh <- 0
> sigmahA <- 1   ##--these are the two assigned variance of the random  
> factor group for each class (A, B)
> sigmahB <- 3
>
> ##----residual random error---
> meanR <- 0
> sigmaR <- 0.2
>
> ##---simulate smaple--------------------------------------------------
> set.seed(55)
> out_A   <- NA
> out_B   <- NA
> raneff_groupA <- rnorm(ngroup, meanh, sigmahA)
> raneff_groupB <- rnorm(ngroup, meanh, sigmahB)
> groupID      <- seq(1, ngroup)
> for (i in 1:ngroup) {
>     raneff_red_A <- rnorm(nrep, meanR, sigmaR)
>     temp_A       <- 0 + raneff_groupA[i] + raneff_red_A
>     raneff_red_B   <- rnorm(nrep, meanR, sigmaR)
>     temp_B   <- 1 + raneff_groupB[i] + raneff_red_B
>     out_A <- rbind(out_A, cbind(groupID[i], raneff_groupA[i],  
> raneff_red_A, temp_A))
>     out_B   <- rbind(out_B, cbind(groupID[i], raneff_groupB[i],  
> raneff_red_B, temp_B))
> }
> out_A <- out_A[-1,]
> out_B   <- out_B[-1,]
>
> dat <- data.frame(groupID=c(out_A[,1], out_B[,1]),  
> prog=c(rep("A",nrow(out_A)), rep("B",nrow(out_B))),  
> out=c(out_A[,"temp_A"], out_B[,"temp_B"]), raneff_haul=c(out_A[,2],  
> out_B[,2]), error=c(out_A[,3], out_B[,3]))
>
> library(lattice)
> bwplot(out~paste(prog,groupID), data=dat)
>
> ##--apply model-----
> fit2 <- lme(out ~ prog-1, random = list(groupID = pdDiag(~prog)), method  
> = "REML", data = dat)
> summary(fit2)
>
> This code gives the expected model output: as estimated variance of  
> groups for A and B are:  1.07957 2.835139, similar to assigned value
> ###############
> ###############
> Model output:
> Linear mixed-effects model fit by REML
> Data: dat
>        AIC      BIC    logLik
>   202.7263 220.2934 -96.36314
>
> Random effects:
> Formula: ~prog | groupID
> Structure: Diagonal
>         (Intercept)    progB  Residual
> StdDev:     1.07957 2.835139 0.1959905
>
> Fixed effects: out ~ prog - 1
>           Value Std.Error  DF   t-value p-value
> progA 0.0316535 0.2166244 224 0.1461215  0.8840
> progB 1.1446339 0.6069982 224 1.8857287  0.0606
> Correlation:
>       progA
> progB 0.355
>
> Standardized Within-Group Residuals:
>         Min          Q1         Med          Q3         Max
> -2.29868033 -0.63163428  0.03981236  0.64863494  2.63613133
>
> Number of Observations: 250
> Number of Groups: 25
> #########################
>
>
>
> However, it I dont make any other changes, but simply to switch the  
> order of the simulation variance:
>
> sigmahA <- 3
> sigmahB <- 1
>
> , generate the data and apply the model, the model gives very confusing  
> result:
> ###############
> ###############
> Model output:
> Linear mixed-effects model fit by REML
> Data: dat
>       AIC      BIC    logLik
>   261.763 279.3301 -125.8815
>
> Random effects:
> Formula: ~prog | groupID
> Structure: Diagonal
>         (Intercept)    progB  Residual
> StdDev:    3.187986 3.295383 0.1959864
>
> Fixed effects: out ~ prog - 1
>           Value Std.Error  DF   t-value p-value
> progA 0.1455012 0.6378382 224 0.2281162  0.8198
> progB 1.0590270 0.9171802 224 1.1546553  0.2495
> Correlation:
>       progA
> progB 0.695
>
> Standardized Within-Group Residuals:
>         Min          Q1         Med          Q3         Max
> -2.30847626 -0.63592816  0.03830767  0.64349268  2.63737406
>
> Number of Observations: 250
> Number of Groups: 25
> #########################
>
>
>
> The model output gives equal variance. The same situation happens if I  
> run the same code in the article (link above) and also simply change the  
> order of the variance assigned to the groups.
>
> Can someone help me to figure out why the order of assigning the  
> variance matters? And why the second case does not gives the expected  
> results?
>
> Thanks
>
> Regards,
> Chun Chen
>
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models


-- 
Paul V. Debes
DFG Research Fellow

Division of Genetics and Physiology
Department of Biology
University of Turku
20014 Finland

Email: paul.debes at utu.fi