[R-sig-ME] How to specify a correlation between cross-nested effects

[Avraham N Kluger]

Hi,



I am trying to model round-robin data, typically analyzed with an ANOVA-based Social Relations Model (SRM), with a mixed model with cross-nested effects.

First, I give a brief background to clarify the motivation behind my mixed model question. One way to analyze SRM data in R is with tripleR package. I use the example data of that package. The data reflect the liking rating each person in a group of ten gave each other. A record with NA is needed in tripleR for self-rating (liking that one would have given to oneself if measured). The model estimates three variances: perceiver (actor), target (partner), and error, and two covariances: perceiver-target and dyadic (error + dyadic). [The perceiver variance reflects differences in liking others; the target variance reflects differences in being liked by others; and the perceiver-target covariance reflects the degree to which people who like others are liked by others. For my purpose, I ignore the dyadic covariance in my question].

Below is the code for the tripleR example (using Group 1 data only for simplicity)

library(TripleR)

library(glmmTMB)

data("multiLikingLong")

df <- multiLikingLong

df <- df[df$group.id ==1, ]

RR(liking_a ~ perceiver.id*target.id, data = df)



It produces the following output



Round-Robin object ('RR'), calculated by TripleR

------------------------------------------------

Univariate analysis of one round robin variable



Univariate analyses for: liking_a

---------

Round robin analysis for a single group; using the formula of Lashley & Bond (1997).



                         estimate standardized    se t.value p.value

actor variance              0.228        0.267 0.126   1.814   0.052

partner variance            0.067        0.079 0.058   1.161   0.138

relationship variance       0.558        0.654 0.098   5.684   0.000

error variance                 NA           NA    NA      NA      NA

actor-partner covariance    0.041        0.329 0.067   0.604   0.561

relationship covariance     0.216        0.388 0.098   2.203   0.028

Actor effect reliability: .777

Partner effect reliability: .506



Next, I get rid of the rows with NA, needed by tripleR, and run a mixed model.



x <- na.omit(df)

m <- glmmTMB(liking_a ~ 1 + (1 |perceiver.id) + (1 |target.id), data = x,

             family = gaussian(link = "identity"))

summary(m)



It yields the following results:



Family: gaussian  ( identity )

Formula:          liking_a ~ 1 + (1 | perceiver.id) + (1 | target.id)

Data: x



     AIC      BIC   logLik deviance df.resid

   231.7    241.7   -111.8    223.7       86



Random effects:



Conditional model:

Groups       Name        Variance Std.Dev.

perceiver.id (Intercept) 0.19986  0.4471

target.id    (Intercept) 0.06102  0.2470

Residual                 0.55838  0.7472

Number of obs: 90, groups:  perceiver.id, 10; target.id, 10



Dispersion estimate for gaussian family (sigma^2): 0.558



Conditional model:

            Estimate Std. Error z value Pr(>|z|)

(Intercept)   3.3556     0.1797   18.67   <2e-16 ***

---

Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1



The results for the perceiver and target variances are similar but not identical. The results for the residuals are identical. For example, tripleR estimates the perceiver variance as .228 but glmmTMB, .199. Yet, tripleR also estimates the actor-partner covariance and correlation; the standardized correlation is .329. My question is, how can I specify this correlation with glmmTMB?



I found an approximation for the correlation:



cor(

    ranef(m)$cond$perceiver.id,

    ranef(m)$cond$target.id

    )



            (Intercept)

(Intercept)   0.2959082



But, this approximation is based on the random effects whose variances differ from the model's.

Thank you,



Avi Kluger

http://avikluger.wix.com/avi-kluger


	[[alternative HTML version deleted]]