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

[Avraham N Kluger]

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]]