[R-sig-ME] Correlations among random variables

[Avraham Kluger]

Dear Wolfgang,

Your CIs seemed to have worked in the past.  I tried to replicated it (attached), but I get 

Warning messages:
1: In confint.rma.mv(res3, rho = 1, verbose = TRUE) :
  Cannot obtain lower bound of profile likelihood CI due to convergence problems.
2: In confint.rma.mv(res3, rho = 1, verbose = TRUE) :
  Cannot obtain upper bound of profile likelihood CI due to convergence problems.

I made sure that I have the most recent version of metafor with
devtools::install_github("wviechtb/metafor", force = TRUE)

to no avail.   I also tried to understand how to use control for confint.

I found out yesterday about your keynote in Dubrovnik on May 31st.  I will see if I can attend.

Best,

Avi




-----Original Message-----
From: Viechtbauer, Wolfgang (SP) [mailto:wolfgang.viechtbauer using maastrichtuniversity.nl] 
Sent: Thursday, January 24, 2019 11:18 PM
To: Avraham Kluger <avik using savion.huji.ac.il>; r-sig-mixed-models using r-project.org
Subject: RE: [R-sig-ME] Correlations among random variables

Happy to hear that.

Yes, you can get profile likelihood CIs with:

ci.rho <- confint(res3, rho=1, verbose=TRUE) ci.rho ci.phi <- confint(res3, phi=1, verbose=TRUE) ci.phi

But this will take quite some time, so you might want to run this when you don't need your computer for a while. Here is what I get:

    estimate  ci.lb  ci.ub 
rho   0.6988 0.3660 1.0000

    estimate  ci.lb  ci.ub 
phi   0.2390 0.1465 0.3274

Best,
Wolfgang

-----Original Message-----
From: Avraham Kluger [mailto:avik using savion.huji.ac.il]
Sent: Thursday, 24 January, 2019 20:24
To: r-sig-mixed-models using r-project.org
Cc: Viechtbauer, Wolfgang (SP)
Subject: Re: [R-sig-ME] Correlations among random variables

Dear Wolfgang,

Your metafor solution beautifully replicates, to the dot, results from SPSS and lavann.  Can you obtain CI around the estimates, rho, and phi?

Avi
--------------------------------------------------

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