APIM for Distinguishable Dyads
Load the Data
Load packages and read in the pairwise dataset
library(dyadr)
library(dplyr)
library(nlme)
library(ppcor)
riggsp <- read.csv("riggsp.csv")Distingushable Dyads
A run examining the effect of own Attachment Anxiety (Anxiety_A) and the effect of partner’s Attachment Anxiety (Anxiety_P) on satisfaction (Sat_A) separately for men and women (Gender_A or Genderstring), allowing for nonindependence by correlating the errors of the two members.
There are two approaches to running the distinguisable dyads APIM:
- Interaction Approach
- Two Intercept Approach
Interaction Approach
Steps for using the interaction approach:
- Add distinguishing variable as a covariate. Note it should be a integer, number, or double type variable. Check any variable’s type by running the
class()function. For example,class(riggsp$Gender_A). Here we have effects coded Gender_A (men = -1 and women = 1), but you could also dummy code, just be sure to note how you’ve done it.
- Have the distinguishing variable interact with the actor and the partner effects. You can do this using the
*symbol, or alternatively, the:symbol. Using*will force R to add the main effects, while:will only give you that specific interaction you’re asking for.
- These interactions evaluate whether actor and partner effects are the same for the two people.
- Add a
weights =option the thegls()function to allow for different error variances for the two members. We are using the factor variable,Genderstring, here to make the output easier to read.
apim_di_int <- gls(Sat_A ~ Anxiety_A + Anxiety_P + Gender_A
+ Anxiety_A*Gender_A + Anxiety_P*Gender_A,
data = riggsp,
correlation = corCompSymm(form=~1|Dyad),
weights = varIdent(form=~1|Genderstring),
na.action = na.omit)
#This model is equivalent to the one above.
apim_di_int <- gls(Sat_A ~ Anxiety_A*Gender_A
+ Anxiety_P*Gender_A,
data = riggsp,
correlation = corCompSymm(form=~1|Dyad),
weights = varIdent(form=~1|Genderstring),
na.action = na.omit)
smallsummary(apim_di_int)## Correlation structure of class corCompSymm representing
## Rho
## 0.6157989
##
## Variance function structure of class varIdent representing
## Man Woman
## 1.000000 1.028966
## Residual standard error: 6.638667
##
## Value Std.Error t-value p-value
## (Intercept) 53.3577 1.5079 35.3850 0.0000
## Anxiety_A -1.7994 0.2974 -6.0504 0.0000
## Gender_A 0.4939 0.7356 0.6714 0.5025
## Anxiety_P -1.2696 0.2969 -4.2754 0.0000
## Anxiety_A:Gender_A -0.1189 0.3356 -0.3544 0.7233
## Gender_A:Anxiety_P -0.0457 0.3352 -0.1365 0.8915## 2.5 % 97.5 %
## (Intercept) 50.4022 56.3132
## Anxiety_A -2.3823 -1.2165
## Gender_A -0.9479 1.9356
## Anxiety_P -1.8516 -0.6876
## Anxiety_A:Gender_A -0.7767 0.5388
## Gender_A:Anxiety_P -0.7027 0.6112Interpretation of Effects
Intercept = 53.36—The predicted score for men and women who have a 0 on Anxious Attachment and their partners also have 0 for Anxious Attachment. (We should have centered!)
Gender_A = .49—Women are slightly more satisfied (about 1 point more) than men when you control for both members’ anxious attachment levels. (Recall women are 1 on Gender_A and men are -1; the difference between men and women is then 2x the effect of Gender_A.)
Axiety_A = -1.80—Actor Effect: The more anxiously attached you are, the less satisfied you are (controling for partner’s anxious attachment).
Anxiety_P = -1.27—Partner Effect: The more anxiously attached your partner is, the less satisfied you are (controling for your own anxious attachment).
Anxiety_A:Gender_A = -0.12—The actor effect is more strongly negative for women than for men (although not significantly so).
Gender_A:Anxiety_P = -0.05—The partner effect is more strongly negative for men->women than for women->men (although not significantly so).
When interpreting it may be helpful to compute the separate actor and partner effects for men and women.
Actor Effect for Women = -1.80 + (-0.12) = -1.92
Actor Effect for Men = -1.80 - (-0.12) = -1.68
Partner Effect for M -> W = -1.27 + (-0.05) = -1.32
Partner Effect for W -> M = -1.27 - (-0.05) = -1.22
Or have R do it like so:
sep_coef <- data.frame(
anx_a_men = coef(apim_di_int)[2] + -1*coef(apim_di_int)[5],
anx_a_wom = coef(apim_di_int)[2] + 1*coef(apim_di_int)[5],
anx_p_men = coef(apim_di_int)[4] + -1*coef(apim_di_int)[6],
anx_p_wom = coef(apim_di_int)[4] + 1*coef(apim_di_int)[6],
row.names = "Randi's nums")
sep_coef## anx_a_men anx_a_wom anx_p_men anx_p_wom
## Randi's nums -1.680476 -1.91837 -1.223819 -1.315314Getting Separate Error Variances
Notice is the model summary that we now get weights for the error variances. That is, we actually have two separate error variances. We can use the function getevd in the dyadr package to retrieve them.
error_vars <- data.frame(ev_wom = getevd(apim_di_int)[1],
ev_men = getevd(apim_di_int)[2])
error_vars## ev_wom ev_men
## Group 1 6.830963 6.6386676.83 is the error variance for women, and 6.64 is the error variance for men.
Two-Intercept Approach
This involves a trick by which one equation becomes two. We create two dummy variables: \(M_{ij}\) which equals 1 for men and 0 for women and \(W_{ij}\) which equals 1 for women and zero for men. We then estimate the following equation:
\[Y_{ij} = b_MM_{ij} + a_MM_{ij}A_{ij} + p_MM_{ij}P_{ij} + M_{ij}e_{ij} + b_WW_{ij} + a_WW_{ij}A_{ij} + p_WW_{ij}P_{ij} + W_{ij}e_{ij}\]
Note that the equation has no ordinary intercept, but rather, in some sense, two intercepts, \(b_M\) and \(b_W\). Note that when \(M_{ij} = 1\) and \(W_{ij} = 0\), the above becomes
\[Y_{ij} = b_M + a_MA_{ij} + p_MP_{ij} + e_{ij}\]
and when \(M_{ij} = 0\) and \(W_{ij} = 1\), the above becomes
\[Y_{ij} = b_W + a_WA_{ij} + p_WP_{ij} + e_{ij}\]
Thus, one equation becomes two and we have actor and partner for both members.
To implement this in R, we do the following:
- Add a distinguishing variable that is a factor type variable. We’ll use
Genderstringto make the output nice. - Have no intercept in the fixed model by adding
-1to the end of the formula. - Have the distinguishing variable (
Genderstring) interact with actor and partner effect, but no actor and partner main effects. We need to use:for this instead of*. Separate actor and partner effects will be estimated for each member type. - Keep the
weights =argument to allow for different error variances for the two members.
apim_di_two <- gls(Sat_A ~ Genderstring + Anxiety_A:Genderstring
+ Anxiety_P:Genderstring -1,
data = riggsp,
correlation = corCompSymm(form=~1|Dyad),
weights = varIdent(form=~1|Genderstring),
na.action = na.omit)
smallsummary(apim_di_two)## Correlation structure of class corCompSymm representing
## Rho
## 0.6157989
##
## Variance function structure of class varIdent representing
## Man Woman
## 1.000000 1.028966
## Residual standard error: 6.638667
##
## Value Std.Error t-value p-value
## GenderstringMan 52.8639 1.6536 31.9680 0.0000
## GenderstringWoman 53.8516 1.7015 31.6486 0.0000
## GenderstringMan:Anxiety_A -1.6805 0.4310 -3.8989 0.0001
## GenderstringWoman:Anxiety_A -1.9184 0.4652 -4.1241 0.0000
## GenderstringMan:Anxiety_P -1.2238 0.4521 -2.7072 0.0072
## GenderstringWoman:Anxiety_P -1.3153 0.4435 -2.9658 0.0033## 2.5 % 97.5 %
## GenderstringMan 49.6228 56.1049
## GenderstringWoman 50.5166 57.1865
## GenderstringMan:Anxiety_A -2.5252 -0.8357
## GenderstringWoman:Anxiety_A -2.8301 -1.0067
## GenderstringMan:Anxiety_P -2.1099 -0.3378
## GenderstringWoman:Anxiety_P -2.1845 -0.4461Notice that these actor and partner estimates are the same as what we calculated above. But now we have p-values for them. So that’s cool.
rbind(sep_coef,
model = coef(summary(apim_di_two))[3:6,1],
p_values = coef(summary(apim_di_two))[3:6,4])## anx_a_men anx_a_wom anx_p_men anx_p_wom
## Randi's nums -1.680475614 -1.918370e+00 -1.223819012 -1.315314348
## model -1.680475614 -1.918370e+00 -1.223819012 -1.315314348
## p_values 0.000118932 4.807799e-05 0.007169707 0.003258215Effect sizes
Actor effect size for Men
riggsp %>%
filter(Genderstring == "Man") %$%
pcor.test(Anxiety_A, Sat_A, Anxiety_P)[1] %>%
as.numeric()## [1] -0.3015266Actor effect size for women
riggsp %>%
filter(Genderstring == "Woman") %$%
pcor.test(Anxiety_A, Sat_A, Anxiety_P)[1] %>%
as.numeric()## [1] -0.3172312Partner effect size for Men
riggsp %>%
filter(Genderstring == "Man") %$%
pcor.test(Anxiety_P, Sat_A, Anxiety_A)[1] %>%
as.numeric()## [1] -0.214471Partner effect size for women
riggsp %>%
filter(Genderstring == "Woman") %$%
pcor.test(Anxiety_P, Sat_A, Anxiety_A)[1] %>%
as.numeric()## [1] -0.2338862R2 for Men and Women
We could also get R2 for men and women separately.
apim_di_empty <- gls(Sat_A ~ Gender_A,
data = riggsp,
correlation = corCompSymm(form=~1|Dyad),
weights = varIdent(form=~1|Genderstring),
na.action = na.omit)
R_sqs <- data.frame(R_sq_wom = crsp(getevd(apim_di_two)[1],
getevd(apim_di_empty)[1]),
R_sq_men = crsp(getevd(apim_di_two)[2],
getevd(apim_di_empty)[2]))
R_sqs## R_sq_wom R_sq_men
## Group 1 0.1621929 0.1431408R2 for women is 0.16, or 16% of the variance in satisfaction for women is explained by the model.
R2 for men is 0.14, or 14% of the variance in satisfaction for men is explained by the model.
Descriptive Statistics by Gender
There are MANY ways to do this. I use the dyplr package.
riggsp %>%
select(Genderstring, Anxiety_A, Sat_A) %>%
na.omit() %>% #omits only if missing on the variables included in select()
group_by(Genderstring) %>%
summarise(Anxiety_A_mean = mean(Anxiety_A),
Anxiety_A_sd = sd(Anxiety_A),
Sat_A_mean = mean(Sat_A),
Sat_A_sd = sd(Sat_A),
n = n())## # A tibble: 2 x 6
## Genderstring Anxiety_A_mean Anxiety_A_sd Sat_A_mean Sat_A_sd n
## <fct> <dbl> <dbl> <dbl> <dbl> <int>
## 1 Man 2.63 1.27 44.9 7.17 155
## 2 Woman 2.93 1.21 44.8 7.46 155Model Diagnostics
We should also check that our assumptions are not violated.
plot(apim_di_int)
qqnorm(apim_di_int$residuals)
Looks pretty good to me!