[R-sig-ME] Question about flexible LME4 Variance-Covariance Structure

Phillip Alday me @end|ng |rom ph||||p@|d@y@com
Thu Aug 12 02:26:31 CEST 2021


The problem is

(0 | BlockingName)

The left-hand side of | doesn't have any coefficients (because 0 is
suppressing the intercept) and that's what's causing the error in lme4.

I don't have time to think about the statistical problem at the moment,
but maybe knowing the software issue already helps a bit.

Phillip

On 11/8/21 6:19 pm, Jared Anderson wrote:
> Hello All,
> 
> I am a Data Scientist with Bayer Crop Science and I deal often with generalized linear mixed models in my role. I am new to this subscription group and wanted to reach out because I am trying to fit two different kinds of models on the same data and use their residual variances as a statistic for assessing data quality.
> For our variance-ratio statistic, two slightly different linear models are used. Both models are simple single treatment factor randomized complete block design models. The difference is in the assumed covariance structure of the residuals.
> The first model makes the standard assumption that the residuals all have the same variance, while the second model assumes that the residuals variances differ from rep to rep.
> Model 1: y = m + T + R + e, where m = overall mean, T = treatment effect, R = blocking effect and e = random error which is assumed to be homoschedastic: cov(e) = common variance*I = Diag(common variance).
> Model 2: y = m + T + R + e, where m = overall mean, T = treatment effect, R = blocking effect and e = random error which is assumed to be heteroschedastic: cov(e) = D(var[i]) i=1, …, number of reps. The variance of the residuals is different for each rep.
> VR[i] = var[i]/common variance, one VR value for each rep.
> 
> I have been able to build the desired Model1 in LME4
> Formula_LME4 <- as.formula(paste0(zParams$value,” ~ 1 + (1|”,zParams$factorLevelId,”) + (1|”BlockingName,”)”))
> HomoS_LME4 <- lmer(Formula_LME4,data=zdf) #covariance structure is the same for all reps
> (Treating the fixed effects (factorLevelId) as random is a way to trick the routines into calculating the needed variances.)
> 
> 
> Where I am running into a problem and was hoping to get your advice: I cannot seem to implement Model 2 successfully
> I think the LME4 Heteroscedastic formula is supposed to be this, if I understand it correctly:
> lmer(NUM_VALUE ~ 1 + (1 | TREATMENT) + (1 | BlockingName) + (0 | BlockingName),data=zdf)
> But it is throwing an error:
> Linear mixed model fit by REML ['lmerMod']
> Formula: NUM_VALUE ~ 1 + (1 | TREATMENT) + (1 | BlockingName) + (0 | BlockingName)
>    Data: zdf
> REML criterion at convergence: 1384.89
> Error in thl[[i]] : subscript out of bounds
> And I was wondering if you’ve ever come across the error above?
> 
> 
> If it helps, back when we leveraged Asreml licenses for our work, I was able to generate Model 2 this way:
> 
> fixedFormula_ASREML <- as.formula(paste0(zParams$value," ~ 1"))
> randomFormula_ASREML <- as.formula(paste0(" ~ ",zParams$factorLevelId," + ", BlockingName))
> covFormula_ASREML    <- as.formula(paste0(" ~ at(", BlockingName, "):units"))
> 
> HeteroS_ASREML <- asreml(fixed  = fixedFormula_ASREML,
>                   random = randomFormula_ASREML,
>                   rcov   = covFormula_ASREML,
>                   data   = zdf)
> 
> 
> Best,
> Jared A.
> 
> 
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>



More information about the R-sig-mixed-models mailing list