Person Perception: Building the two-level crossed model

# load in our example data
rm_df <- read.csv("/Users/kareenadelrosario/Downloads/DyadicLongitudinal_ExampleDataset_RM3.csv")

Dataset overview

482 stranger dyads interacted during four weekly guided conversations. People reported on their intellectual humility (IH) at three time points: pre (before the first conversation), post (after the last conversation), and follow-up (one month after the last conversation). Three measures of IH per person within each dyad. We are going to use these data to ask whether IH changes across the three time points.

Here are the variables we’re interested in:

• IH: intellectual humility, reported at three time points
• time: time, pre = 1; post = 2; follow-up = 3
• dyad: unique ID number for each dyad
• id: unique ID number for each person
• timecontrast1_pre_post: pre = 0, post = 1, follow-up = 0
• timecontrast2_pre_followup: pre = 0, post = 0, follow-up = 1
• party: -1 republican, 1 democrat (fyi totally made up for the sake of this demonstration)
• RM: a variable coded 1 through 6 for the two dyad members at each time point (see below)

Democrats • 1: pre • 2: post • 3: follow-up

Republicans • 4: pre • 5: post • 6: follow-up

---

Analysis steps

Before we jump into the analysis, let’s take a step back and go through some pre-analyses to figure out:

1. whether we need to account for nonindependence in the first place
2. determine if our dyads are distinguishable
3. build our model

---

Intraclass Correlation Coefficient (ICC)

Sometimes used to “justify” dyadic analyses. ICC indexes the similarity between dyad members for a particular outcome. Technically, it represents the proportion of the total variance in an outcome that is explained by between-dyad variability in mean levels.

Calculating an ICC:

1. Run an MLM with no predictors (i.e., a null model)
2. Use dyad as a random factor
3. Using the output, calculate the ICC as the following:

Run null model with nlme

library(haven) # read SPSS and SAS files
library(reghelper) # for ICC command


null_model_lme <- lme(IH ~ 1,
             data = rm_df,
             random = ~ 1 | dyad,
             method = "REML",
             na.action = na.omit)

summary(null_model_lme)

## Linear mixed-effects model fit by REML
##   Data: rm_df 
##        AIC     BIC    logLik
##   5971.786 5988.81 -2982.893
## 
## Random effects:
##  Formula: ~1 | dyad
##         (Intercept)  Residual
## StdDev:   0.4199118 0.8951192
## 
## Fixed effects:  IH ~ 1 
##               Value  Std.Error   DF  t-value p-value
## (Intercept) 5.71142 0.02741095 1672 208.3627       0
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -5.1820286 -0.4226934  0.1584980  0.6290993  2.0900418 
## 
## Number of Observations: 2154
## Number of Groups: 482

ICC(null_model_lme)

## [1] 0.1803728

Run null model with lme4

library(lme4)

null_model_lmer <- lmer(IH ~ 1 + 
                    (1 | dyad), 
                  data=rm_df)

summary(null_model_lmer)

## Linear mixed model fit by REML ['lmerMod']
## Formula: IH ~ 1 + (1 | dyad)
##    Data: rm_df
## 
## REML criterion at convergence: 5965.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.1820 -0.4227  0.1585  0.6291  2.0900 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  dyad     (Intercept) 0.1763   0.4199  
##  Residual             0.8012   0.8951  
## Number of obs: 2154, groups:  dyad, 482
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  5.71142    0.02741   208.4

Here’s the ICC equation:

\[ICC = \frac {\ between-dyad\ variance\ (random\ intercept\ variance)} {total\ variance\ (random \ intercept\ variance\ +\ residual\ variance)} \]

Dyad_variance <- 0.1763
residual <- 0.8012

# calculating ICC
Dyad_variance / (Dyad_variance + residual)

## [1] 0.1803581

What is our takeaway from the ICC?

---

Testing distinguishability

Now, we need to determine whether our dyads are distinguishable or indistinguishable. We can empirically test whether our ‘distinguishing factor’ (e.g., gender, role) matters. To do this, we’re going to use an omnibus test of distinguishability, which tests whether the distinguishing factor matters across all components of your model.

What’s the difference?

Our predictors are timecontrast1_pre_post, timecontrast2_pre_followup, and party.

If we treat the dyad as distinguishable, we’re estimating 8 parameters:

1. Predictors: timecontrast1_pre_post, timecontrast2_pre_followup, party, timecontrast1_pre_post x party, timecontrast2_pre_followup x party
2. Error variance for democrats, error variance for republicans, and the covariance between errors

If we treat the dyad as indistinguishable, we’re estimating 4 parameters:

1. Predictors: timecontrast1_pre_post, timecontrast2_pre_followup
2. One error variance for both dyad members and the covariance between errors

Now onto the test

1. Run the two models using ML, instead of the default REML
2. Conduct a χ2 difference test
3. Subtract the deviances of the two models (called “-2 Log Likelihood”)
4. Degrees of freedom for the test = difference between the number of parameters in the two models
5. Look up p value for the difference with the appropriate df

First, the distinguishable model

distinguishable_model_ML <- gls(IH ~ timecontrast1_pre_post + timecontrast2_pre_followup +
                                  party+
                                  timecontrast1_pre_post:party + 
                                  timecontrast2_pre_followup:party,
                                na.action=na.omit,
                                method="ML",
                                correlation = corCompSymm (form=~1|dyad),
                                weights = varIdent(form=~1|party),
                                data=rm_df)

Now, the indistinguishable model

indistinguishable_model_ML <- gls(IH ~ timecontrast1_pre_post + timecontrast2_pre_followup +
                                  party+
                                  timecontrast1_pre_post:party + 
                                  timecontrast2_pre_followup:party,
                                na.action=na.omit,
                                method="ML",
                                correlation = corCompSymm (form=~1|dyad),
                                data=rm_df)

# compare the two models
anova(distinguishable_model_ML, indistinguishable_model_ML)

##                            Model df      AIC      BIC   logLik   Test   L.Ratio
## distinguishable_model_ML       1  9 5878.319 5929.395 -2930.16                 
## indistinguishable_model_ML     2  8 5877.261 5922.662 -2930.63 1 vs 2 0.9416936
##                            p-value
## distinguishable_model_ML          
## indistinguishable_model_ML  0.3318

What’s the takeaway?

---

Breaking down the two-level crossed model

This will be a marginal model– we’re going to adjust for nonindependence by providing a variance-covariance structure to the residuals. This approach is useful for repeated measure with 3 or 4 time points.

In other words, this will model the nonindependence with a covariance (adjust for nonindependence by correlating the residuals between dyad members. This allows for both positive and negative nonindependence).

Output

1. Fixed parameter estimates
2. Error covariance: are the two dyad members’ errors related ot each other?
3. Error variance: is there significant variance that is unexplained by the other components in the model?

Now, in R, we can only run a two-level crossed model with distinguishable dyads. Why you ask? Well, with indistinguishable dyads, we need to pool the variances which requires us to customize the variance-covariance matrix.

MIXED IH WITH party timecontrast1_pre_post timecontrast2_pre_followup
/FIXED = party timecontrast1_pre_post timecontrast2_pre_followup
party*timecontrast1_pre_post party*timecontrast2_pre_followup |
SSTYPE(3)
/METHOD = REML
/PRINT = G TESTCOV SOLUTION
/REPEATED RM| subject(dyad) COVTYPE(UN).

Two-level crossed model

model <- gls(IH ~ timecontrast1_pre_post + timecontrast2_pre_followup + party + 
               party*timecontrast1_pre_post +
                party*timecontrast2_pre_followup,
             data = rm_df,
             correlation = corSymm(form = ~ RM | dyad),
             weights=varIdent(form=~1|RM),
             method = "REML",
             na.action = na.omit)

# omnibus test
anova(model, type = "marginal")

## Denom. DF: 2148 
##                                  numDF   F-value p-value
## (Intercept)                          1 25041.583  <.0001
## timecontrast1_pre_post               1   115.434  <.0001
## timecontrast2_pre_followup           1    67.907  <.0001
## party                                1     0.072  0.7878
## timecontrast1_pre_post:party         1     0.367  0.5446
## timecontrast2_pre_followup:party     1     0.085  0.7702

SAS Omnibus output

Residual structure (Level 1 random effects)

• 3 variances
  • Variance at pre
  • Variance at post
  • Variance at followup
  
• 3 within-dyad covariances at the same time point
  • Within-dyad covariance between pre ratings
  • Within-dyad covariance between post ratings
  • Within-dyad covariance between follow-up ratings
  
• 3 within-person covariances at different time points
  • Within-person covariance between pre and post
  • Within-person covariance between pre and follow-up
  • Within-person covariance between post and follow-up
  
• 3 within-dyad covariances at different time points
  • Within-dyad covariance between pre and post
  • Within-dyad covariance between pre and follow-up
  • Within-dyad covariance between post and follow-up

Print covariance parameter estimates

getVarCov(model)

## Marginal variance covariance matrix
##             [,1]     [,2]     [,3]        [,4]
## [1,]  1.08010000 0.446180 0.073189 -0.00016797
## [2,]  0.44618000 0.961690 0.049031  0.05365500
## [3,]  0.07318900 0.049031 1.041300  0.43985000
## [4,] -0.00016797 0.053655 0.439850  0.89507000
##   Standard Deviations: 1.0393 0.98066 1.0205 0.94608

getVarCov(model, individual = 2)

## Marginal variance covariance matrix
##             [,1]     [,2]       [,3]      [,4]        [,5]      [,6]
## [1,]  1.08010000 0.446180  0.3666300  0.073189 -0.00016797 0.0633330
## [2,]  0.44618000 0.961690  0.3452600  0.049031  0.05365500 0.0865530
## [3,]  0.36663000 0.345260  0.6338600 -0.049843 -0.10210000 0.0034325
## [4,]  0.07318900 0.049031 -0.0498430  1.041300  0.43985000 0.4533000
## [5,] -0.00016797 0.053655 -0.1021000  0.439850  0.89507000 0.4190200
## [6,]  0.06333300 0.086553  0.0034325  0.453300  0.41902000 0.6540700
##   Standard Deviations: 1.0393 0.98066 0.79616 1.0205 0.94608 0.80874

SAS Covariance output

What can we learn from this?

---

Your turn!

Now, we’ll follow the same steps with the Roommate data we’ve been working with. Here’s what you’ll need to do:

1. Are the data nonindependent?
2. Are the dyads distinguishable?
3. Does wishing you had a new roommate depend on whether they tend to be sloppy, whether you thought you guys were similar when you first met, your gender, and day?

As a reminder, here are the variables we’re playing with: awish sloppy_roommateC first_similarNC genderN dayR Dyad RM partnum

Before we start, load the sample data_1.sav file from the Google drive -> week2 -> datafiles. Note this is a different file than what we’ve been using.

library(haven)

roommate_df <- read_sav("/Users/kareenadelrosario/Desktop/everything/2_areas/spring 2024/person perception/sample data_1.sav")

In this analysis, our grouping/subject variable is “Dyad” and our repeated variable is “RM”. For reference, here is the SPSS syntax:

MIXED Awish WITH sloppy_roommateC first_similarNC genderN dayR 
 /FIXED = sloppy_roommateC first_similarNC genderN dayR 
 /PRINT = SOLUTION TESTCOV
 /REPEATED =RM | SUBJECT (Dyad) COVTYPE (UNR).

For R, the equivalent to COVTYPE (UNR) is correlation = corSymm.

For reference, here are the results you should get

Fixed effects

Covariance parameters

• note that in R, we’re getting the covariances, not correlations. Only compare the variances.