[R-sig-ME] Correlations among random variables
Viechtbauer, Wolfgang (SP)
wo||g@ng@v|echtb@uer @end|ng |rom m@@@tr|chtun|ver@|ty@n|
Thu Jan 24 17:42:52 CET 2019
For the data Avi is working with, the default optimizer (nlminb) fails. Switching to 'optim' (with method 'BFGS') works. Here is the fully reproducible code:
df <- read.csv("https://raw.githubusercontent.com/avi-kluger/RCompanion4DDABook/master/Chapter%2010/Chapter10_df.csv")
library(nlme)
df$focalcode <- 1 - df$focalcode
df$partcode <- 1 - df$partcode
### overparamterized model
res1 <- lme(outcome ~ 0 + focalcode + partcode, random = ~ 0 + focalcode + partcode | focalid/dyadid, data = df)
summary(res1)
### contrain sigma to a very small value
res2 <- lme(outcome ~ 0 + focalcode + partcode, random = ~ 0 + focalcode + partcode | focalid/dyadid, data = df, control=list(sigma=1e-8, opt="optim"))
summary(res2)
Just for fun, I also fitted the same model using 'metafor'. While it was not really made for analyzing raw data like this, it can be used to fit the same model (with the devel version) and then sigma can be constrained exactly to 0:
devtools::install_github("wviechtb/metafor")
library(metafor)
df$dyadid.in.focalid <- interaction(df$focalid, df$dyadid)
res3 <- rma.mv(outcome ~ 0 + focalcode + partcode, V=0, random = list(~ 0 + focalcode + partcode | focalid, ~ 0 + focalcode + partcode | dyadid.in.focalid), struct="GEN", data = df, sparse=TRUE)
res3
(note that 'focalid/dyadid' doesn't work at the moment, so you have to create the nested factor manually first; also, model fitting can be slow with rma.mv(), so you might have to wait a bit for it to converge)
The results for res2 and res3 are quite close.
Best,
Wolfgang
-----Original Message-----
From: R-sig-mixed-models [mailto:r-sig-mixed-models-bounces using r-project.org] On Behalf Of Ben Bolker
Sent: Thursday, 24 January, 2019 16:30
To: r-sig-mixed-models using r-project.org
Subject: Re: [R-sig-ME] Correlations among random variables
Not in lme4:
http://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#setting-residual-variances-to-a-fixed-value-zero-or-other
In lme you can set the residual std dev to a fixed small value, but
not to exactly zero:
lme(Reaction~Days,random=~1|Subject,sleepstudy,control=list(sigma=1e-8))
If you're trying to test that the covariance is significantly
positive, I think getting the standard error is the wrong approach; the
Wald (quadratic) approximation is often very bad for random-effects
variances and covariances. I would suggest profile confidence intervals
or a likelihood ratio test.
On 2019-01-23 11:39 p.m., Avraham Kluger wrote:
> Hi again,
>
> A reminder, I am trying to get a stable estimate of the correlation between two random errors in a mode like:
>
> lme(outcome ~ 0 + focalcode + partcode, random = ~ 0 + focalcode + partcode|focalid/dyadid, data = x)
>
> I figured out that five different methods yield identical covariance estimate, as I show below. My new question is CAN YOU CONSTRAIN RESIDUAL TO ZERO in either lme4 or nlme? This may allow convergence and hence estimate of standard error for the covariance.
> -------------------------------------------------------------------------------------------------------------------------------------------------
> Method Var Focal Var Partner Residual Cov Corr
> --------------------------------------------------------------------------------------------------------------------------------------------------
> SPSS .549 .423 .239
> lavaan (SEM) .552 .425 .116 .239
>
> lme4 .326 .200 .223 .115 .451
> nlme .376 .502 .046 .265
> lavaan (SEM) with latent residuals .331 .457 .094 .116 .298
> -------------------------------------------------------------------------------------------------------------------------------------------------
> Some methods printout the covariance and some do not (or fail to converge properly), but all suggest that the covariance is .115 or .116. For MLM in R:
>
> lme4, cov = .451 * sqrt(.326 * .200) = .115
> nlme, cov = .265 * sqrt(.376 * .502) = .115
>
> for SPSS that runs without any warning, cov = .239 * sqrt(.549 * .423) = .115
>
> I appreciate your patience with a novice,
>
> Best,
>
> Avi
More information about the R-sig-mixed-models
mailing list