Variability in within-participant avoidance, anxiety, and relationship satisfaction by hierarchy type
Step 1: Look at hierarchy type only
To test whether the within-participant variability in avoidance, anxiety, and relationship satisfaction differed by hierarchy type, we compared (via likelihood ratio test) two linear mixed effects models for each outcome variable: one which assumed a homogeneous variance across hierarchical and non-hierarchical groups and one which assumed heterogeneous variance. Both models had a random intercept and no fixed effects.
Options used in the nlme::lme() function throughout this analysis:
random = list(id = pdSymm(form = ~ 1)), which is the same as random = ~ 1|id, sets the random effects to be structured as a positive-definite symmetric matrix; it assumes homogeneous between-subject variance.
random = list(id = pdDiag(form = ~ iv)) sets the random effects to be structured as a diagonal matrix; it assumes heterogeneous between-subject variance across the levels of iv (hierarchical and non-hierarchical groups).
weights = varIdent(form = ~ 1 | iv) sets heterogeneous within-participant variance for levels of iv (hierarchical and non-hierarchical groups)
method = "ML" Maximum Likelihood is used to compare models and method = "REML" Restricted ML on final model.
wph1 <- lme(dv1 ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~ 1 | iv),
method="ML")
wph2 <- lme(dv2 ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~ 1 | iv),
method="ML")
wph3 <- lme(dv3 ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~ 1 | iv),
method="ML")
a1 <- anova(ri1,wph1); print(a1)
## Model df AIC BIC logLik Test L.Ratio p-value
## ri1 1 3 1698.657 1711.663 -846.3287
## wph1 2 4 1689.768 1707.109 -840.8843 1 vs 2 10.8889 0.001
a2 <- anova(ri2,wph2); print(a2)
## Model df AIC BIC logLik Test L.Ratio p-value
## ri2 1 3 1921.240 1934.245 -957.6200
## wph2 2 4 1912.517 1929.858 -952.2586 1 vs 2 10.72264 0.0011
a3 <- anova(ri3,wph3); print(a3)
## Model df AIC BIC logLik Test L.Ratio p-value
## ri3 1 3 2971.427 2984.433 -1482.714
## wph3 2 4 2969.064 2986.405 -1480.532 1 vs 2 4.363152 0.0367
pvals <- c(pvals,
a1$`p-value`, a2$`p-value`, a3$`p-value`)
Hierarchical relationships are associated with higher within-participant variability in avoidance, anxiety, and relationship satisfaction. However, after adjusting for multiple testing, the effect for relationship satisfaction is not significant.
summary(wph1)
## Linear mixed-effects model fit by maximum likelihood
## Data: NULL
## AIC BIC logLik
## 1689.769 1707.109 -840.8843
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 0.39479 0.9110345
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | iv
## Parameter estimates:
## NonHierarchical Hierarchical
## 1.000000 1.251505
## Fixed effects: dv1 ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) 2.06929 0.05016391 339 41.25057 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.7100077 -0.7452329 -0.2186775 0.5639658 4.6224882
##
## Number of Observations: 564
## Number of Groups: 225
summary(wph2)
## Linear mixed-effects model fit by maximum likelihood
## Data: NULL
## AIC BIC logLik
## 1912.517 1929.858 -952.2586
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 0.770111 1.009911
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | iv
## Parameter estimates:
## NonHierarchical Hierarchical
## 1.00000 1.26607
## Fixed effects: dv2 ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) 2.227532 0.07048204 339 31.60425 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.4722788 -0.5619507 -0.3110696 0.5062236 3.4688398
##
## Number of Observations: 564
## Number of Groups: 225
summary(wph3)
## Linear mixed-effects model fit by maximum likelihood
## Data: NULL
## AIC BIC logLik
## 2969.064 2986.405 -1480.532
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 1.709768 2.775777
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | iv
## Parameter estimates:
## NonHierarchical Hierarchical
## 1.000000 1.158351
## Fixed effects: dv3 ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) 16.04599 0.170797 339 93.94772 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.6383066 -0.5963369 0.1176959 0.6894816 1.9773680
##
## Number of Observations: 564
## Number of Groups: 225
The within-participant standard deviations for avoidance, anxiety, and relationship satisfaction were higher in the hierarchical group by 25%, 26%, and 16% respectively.
Estimated standard deviations for non-hierarchical and hierarchical groups:
summary(wph1)$sigma # avoidance, non-hierarchical
## [1] 0.9110345
summary(wph1)$sigma*coef(wph1$modelStruct$varStruct, uncons=FALSE)
## Hierarchical
## 1.140164
summary(wph2)$sigma # anxiety, non-hierarchical
## [1] 1.009911
summary(wph2)$sigma*coef(wph2$modelStruct$varStruct, uncons=FALSE)
## Hierarchical
## 1.278618
summary(wph3)$sigma # satisfaction, non-hierarchical
## [1] 2.775777
summary(wph3)$sigma*coef(wph3$modelStruct$varStruct, uncons=FALSE)
## Hierarchical
## 3.215325
Confidence intervals:
intervals(wph1, level = 0.95, which = "all")
## Approximate 95% confidence intervals
##
## Fixed effects:
## lower est. upper
## (Intercept) 1.970705 2.06929 2.167874
## attr(,"label")
## [1] "Fixed effects:"
##
## Random Effects:
## Level: id
## lower est. upper
## sd((Intercept)) 0.2832421 0.39479 0.5502682
##
## Variance function:
## lower est. upper
## Hierarchical 1.095998 1.251505 1.429076
## attr(,"label")
## [1] "Variance function:"
##
## Within-group standard error:
## lower est. upper
## 0.8316119 0.9110345 0.9980424
intervals(wph2, level = 0.95, which = "all")
## Approximate 95% confidence intervals
##
## Fixed effects:
## lower est. upper
## (Intercept) 2.089018 2.227532 2.366046
## attr(,"label")
## [1] "Fixed effects:"
##
## Random Effects:
## Level: id
## lower est. upper
## sd((Intercept)) 0.6433183 0.770111 0.9218935
##
## Variance function:
## lower est. upper
## Hierarchical 1.10037 1.26607 1.456722
## attr(,"label")
## [1] "Variance function:"
##
## Within-group standard error:
## lower est. upper
## 0.915823 1.009911 1.113665
intervals(wph3, level = 0.95, which = "all")
## Approximate 95% confidence intervals
##
## Fixed effects:
## lower est. upper
## (Intercept) 15.71034 16.04599 16.38165
## attr(,"label")
## [1] "Fixed effects:"
##
## Random Effects:
## Level: id
## lower est. upper
## sd((Intercept)) 1.387745 1.709768 2.106515
##
## Variance function:
## lower est. upper
## Hierarchical 1.008444 1.158351 1.330542
## attr(,"label")
## [1] "Variance function:"
##
## Within-group standard error:
## lower est. upper
## 2.521170 2.775777 3.056096
Exploratory - checked whether a heterogeneous between-subjects variance model improved the heterogeneous within-subjects variance model:
bph1 <- lme(dv1 ~ 1,
random = list(id = pdDiag(form = ~ iv)),
weights = varIdent(form = ~ 1 | iv),
method="ML")
bph2 <- lme(dv2 ~ 1,
random = list(id = pdDiag(form = ~ iv)),
weights = varIdent(form = ~ 1 | iv),
method="ML")
bph3 <- lme(dv3 ~ 1,
random = list(id = pdDiag(form = ~ iv)),
weights = varIdent(form = ~ 1 | iv),
method="ML")
a1 <- anova(wph1,bph1); print(a1)
## Model df AIC BIC logLik Test L.Ratio p-value
## wph1 1 4 1689.768 1707.109 -840.8843
## bph1 2 5 1691.768 1713.444 -840.8843 1 vs 2 3.92165e-07 0.9995
a2 <- anova(wph2,bph2); print(a2)
## Model df AIC BIC logLik Test L.Ratio p-value
## wph2 1 4 1912.517 1929.858 -952.2586
## bph2 2 5 1914.517 1936.192 -952.2584 1 vs 2 0.0005435428 0.9814
a3 <- anova(wph3,bph3); print(a3)
## Model df AIC BIC logLik Test L.Ratio p-value
## wph3 1 4 2969.064 2986.405 -1480.532
## bph3 2 5 2971.064 2992.740 -1480.532 1 vs 2 5.476795e-07 0.9994
pvals <- c(pvals, a1$`p-value`, a2$`p-value`, a3$`p-value`)
No significant improvement.
Step 2: add control variables
We first checked which potential control variables are related to the within-participant variability in the outcomes.
# numerical vars
cvs.num <- c("rel.length.c", "age.c", "incomeperPPL.c")
# categorical vars
cvs.cat <- c("N.Partners.3levels", "gender", "orientation.binary", "race.binary", "children.binary", "education.3levels", "marital.status", "cohab", "coparent")
cvs <- c(cvs.num, cvs.cat)
For numerical variables, we considered 3 options:
weights = varFixed(~cv), which assumes that the variance is proportional to cv;
weights = varPower(1,~cv), which assumes that the variance is a power of cv;
weights = varExp(1,~cv), which assumes that the variance is an exponential of cv.
The option, if any, that improved the model the most was selected.
for(x in cvs.num){ # numerical vars
data <- na.omit(study1.long[c(vars,x)])
print(x)
cv <- data[[x]]
dv1 <- data$avoidance
dv2 <- data$anxiety
dv3 <- data$satisfaction
iv <- data$hierarchy
id <- data$Respondent.ID
m1 <- lme(dv1 ~ 1,
random = list(id = pdSymm(form = ~ 1)),
method="ML")
m2 <- lme(dv2 ~ 1,
random = list(id = pdSymm(form = ~ 1)),
method="ML")
m3 <- lme(dv3 ~ 1,
random = list(id = pdSymm(form = ~ 1)),
method="ML")
mcv1 <- lme(dv1 ~ 1,
random = list(id = pdSymm(form = ~ 1)),
#weights = varFixed(~cv),
#weights = varPower(1,~cv),
weights = varExp(1,~cv),
method="ML",
control=lmeControl(opt = "optim"))
mcv2 <- lme(dv2 ~ 1,
random = list(id = pdSymm(form = ~ 1)),
#weights = varFixed(~cv),
#weights = varPower(1,~cv),
weights = varExp(1,~cv),
method="ML",
control=lmeControl(opt = "optim"))
mcv3 <- lme(dv3 ~ 1,
random = list(id = pdSymm(form = ~ 1)),
#weights = varFixed(~cv),
#weights = varPower(1,~cv),
weights = varExp(1,~cv),
method="ML",
control=lmeControl(opt = "optim"))
a1 <- anova(m1,mcv1);
a2 <- anova(m2,mcv2);
a3 <- anova(m3,mcv3);
print(a1); print(a2); print(a3);
pvals <- c(pvals,
a1$`p-value`, a2$`p-value`, a3$`p-value`)
}
## [1] "rel.length.c"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 3 1682.265 1695.238 -838.1323
## mcv1 2 4 1666.686 1683.983 -829.3430 1 vs 2 17.57869 <.0001
## Model df AIC BIC logLik Test L.Ratio p-value
## m2 1 3 1897.491 1910.464 -945.7455
## mcv2 2 4 1894.945 1912.243 -943.4726 1 vs 2 4.545705 0.033
## Model df AIC BIC logLik Test L.Ratio p-value
## m3 1 3 2943.723 2956.696 -1468.861
## mcv3 2 4 2943.578 2960.876 -1467.789 1 vs 2 2.14462 0.1431
## [1] "age.c"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 3 1698.657 1711.663 -846.3287
## mcv1 2 4 1713.171 1730.511 -852.5854 1 vs 2 12.5134 4e-04
## Model df AIC BIC logLik Test L.Ratio p-value
## m2 1 3 1921.240 1934.245 -957.6200
## mcv2 2 4 1967.276 1984.616 -979.6378 1 vs 2 44.0357 <.0001
## Model df AIC BIC logLik Test L.Ratio p-value
## m3 1 3 2971.427 2984.433 -1482.714
## mcv3 2 4 3007.350 3024.690 -1499.675 1 vs 2 33.92267 <.0001
## [1] "incomeperPPL.c"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 3 1328.222 1340.455 -661.1110
## mcv1 2 4 1336.862 1353.173 -664.4312 1 vs 2 6.640512 0.01
## Model df AIC BIC logLik Test L.Ratio p-value
## m2 1 3 1491.388 1503.620 -742.6938
## mcv2 2 4 1487.526 1503.836 -739.7628 1 vs 2 5.861914 0.0155
## Model df AIC BIC logLik Test L.Ratio p-value
## m3 1 3 2305.779 2318.012 -1149.890
## mcv3 2 4 2334.305 2350.616 -1163.152 1 vs 2 26.5259 <.0001
Categorical variables:
for(x in cvs.cat){ # categorical vars
data <- na.omit(study1.long[c(vars,x)])
print(x)
cv <- data[[x]]
dv1 <- data$avoidance
dv2 <- data$anxiety
dv3 <- data$satisfaction
iv <- data$hierarchy
id <- data$Respondent.ID
m1 <- lme(dv1 ~ 1,
random = list(id = pdSymm(form = ~ 1)),
method="ML")
m2 <- lme(dv2 ~ 1,
random = list(id = pdSymm(form = ~ 1)),
method="ML")
m3 <- lme(dv3 ~ 1,
random = list(id = pdSymm(form = ~ 1)),
method="ML")
mcv1 <- lme(dv1 ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~1|cv),
method="ML")
mcv2 <- lme(dv2 ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~1|cv),
method="ML")
mcv3 <- lme(dv3 ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~1|cv),
method="ML")
a1 <- anova(m1,mcv1)
a2 <- anova(m2,mcv2)
a3 <- anova(m3,mcv3)
print(a1); print(a2); print(a3);
pvals <- c(pvals,
a1$`p-value`, a2$`p-value`, a3$`p-value`)
}
## [1] "N.Partners.3levels"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 3 1698.657 1711.663 -846.3287
## mcv1 2 5 1679.752 1701.427 -834.8760 1 vs 2 22.90553 <.0001
## Model df AIC BIC logLik Test L.Ratio p-value
## m2 1 3 1921.240 1934.245 -957.6200
## mcv2 2 5 1925.191 1946.866 -957.5955 1 vs 2 0.0489621 0.9758
## Model df AIC BIC logLik Test L.Ratio p-value
## m3 1 3 2971.427 2984.433 -1482.714
## mcv3 2 5 2965.826 2987.501 -1477.913 1 vs 2 9.601514 0.0082
## [1] "gender"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 3 1406.073 1418.505 -700.0365
## mcv1 2 5 1396.273 1416.994 -693.1366 1 vs 2 13.79972 0.001
## Model df AIC BIC logLik Test L.Ratio p-value
## m2 1 3 1587.874 1600.307 -790.9372
## mcv2 2 5 1590.274 1610.995 -790.1371 1 vs 2 1.600361 0.4492
## Model df AIC BIC logLik Test L.Ratio p-value
## m3 1 3 2460.656 2473.088 -1227.328
## mcv3 2 5 2462.291 2483.012 -1226.145 1 vs 2 2.365021 0.3065
## [1] "orientation.binary"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 3 1389.050 1401.450 -691.5247
## mcv1 2 4 1386.429 1402.963 -689.2146 1 vs 2 4.620197 0.0316
## Model df AIC BIC logLik Test L.Ratio p-value
## m2 1 3 1562.749 1575.149 -778.3744
## mcv2 2 4 1564.359 1580.893 -778.1797 1 vs 2 0.3895427 0.5325
## Model df AIC BIC logLik Test L.Ratio p-value
## m3 1 3 2429.558 2441.958 -1211.779
## mcv3 2 4 2431.553 2448.087 -1211.777 1 vs 2 0.004521638 0.9464
## [1] "race.binary"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 3 1401.881 1414.300 -697.9404
## mcv1 2 4 1403.879 1420.439 -697.9396 1 vs 2 0.001567308 0.9684
## Model df AIC BIC logLik Test L.Ratio p-value
## m2 1 3 1575.462 1587.882 -784.7312
## mcv2 2 4 1576.899 1593.459 -784.4496 1 vs 2 0.5630493 0.453
## Model df AIC BIC logLik Test L.Ratio p-value
## m3 1 3 2447.338 2459.758 -1220.669
## mcv3 2 4 2448.683 2465.242 -1220.341 1 vs 2 0.6553933 0.4182
## [1] "children.binary"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 3 1406.073 1418.505 -700.0365
## mcv1 2 4 1406.876 1423.453 -699.4380 1 vs 2 1.197015 0.2739
## Model df AIC BIC logLik Test L.Ratio p-value
## m2 1 3 1587.874 1600.307 -790.9372
## mcv2 2 4 1589.062 1605.638 -790.5307 1 vs 2 0.8129924 0.3672
## Model df AIC BIC logLik Test L.Ratio p-value
## m3 1 3 2460.656 2473.088 -1227.328
## mcv3 2 4 2461.101 2477.678 -1226.551 1 vs 2 1.554649 0.2125
## [1] "education.3levels"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 3 1385.941 1398.328 -689.9705
## mcv1 2 5 1389.896 1410.541 -689.9478 1 vs 2 0.04538158 0.9776
## Model df AIC BIC logLik Test L.Ratio p-value
## m2 1 3 1566.216 1578.603 -780.108
## mcv2 2 5 1569.048 1589.693 -779.524 1 vs 2 1.167952 0.5577
## Model df AIC BIC logLik Test L.Ratio p-value
## m3 1 3 2418.121 2430.508 -1206.060
## mcv3 2 5 2421.047 2441.692 -1205.523 1 vs 2 1.074016 0.5845
## [1] "marital.status"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 3 1698.657 1711.663 -846.3287
## mcv1 2 4 1696.037 1713.377 -844.0186 1 vs 2 4.620259 0.0316
## Model df AIC BIC logLik Test L.Ratio p-value
## m2 1 3 1921.240 1934.245 -957.6200
## mcv2 2 4 1919.955 1937.295 -955.9774 1 vs 2 3.285205 0.0699
## Model df AIC BIC logLik Test L.Ratio p-value
## m3 1 3 2971.427 2984.433 -1482.714
## mcv3 2 4 2972.554 2989.894 -1482.277 1 vs 2 0.8738922 0.3499
## [1] "cohab"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 3 1698.657 1711.663 -846.3287
## mcv1 2 5 1680.353 1702.029 -835.1766 1 vs 2 22.30414 <.0001
## Model df AIC BIC logLik Test L.Ratio p-value
## m2 1 3 1921.240 1934.245 -957.6200
## mcv2 2 5 1915.071 1936.746 -952.5355 1 vs 2 10.16901 0.0062
## Model df AIC BIC logLik Test L.Ratio p-value
## m3 1 3 2971.427 2984.433 -1482.714
## mcv3 2 5 2974.534 2996.209 -1482.267 1 vs 2 0.8933916 0.6397
## [1] "coparent"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 3 1698.657 1711.663 -846.3287
## mcv1 2 5 1698.474 1720.149 -844.2368 1 vs 2 4.183865 0.1234
## Model df AIC BIC logLik Test L.Ratio p-value
## m2 1 3 1921.240 1934.245 -957.6200
## mcv2 2 5 1920.436 1942.112 -955.2181 1 vs 2 4.803633 0.0906
## Model df AIC BIC logLik Test L.Ratio p-value
## m3 1 3 2971.427 2984.433 -1482.714
## mcv3 2 5 2974.812 2996.487 -1482.406 1 vs 2 0.6154278 0.7351
After adjusting p-values for multiple comparisons, the following variables were related to within-participant variance:
Avoidance - relationship length, number of partners, gender, and cohabitation.
Anxiety - cohabitation
Satisfaction - number of partners
Model selection for within-participant variability in avoidance
Strategy: start with a model that includes hierarchy (variable of interest), relationship length, number of partners, gender, and cohabitation status. Then perform model selection through a backwards elimination process, that is, only variables that significantly improve a model with fewer variables are left in the final model.
Model with 4 control variables
data <- na.omit(study1.long[c(vars,"rel.length.c","N.Partners.3levels","gender","cohab")])
dv <- data$avoidance
iv <- data$hierarchy
id <- data$Respondent.ID
cv1 <- data[["rel.length.c"]]
cv4 <- data[["N.Partners.3levels"]]
cv5 <- data[["gender"]]
cv11 <- data[["cohab"]]
mcv4 <- lme(dv ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
varIdent(form = ~ 1 | iv),
varExp(1, ~cv1), # length
varIdent(form = ~ 1 | cv4), # n partners
varIdent(form = ~ 1 | cv5), # gender
varIdent(form = ~ 1 | cv11) # cohab
),
method="ML")
aic3 <- numeric() # vector that will store AIC values for later comparison
mcv3 <- lme(dv ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
varIdent(form = ~ 1 | iv),
#varExp(1, ~cv1), # length
varIdent(form = ~ 1 | cv4), # n partners
varIdent(form = ~ 1 | cv5), # gender
varIdent(form = ~ 1 | cv11) # cohab
),
method="ML")
a1 <- anova(mcv3,mcv4)
aic3[1] <- a1$AIC[2]-a1$AIC[1] # AIC of full minus reduced model
mcv3 <- lme(dv ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
varIdent(form = ~ 1 | iv),
varExp(1, ~cv1), # 1 length
#varIdent(form = ~ 1 | cv4), # n partners
varIdent(form = ~ 1 | cv5), # gender
varIdent(form = ~ 1 | cv11) # cohab
),
method="ML")
a4 <- anova(mcv3,mcv4)
aic3[2] <- a4$AIC[2]-a4$AIC[1] # AIC of full minus reduced model
mcv3 <- lme(dv ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
varIdent(form = ~ 1 | iv),
varExp(1, ~cv1), # length
varIdent(form = ~ 1 | cv4), # n partners
#varIdent(form = ~ 1 | cv5), # gender
varIdent(form = ~ 1 | cv11) # cohab
),
method="ML")
a5 <- anova(mcv3,mcv4)
aic3[3] <- a5$AIC[2]-a5$AIC[1] # AIC of full minus reduced model
mcv3 <- lme(dv ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
varIdent(form = ~ 1 | iv),
varExp(1, ~cv1), # length
varIdent(form = ~ 1 | cv4), # n partners
varIdent(form = ~ 1 | cv5) # gender
#varIdent(form = ~ 1 | cv11) # cohab
),
method="ML")
a11 <- anova(mcv3,mcv4)
aic3[4] <- a11$AIC[2]-a11$AIC[1] # AIC of full minus reduced model
print(aic3)
## [1] -7.822819 -12.915617 -4.199764 -2.843710
pvals <- c(pvals, a1$`p-value`, a4$`p-value`, a5$`p-value`, a11$`p-value`)
Remove cohabitation - it improves the model with 3 control variables, but not by much (after adjusting for multiple testing, improvement is not significant).
Model with 3 control variables
mcv3 <- lme(dv ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
varIdent(form = ~ 1 | iv),
varExp(1, ~cv1), # length
varIdent(form = ~ 1 | cv4), # n partners
varIdent(form = ~ 1 | cv5) # gender
),
method="ML")
aic2 <- numeric() # vector that will store AIC values for later comparison
mcv2 <- lme(dv ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
varIdent(form = ~ 1 | iv),
#varExp(1, ~cv1), # length
varIdent(form = ~ 1 | cv4), # n partners
varIdent(form = ~ 1 | cv5) # gender
),
control=lmeControl(opt = "optim"
#opt = "nlminb"
#returnObject=TRUE
),
method="ML")
a1 <- anova(mcv2,mcv3)
aic2[1] <- a1$AIC[2]-a1$AIC[1] # AIC of full minus reduced model
mcv2 <- lme(dv ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
varIdent(form = ~ 1 | iv),
varExp(1, ~cv1), # length
#varIdent(form = ~ 1 | cv4), # n partners
varIdent(form = ~ 1 | cv5) # gender
),
method="ML")
a4 <- anova(mcv2,mcv3)
aic2[2] <- a4$AIC[2]-a4$AIC[1] # AIC of full minus reduced model
mcv2 <- lme(dv ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
varIdent(form = ~ 1 | iv),
varExp(1, ~cv1), # length
varIdent(form = ~ 1 | cv4) # n partners
#varIdent(form = ~ 1 | cv5) # gender
),
method="ML")
a5 <- anova(mcv2,mcv3)
aic2[3] <- a5$AIC[2]-a5$AIC[1] # AIC of full minus reduced model
print(aic2)
## [1] -18.377241 -19.939034 -3.794024
pvals <- c(pvals, a1$`p-value`, a4$`p-value`, a5$`p-value`)
Remove gender - only a small improvement, non-significant after correcting p-value for multiple testing.
aic1 <- numeric()
mcv2 <- lme(dv ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
varIdent(form = ~ 1 | iv),
varExp(1, ~cv1), # length
varIdent(form = ~ 1 | cv4) # n partners
),
method="ML")
mcv1 <- lme(dv ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
varIdent(form = ~ 1 | iv),
#varExp(1, ~cv1), # length
varIdent(form = ~ 1 | cv4) # n partners
),
method="ML")
a1 <- anova(mcv1, mcv2)
aic1[1] <- a1$AIC[2] - a1$AIC[1]
print(a1)
## Model df AIC BIC logLik Test L.Ratio p-value
## mcv1 1 6 1359.641 1384.428 -673.8205
## mcv2 2 7 1339.151 1368.069 -662.5753 1 vs 2 22.49039 <.0001
mcv1 <- lme(dv ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
varIdent(form = ~ 1 | iv),
varExp(1, ~cv1) # length
#varIdent(form = ~ 1 | cv4) # n partners
),
method="ML")
a4 <- anova(mcv1, mcv2)
aic1[2] <- a4$AIC[2] - a4$AIC[1]
print(a4)
## Model df AIC BIC logLik Test L.Ratio p-value
## mcv1 1 5 1367.337 1387.993 -678.6685
## mcv2 2 7 1339.151 1368.069 -662.5753 1 vs 2 32.18648 <.0001
pvals <- c(pvals, a1$`p-value`, a4$`p-value`)
Both variables improve the model with hierarchy type only. That is, relationship length and number of partners account for unexplained variability in the model including hierarchy type.
mcv2 <- lme(dv ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
varIdent(form = ~ 1 | iv),
varExp(1, ~cv1), # length
varIdent(form = ~ 1 | cv4) # n partners
),
method="ML")
miv <- lme(dv ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
#varIdent(form = ~ 1 | iv),
varExp(1, ~cv1), # length
varIdent(form = ~ 1 | cv4) # n partners
),
method="ML")
aiv <- anova(miv, mcv2)
print(aiv)
## Model df AIC BIC logLik Test L.Ratio p-value
## miv 1 6 1348.198 1372.985 -668.0990
## mcv2 2 7 1339.151 1368.069 -662.5753 1 vs 2 11.04751 9e-04
pvals <- c(pvals, aiv$`p-value`)
Likewise, hierarchy type accounts for unexplained variability in the model including the control variables.
Final model for within-participant avoidance:
var.avoidance <- lme(dv ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
varIdent(form = ~ 1 | iv),
varExp(1, ~cv1),
varIdent(form = ~ 1 | cv4)
),
method="REML")
summary(var.avoidance)
## Linear mixed-effects model fit by REML
## Data: NULL
## AIC BIC logLik
## 1343.368 1372.272 -664.6842
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 0.3738429 1.12048
##
## Combination of variance functions:
## Structure: Different standard deviations per stratum
## Formula: ~1 | iv
## Parameter estimates:
## NonHierarchical Hierarchical
## 1.000000 1.298593
## Structure: Exponential of variance covariate
## Formula: ~cv1
## Parameter estimates:
## expon
## -0.2451602
## Structure: Different standard deviations per stratum
## Formula: ~1 | cv4
## Parameter estimates:
## >3 3 2
## 1.0000000 0.8949205 0.6129542
## Fixed effects: dv ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) 1.853254 0.04848688 274 38.22175 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.39799141 -0.63526422 -0.05280815 0.77141679 3.64912019
##
## Number of Observations: 460
## Number of Groups: 186
res <- resid(var.avoidance, type = "pearson")
fit <- fitted(var.avoidance, level=0)
plot(ranef(var.avoidance, level = 1))
plot(res ~ iv)
plot(res ~ cv1)
plot(res ~ cv4)
qqnorm(res)
qqline(res)
Estimated standard deviations for non-hierarchical and hierarchical groups in the final model:
summary(var.avoidance)$sigma # avoidance, non-hierarchical
## [1] 1.12048
summary(var.avoidance)$sigma*coef(var.avoidance$modelStruct$varStruct, uncons=FALSE)
## A.Hierarchical B.expon C.3 C.2
## 1.4550482 -0.2746972 1.0027408 0.6868032
Confidence intervals:
intervals(var.avoidance, level = 0.95, which = "all")
## Approximate 95% confidence intervals
##
## Fixed effects:
## lower est. upper
## (Intercept) 1.757799 1.853254 1.948708
## attr(,"label")
## [1] "Fixed effects:"
##
## Random Effects:
## Level: id
## lower est. upper
## sd((Intercept)) 0.2695546 0.3738429 0.5184794
##
## Variance function:
## lower est. upper
## A.Hierarchical 1.1114660 1.2985932 1.5172253
## B.expon -0.3581005 -0.2451602 -0.1322198
## C.3 0.7277208 0.8949205 1.1005355
## C.2 0.5032379 0.6129542 0.7465910
## attr(,"label")
## [1] "Variance function:"
##
## Within-group standard error:
## lower est. upper
## 0.9434688 1.1204804 1.3307025
Model selection for within-participant variability in anxiety
Check whether adding cohabitation to a model with hierarchy type improves the model
data <- na.omit(study1.long[c(vars,"cohab")])
dv <- data$anxiety
iv <- data$hierarchy
id <- data$Respondent.ID
cv <- data[["cohab"]]
m0 <- lme(dv ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~ 1 | iv),
method="ML")
mcv <- lme(dv ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
varIdent(form = ~ 1 | iv),
varIdent(form = ~ 1 | cv) # cohab
),
method="ML")
a <- anova(m0,mcv)
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m0 1 4 1912.517 1929.858 -952.2586
## mcv 2 6 1904.627 1930.637 -946.3134 1 vs 2 11.89047 0.0026
pvals <- c(pvals, a$`p-value`)
In a model with heterogeneous within-participant variance by hierarchy type, cohabitation significantly accounts for unexplained within-participant variability in anxiety.
data <- na.omit(study1.long[c(vars,"cohab")])
dv <- data$anxiety
iv <- data$hierarchy
id <- data$Respondent.ID
cv <- data[["cohab"]]
m0cv <- lme(dv ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~ 1 | cv),
method="ML")
miv <- lme(dv ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
varIdent(form = ~ 1 | iv),
varIdent(form = ~ 1 | cv) # cohab
),
method="ML")
a <- anova(m0cv,miv)
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m0cv 1 5 1915.071 1936.746 -952.5355
## miv 2 6 1904.627 1930.637 -946.3134 1 vs 2 12.4441 4e-04
pvals <- c(pvals, a$`p-value`)
Final model for within-participant variability in anxiety:
var.anxiety <- lme(dv ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(varIdent(form = ~ 1 | iv),
varIdent(form = ~ 1 | cv)),
method="REML")
summary(var.anxiety)
## Linear mixed-effects model fit by REML
## Data: NULL
## AIC BIC logLik
## 1908.142 1934.142 -948.0712
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 0.7648507 1.061122
##
## Combination of variance functions:
## Structure: Different standard deviations per stratum
## Formula: ~1 | iv
## Parameter estimates:
## NonHierarchical Hierarchical
## 1.000000 1.292722
## Structure: Different standard deviations per stratum
## Formula: ~1 | cv
## Parameter estimates:
## Living Separately, Long Distance Living Separately, local
## 1.0000000 1.0514052
## Living Together
## 0.7558595
## Fixed effects: dv ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) 2.135882 0.06886492 339 31.01553 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.1112234 -0.5300775 -0.2287174 0.5985833 3.4346434
##
## Number of Observations: 564
## Number of Groups: 225
res <- resid(var.anxiety, type = "pearson")
fit <- fitted(var.anxiety, level=0)
plot(ranef(var.anxiety, level = 1))
plot(res ~ iv)
qqnorm(res)
qqline(res)
Model selection for within-participant variability in relationship satisfaction
Check whether adding number of partners to a model with hierarchy type improves de model
data <- na.omit(study1.long[c(vars,"N.Partners.3levels")])
dv <- data$satisfaction
iv <- data$hierarchy
id <- data$Respondent.ID
cv <- data[["N.Partners.3levels"]]
m0 <- lme(dv ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~ 1 | iv), # hierarchy
method="ML")
mcv <- lme(dv ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
varIdent(form = ~ 1 | iv),
varIdent(form = ~ 1 | cv) # number of partners
),
method="ML")
a <- anova(m0,mcv)
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m0 1 4 2969.064 2986.405 -1480.532
## mcv 2 6 2963.177 2989.187 -1475.589 1 vs 2 9.887362 0.0071
pvals <- c(pvals, a$`p-value`)
Number of partners significantly accounts for unexplained within-participant variability in satisfaction.
data <- na.omit(study1.long[c(vars,"N.Partners.3levels")])
dv <- data$satisfaction
iv <- data$hierarchy
id <- data$Respondent.ID
cv <- data[["N.Partners.3levels"]]
m0cv <- lme(dv ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~ 1 | cv), # hierarchy
method="ML")
miv <- lme(dv ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
varIdent(form = ~ 1 | iv),
varIdent(form = ~ 1 | cv) # number of partners
),
method="ML")
a <- anova(m0cv,miv)
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m0cv 1 5 2965.826 2987.501 -1477.913
## miv 2 6 2963.177 2989.187 -1475.589 1 vs 2 4.649 0.0311
pvals <- c(pvals, a$`p-value`)
Final model for within-participant variability in relationship satisfaction:
var.satisfaction <- lme(dv ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(varIdent(form = ~ 1 | iv),
varIdent(form = ~ 1 | cv)),
method="REML")
summary(var.satisfaction)
## Linear mixed-effects model fit by REML
## Data: NULL
## AIC BIC logLik
## 2964.907 2990.907 -1476.454
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 1.753299 3.272212
##
## Combination of variance functions:
## Structure: Different standard deviations per stratum
## Formula: ~1 | iv
## Parameter estimates:
## NonHierarchical Hierarchical
## 1.000000 1.166539
## Structure: Different standard deviations per stratum
## Formula: ~1 | cv
## Parameter estimates:
## >3 3 2
## 1.0000000 0.8775893 0.7517948
## Fixed effects: dv ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) 16.09517 0.1681368 339 95.72663 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.27521043 -0.63054468 0.08644708 0.68576033 1.87229014
##
## Number of Observations: 564
## Number of Groups: 225
res <- resid(var.satisfaction, type = "pearson")
fit <- fitted(var.satisfaction, level=0)
plot(ranef(var.satisfaction, level = 1))
plot(res ~ iv)
qqnorm(res)
qqline(res)
Avoidance, Anxiety, and relationship satisfaction by hierarchy type
Here we use the same variance models from the previous analyses and add fixed effects to it.
Avoidance
Step 1: add hierarchy type to model from previous analysis
data <- na.omit(study1.long[c(vars,"rel.length.c","N.Partners.3levels")])
dv <- data$avoidance
iv <- data$hierarchy
id <- data$Respondent.ID
cv1 <- data[["rel.length.c"]]
cv4 <- data[["N.Partners.3levels"]]
m.avoid <- lme(dv ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
varIdent(form = ~ 1 | iv),
varExp(1, ~cv1), # length
varIdent(form = ~ 1 | cv4) # number of partners
),
method="ML")
m.avoid.h <- lme(dv ~ iv,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
varIdent(form = ~ 1 | iv),
varExp(1, ~cv1), # length
varIdent(form = ~ 1 | cv4) # number of partners
),
method="ML")
a <- anova(m.avoid, m.avoid.h); print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m.avoid 1 7 1631.025 1661.296 -808.5127
## m.avoid.h 2 8 1616.309 1650.904 -800.1547 1 vs 2 16.716 <.0001
summary(m.avoid.h)
## Linear mixed-effects model fit by maximum likelihood
## Data: NULL
## AIC BIC logLik
## 1616.309 1650.904 -800.1547
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 0.3007553 1.066648
##
## Combination of variance functions:
## Structure: Different standard deviations per stratum
## Formula: ~1 | iv
## Parameter estimates:
## NonHierarchical Hierarchical
## 1.000000 1.334169
## Structure: Exponential of variance covariate
## Formula: ~cv1
## Parameter estimates:
## expon
## -0.1589828
## Structure: Different standard deviations per stratum
## Formula: ~1 | cv4
## Parameter estimates:
## >3 3 2
## 1.0000000 0.9096141 0.6829270
## Fixed effects: dv ~ iv
## Value Std.Error DF t-value p-value
## (Intercept) 2.1804693 0.07532642 333 28.946939 0
## ivNonHierarchical -0.3996034 0.09211576 223 -4.338057 0
## Correlation:
## (Intr)
## ivNonHierarchical -0.818
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.52140835 -0.68782754 -0.09703162 0.68923546 4.13118995
##
## Number of Observations: 558
## Number of Groups: 225
pvals <- c(pvals, a$`p-value`)
Adding hierarchy type as a fixed effect improved the variance-only model for avoidance.
Step 2: control variables
Checking which potential control variables are related to avoidance:
for(x in cvs){
if(x != "rel.length.c" & x != "N.Partners.3levels"){
data <- na.omit(study1.long[c(vars,"rel.length.c","N.Partners.3levels",x)])
} else {
data <- na.omit(study1.long[c(vars,"rel.length.c","N.Partners.3levels")])
}
print(x)
dv <- data$avoidance
iv <- data$hierarchy
id <- data$Respondent.ID
cv1 <- data$rel.length.c
cv4 <- data$N.Partners.3levels
cv <- data[[x]]
m1 <- lme(dv ~ 1, random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
varIdent(form = ~ 1 | iv),
varExp(1, ~cv1), # length
varIdent(form = ~ 1 | cv4) # number of partners
),
method="ML")
mcv <- lme(dv ~ cv, random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
varIdent(form = ~ 1 | iv),
varExp(1, ~cv1), # length
varIdent(form = ~ 1 | cv4) # number of partners
),
method="ML")
a1 <- anova(m1,mcv); print(a1)
pvals <- c(pvals, a1$`p-value`)
}
## [1] "rel.length.c"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 7 1631.025 1661.296 -808.5127
## mcv 2 8 1571.818 1606.413 -777.9089 1 vs 2 61.20755 <.0001
## [1] "age.c"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 7 1631.025 1661.296 -808.5127
## mcv 2 8 1630.423 1665.018 -807.2116 1 vs 2 2.602187 0.1067
## [1] "incomeperPPL.c"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 7 1264.281 1292.727 -625.1405
## mcv 2 8 1265.613 1298.123 -624.8065 1 vs 2 0.668069 0.4137
## [1] "N.Partners.3levels"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 7 1631.025 1661.296 -808.5127
## mcv 2 9 1621.103 1660.023 -801.5516 1 vs 2 13.92204 9e-04
## [1] "gender"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 7 1339.151 1368.069 -662.5753
## mcv 2 9 1341.357 1378.538 -661.6783 1 vs 2 1.793924 0.4078
## [1] "orientation.binary"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 7 1327.215 1356.073 -656.6077
## mcv 2 8 1329.213 1362.193 -656.6064 1 vs 2 0.002461744 0.9604
## [1] "race.binary"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 7 1335.926 1364.814 -660.9630
## mcv 2 8 1337.913 1370.928 -660.9564 1 vs 2 0.01326153 0.9083
## [1] "children.binary"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 7 1339.151 1368.069 -662.5753
## mcv 2 8 1335.830 1368.879 -659.9147 1 vs 2 5.321132 0.0211
## [1] "education.3levels"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 7 1318.133 1346.944 -652.0666
## mcv 2 9 1319.947 1356.991 -650.9737 1 vs 2 2.185765 0.3352
## [1] "marital.status"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 7 1631.025 1661.296 -808.5127
## mcv 2 8 1574.408 1609.003 -779.2040 1 vs 2 58.61739 <.0001
## [1] "cohab"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 7 1631.025 1661.296 -808.5127
## mcv 2 9 1571.167 1610.086 -776.5835 1 vs 2 63.8583 <.0001
## [1] "coparent"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 7 1631.025 1661.296 -808.5127
## mcv 2 9 1605.653 1644.572 -793.8263 1 vs 2 29.3727 <.0001
The following variables were related to avoidance: relationship length, number of partners, marital status, cohabitation, and coparenting. PS: After adjusting for multiple testing, children (p=0.02) is not a significant predictor.
Model selection strategy: start with a model that includes hierarchy (iv of interest), relationship length, number of partners, marital status, cohabitation status, and coparenting status, all as fixed effects variables. Then perform model selection through a backwards elimination process, that is, only variables that significantly improve a model with fewer variables are left in the final model.
cvs.avoid.fix <- c("rel.length.c", "N.Partners.3levels",
"marital.status", "cohab", "coparent")
data <- na.omit(study1.long[c(vars,cvs.avoid.fix)])
dv <- data$avoidance
iv <- data$hierarchy
id <- data$Respondent.ID
cv1 <- data$rel.length.c
cv4 <- data$N.Partners.3levels
cv10 <- data$marital.status
cv11 <- data$cohab
cv12 <- data$coparent
aic0 <- numeric()
m0 <- lme(dv ~ iv + cv1 + cv4 + cv10 + cv11 + cv12,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
varIdent(form = ~ 1 | iv),
varExp(1, ~cv1), # length
varIdent(form = ~ 1 | cv4) # number of partners
),
method="ML")
print(cvs.avoid.fix[1])
## [1] "rel.length.c"
m <- update(m0, .~. -cv1)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[1] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 15 1542.469 1607.334 -756.2343
## m0 2 16 1534.467 1603.656 -751.2333 1 vs 2 10.00203 0.0016
pvals <- c(pvals, a$`p-value`)
print(cvs.avoid.fix[2])
## [1] "N.Partners.3levels"
m <- update(m0, .~. -cv4)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[2] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 14 1534.653 1595.194 -753.3266
## m0 2 16 1534.467 1603.656 -751.2333 1 vs 2 4.186693 0.1233
pvals <- c(pvals, a$`p-value`)
print(cvs.avoid.fix[3])
## [1] "marital.status"
m <- update(m0, .~. -cv10)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[3] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 15 1533.913 1598.779 -751.9567
## m0 2 16 1534.467 1603.656 -751.2333 1 vs 2 1.446789 0.229
pvals <- c(pvals, a$`p-value`)
print(cvs.avoid.fix[4])
## [1] "cohab"
m <- update(m0, .~. -cv11)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[4] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 14 1540.894 1601.436 -756.4472
## m0 2 16 1534.467 1603.656 -751.2333 1 vs 2 10.42784 0.0054
pvals <- c(pvals, a$`p-value`)
print(cvs.avoid.fix[5])
## [1] "coparent"
m <- update(m0, .~. -cv12)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[5] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 14 1530.730 1591.271 -751.3651
## m0 2 16 1534.467 1603.656 -751.2333 1 vs 2 0.2637 0.8765
pvals <- c(pvals, a$`p-value`)
print(aic0)
## [1] -8.0020353 -0.1866934 0.5532115 -6.4278388 3.7363000
Remove coparenting - it does not improve model with 4 control variables and increases AIC the most.
cvs.avoid.fix <- c("rel.length.c", "N.Partners.3levels",
"marital.status", "cohab")
aic0 <- numeric()
data <- na.omit(study1.long[c(vars,cvs.avoid.fix)])
dv <- data$avoidance
iv <- data$hierarchy
id <- data$Respondent.ID
cv1 <- data$rel.length.c
cv4 <- data$N.Partners.3levels
cv10 <- data$marital.status
cv11 <- data$cohab
aic1 <- numeric()
m0 <- lme(dv ~ iv + cv1 + cv4 + cv10 + cv11,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
varIdent(form = ~ 1 | iv),
varExp(1, ~cv1), # length
varIdent(form = ~ 1 | cv4) # number of partners
),
method="ML")
print(cvs.avoid.fix[1])
## [1] "rel.length.c"
m <- update(m0, .~. -cv1)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[1] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 13 1540.177 1596.394 -757.0886
## m0 2 14 1530.730 1591.271 -751.3651 1 vs 2 11.44684 7e-04
pvals <- c(pvals, a$`p-value`)
print(cvs.avoid.fix[2])
## [1] "N.Partners.3levels"
m <- update(m0, .~. -cv4)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[2] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 12 1530.774 1582.667 -753.3873
## m0 2 14 1530.730 1591.271 -751.3651 1 vs 2 4.044229 0.1324
pvals <- c(pvals, a$`p-value`)
print(cvs.avoid.fix[3])
## [1] "marital.status"
m <- update(m0, .~. -cv10)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[3] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 13 1530.151 1586.368 -752.0755
## m0 2 14 1530.730 1591.271 -751.3651 1 vs 2 1.420667 0.2333
pvals <- c(pvals, a$`p-value`)
print(cvs.avoid.fix[4])
## [1] "cohab"
m <- update(m0, .~. -cv11)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[4] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 12 1537.077 1588.970 -756.5386
## m0 2 14 1530.730 1591.271 -751.3651 1 vs 2 10.34701 0.0057
pvals <- c(pvals, a$`p-value`)
print(aic0)
## [1] -9.44684043 -0.04422927 0.57933324 -6.34700617
Remove marital status - it does not improve model with 3 control variables and increases AIC the most.
cvs.avoid.fix <- c("rel.length.c", "N.Partners.3levels", "cohab")
aic0 <- numeric()
data <- na.omit(study1.long[c(vars,cvs.avoid.fix)])
dv <- data$avoidance
iv <- data$hierarchy
id <- data$Respondent.ID
cv1 <- data$rel.length.c
cv4 <- data$N.Partners.3levels
cv11 <- data$cohab
aic1 <- numeric()
m0 <- lme(dv ~ iv + cv1 + cv4 + cv11,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
varIdent(form = ~ 1 | iv),
varExp(1, ~cv1), # length
varIdent(form = ~ 1 | cv4) # number of partners
),
method="ML")
print(cvs.avoid.fix[1])
## [1] "rel.length.c"
m <- update(m0, .~. -cv1)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[1] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 12 1554.472 1606.364 -765.2361
## m0 2 13 1530.151 1586.368 -752.0755 1 vs 2 26.32116 <.0001
pvals <- c(pvals, a$`p-value`)
print(cvs.avoid.fix[2])
## [1] "N.Partners.3levels"
m <- update(m0, .~. -cv4)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[2] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 11 1530.866 1578.434 -754.4328
## m0 2 13 1530.151 1586.368 -752.0755 1 vs 2 4.714721 0.0947
pvals <- c(pvals, a$`p-value`)
print(cvs.avoid.fix[3])
## [1] "cohab"
m <- update(m0, .~. -cv11)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[3] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 11 1542.581 1590.148 -760.2902
## m0 2 13 1530.151 1586.368 -752.0755 1 vs 2 16.42951 3e-04
pvals <- c(pvals, a$`p-value`)
print(aic0)
## [1] -24.321164 -0.714721 -12.429509
Remove number of partners - it has a small, non-significant improvement in the model.
cvs.avoid.fix <- c("rel.length.c", "cohab")
aic0 <- numeric()
data <- na.omit(study1.long[c(vars,cvs.avoid.fix)])
dv <- data$avoidance
iv <- data$hierarchy
id <- data$Respondent.ID
cv1 <- data$rel.length.c
cv11 <- data$cohab
aic1 <- numeric()
m0 <- lme(dv ~ iv + cv1 + cv11,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
varIdent(form = ~ 1 | iv),
varExp(1, ~cv1), # length
varIdent(form = ~ 1 | cv4) # number of partners
),
method="ML")
print(cvs.avoid.fix[1])
## [1] "rel.length.c"
m <- update(m0, .~. -cv1)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[1] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 10 1554.021 1597.265 -767.0105
## m0 2 11 1530.866 1578.434 -754.4328 1 vs 2 25.15542 <.0001
pvals <- c(pvals, a$`p-value`)
print(cvs.avoid.fix[2])
## [1] "cohab"
m <- update(m0, .~. -cv11)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[2] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 9 1547.867 1586.786 -764.9335
## m0 2 11 1530.866 1578.434 -754.4328 1 vs 2 21.00124 <.0001
pvals <- c(pvals, a$`p-value`)
print(aic0)
## [1] -23.15542 -17.00124
Keep length of relationship and cohabitation status - both significantly improve the model
cvs.avoid.fix <- c("rel.length.c", "cohab")
aic0 <- numeric()
data <- na.omit(study1.long[c(vars,cvs.avoid.fix)])
dv <- data$avoidance
iv <- data$hierarchy
id <- data$Respondent.ID
cv1 <- data$rel.length.c
cv11 <- data$cohab
aic1 <- numeric()
m0 <- lme(dv ~ iv + cv1 + cv11,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
varIdent(form = ~ 1 | iv),
varExp(1, ~cv1), # length
varIdent(form = ~ 1 | cv4) # number of partners
),
method="ML")
print(cvs.avoid.fix[1])
## [1] "rel.length.c"
m <- update(m0, .~. -iv)
a <- anova(m,m0)
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 10 1554.332 1597.576 -767.1660
## m0 2 11 1530.866 1578.434 -754.4328 1 vs 2 25.46638 <.0001
pvals <- c(pvals, a$`p-value`)
Hierarchy type improves a model with the CVs.
Residual plots:
summary(m0)
## Linear mixed-effects model fit by maximum likelihood
## Data: NULL
## AIC BIC logLik
## 1530.866 1578.434 -754.4328
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 0.2857844 0.9621063
##
## Combination of variance functions:
## Structure: Different standard deviations per stratum
## Formula: ~1 | iv
## Parameter estimates:
## NonHierarchical Hierarchical
## 1.000000 1.232694
## Structure: Exponential of variance covariate
## Formula: ~cv1
## Parameter estimates:
## expon
## -0.2443413
## Structure: Different standard deviations per stratum
## Formula: ~1 | cv4
## Parameter estimates:
## >3 3 2
## 1.0000000 0.9252155 0.7461111
## Fixed effects: dv ~ iv + cv1 + cv11
## Value Std.Error DF t-value
## (Intercept) 2.4571677 0.08247076 330 29.794411
## ivNonHierarchical -0.4326165 0.08266449 223 -5.233402
## cv1 -0.1579385 0.03066157 330 -5.151024
## cv11Living Separately, Long Distance -0.0129656 0.11883148 330 -0.109109
## cv11Living Together -0.4098462 0.09478961 330 -4.323746
## p-value
## (Intercept) 0.0000
## ivNonHierarchical 0.0000
## cv1 0.0000
## cv11Living Separately, Long Distance 0.9132
## cv11Living Together 0.0000
## Correlation:
## (Intr) ivNnHr cv1 c11SLD
## ivNonHierarchical -0.620
## cv1 0.064 0.125
## cv11Living Separately, Long Distance -0.405 -0.022 -0.204
## cv11Living Together -0.529 -0.047 -0.552 0.445
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.5885507 -0.7113122 -0.1979085 0.5684950 3.9868007
##
## Number of Observations: 558
## Number of Groups: 225
res <- resid(m0, type = "pearson")
fit <- fitted(m0, level=0)
plot(ranef(m0, level = 1))
plot(res ~ iv)
plot(res ~ cv1)
plot(res ~ cv11)
qqnorm(res)
qqline(res)
Differences in avoidance between non-hierarchical relationships and 3 levels of hierarchy (primary, secondary, and tertiary partners)
Checking whether avoidance in non-hierarchical relationships differ from those in primary, secondary, and tertiary relationships:
data <- na.omit(study1.long[c(vars,"hierarchy2","N.Partners.3levels", "rel.length.c","cohab")])
dv <- data$avoidance
iv <- data$hierarchy # 2 levels: hierarchical & non-hierarchical
iv2 <- data$hierarchy2 # 4 leveles: non-hierarchical, primary, secondary, and tertiary
id <- data$Respondent.ID
cv1 <- data$rel.length.c
cv4 <- data$N.Partners.3levels
cv11 <- data$cohab
m.avoid <- lme(dv ~ 1,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
varIdent(form = ~ 1 | iv),
varExp(1, ~cv1), # length
varIdent(form = ~ 1 | cv4) # number of partners
),
method="ML")
m.avoid.h <- lme(dv ~ iv2,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
varIdent(form = ~ 1 | iv),
varExp(1, ~cv1), # length
varIdent(form = ~ 1 | cv4) # number of partners
),
method="ML")
m.avoid.hcv <- lme(dv ~ iv2 + cv1 + cv11,
random = list(id = pdSymm(form = ~ 1)),
weights = varComb(
varIdent(form = ~ 1 | iv),
varExp(1, ~cv1), # length
varIdent(form = ~ 1 | cv4) # number of partners
),
method="ML")
print(anova(m.avoid, m.avoid.h, m.avoid.hcv))
## Model df AIC BIC logLik Test L.Ratio p-value
## m.avoid 1 7 1631.025 1661.296 -808.5127
## m.avoid.h 2 11 1517.805 1565.373 -747.9023 1 vs 2 121.22067 <.0001
## m.avoid.hcv 3 14 1487.905 1548.446 -729.9523 2 vs 3 35.90009 <.0001
summary(m.avoid)
## Linear mixed-effects model fit by maximum likelihood
## Data: NULL
## AIC BIC logLik
## 1631.025 1661.296 -808.5127
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 0.3729916 1.044518
##
## Combination of variance functions:
## Structure: Different standard deviations per stratum
## Formula: ~1 | iv
## Parameter estimates:
## NonHierarchical Hierarchical
## 1.000000 1.387752
## Structure: Exponential of variance covariate
## Formula: ~cv1
## Parameter estimates:
## expon
## -0.2166751
## Structure: Different standard deviations per stratum
## Formula: ~1 | cv4
## Parameter estimates:
## >3 3 2
## 1.0000000 0.9022662 0.6753020
## Fixed effects: dv ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) 1.884467 0.04491353 333 41.95767 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.52312184 -0.64435785 -0.06007659 0.75089252 3.88281866
##
## Number of Observations: 558
## Number of Groups: 225
summary(m.avoid.h)
## Linear mixed-effects model fit by maximum likelihood
## Data: NULL
## AIC BIC logLik
## 1517.805 1565.373 -747.9023
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 0.3824019 1.051132
##
## Combination of variance functions:
## Structure: Different standard deviations per stratum
## Formula: ~1 | iv
## Parameter estimates:
## NonHierarchical Hierarchical
## 1.00000 1.05426
## Structure: Exponential of variance covariate
## Formula: ~cv1
## Parameter estimates:
## expon
## -0.2173105
## Structure: Different standard deviations per stratum
## Formula: ~1 | cv4
## Parameter estimates:
## >3 3 2
## 1.0000000 0.8861006 0.6673136
## Fixed effects: dv ~ iv2
## Value Std.Error DF t-value p-value
## (Intercept) 1.7633776 0.05569555 329 31.66101 0.0000
## iv2Other 1.0512247 0.26335131 329 3.99172 0.0001
## iv2Primary -0.1287850 0.09598497 329 -1.34172 0.1806
## iv2Secondary 0.9270025 0.11643425 329 7.96160 0.0000
## iv2Tertiary 1.7063925 0.22573872 329 7.55915 0.0000
## Correlation:
## (Intr) iv2Oth iv2Prm iv2Scn
## iv2Other -0.211
## iv2Primary -0.580 0.168
## iv2Secondary -0.478 0.127 0.421
## iv2Tertiary -0.247 0.060 0.195 0.152
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.0144395 -0.6164220 -0.1136943 0.5695850 3.9809269
##
## Number of Observations: 558
## Number of Groups: 225
summary(m.avoid.hcv)
## Linear mixed-effects model fit by maximum likelihood
## Data: NULL
## AIC BIC logLik
## 1487.905 1548.446 -729.9523
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 0.3309848 0.9593045
##
## Combination of variance functions:
## Structure: Different standard deviations per stratum
## Formula: ~1 | iv
## Parameter estimates:
## NonHierarchical Hierarchical
## 1.000000 1.116213
## Structure: Exponential of variance covariate
## Formula: ~cv1
## Parameter estimates:
## expon
## -0.2523947
## Structure: Different standard deviations per stratum
## Formula: ~1 | cv4
## Parameter estimates:
## >3 3 2
## 1.0000000 0.9200493 0.7185425
## Fixed effects: dv ~ iv2 + cv1 + cv11
## Value Std.Error DF t-value
## (Intercept) 1.9126557 0.07416250 326 25.790066
## iv2Other 0.9321925 0.25859749 326 3.604801
## iv2Primary 0.0856490 0.09602067 326 0.891985
## iv2Secondary 0.7589260 0.12083781 326 6.280534
## iv2Tertiary 1.5239990 0.22640467 326 6.731306
## cv1 -0.1177103 0.03069555 326 -3.834768
## cv11Living Separately, Long Distance 0.0276351 0.11540364 326 0.239465
## cv11Living Together -0.2368111 0.09632493 326 -2.458461
## p-value
## (Intercept) 0.0000
## iv2Other 0.0004
## iv2Primary 0.3731
## iv2Secondary 0.0000
## iv2Tertiary 0.0000
## cv1 0.0002
## cv11Living Separately, Long Distance 0.8109
## cv11Living Together 0.0145
## Correlation:
## (Intr) iv2Oth iv2Prm iv2Scn iv2Trt
## iv2Other -0.178
## iv2Primary -0.244 0.123
## iv2Secondary -0.472 0.124 0.248
## iv2Tertiary -0.261 0.061 0.118 0.171
## cv1 0.162 0.006 -0.204 0.044 0.003
## cv11Living Separately, Long Distance -0.475 0.014 -0.013 0.058 0.045
## cv11Living Together -0.671 0.053 -0.140 0.234 0.140
## cv1 c11SLD
## iv2Other
## iv2Primary
## iv2Secondary
## iv2Tertiary
## cv1
## cv11Living Separately, Long Distance -0.200
## cv11Living Together -0.455 0.444
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.9578732 -0.6828616 -0.1928818 0.5702773 4.0908065
##
## Number of Observations: 558
## Number of Groups: 225
After controling for relationship length and cohabitation, avoidance in non-hierarchical relationships is lower, on average, than avoidance in secondary and tertiary relationships, but is at a similar level as primary ones.
Anxiety
Step 1: add hierarchy type to model from previous analysis
data <- na.omit(study1.long[vars])
dv <- data$anxiety
iv <- data$hierarchy
id <- data$Respondent.ID
m.anx <- lme(dv ~ 1, random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~ 1 | iv),
method="ML")
m.anx.h <- lme(dv ~ iv, random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~ 1 | iv),
method="ML")
summary(m.anx.h)
## Linear mixed-effects model fit by maximum likelihood
## Data: NULL
## AIC BIC logLik
## 1900.971 1922.647 -945.4857
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 0.7270242 1.012684
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | iv
## Parameter estimates:
## NonHierarchical Hierarchical
## 1.000000 1.258141
## Fixed effects: dv ~ iv
## Value Std.Error DF t-value p-value
## (Intercept) 2.5509474 0.1104378 339 23.098493 0e+00
## ivNonHierarchical -0.5251889 0.1406963 223 -3.732784 2e-04
## Correlation:
## (Intr)
## ivNonHierarchical -0.785
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.3699878 -0.5642526 -0.3003285 0.5140620 3.5548607
##
## Number of Observations: 564
## Number of Groups: 225
a <- anova(m.anx, m.anx.h); print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m.anx 1 4 1912.517 1929.858 -952.2586
## m.anx.h 2 5 1900.972 1922.647 -945.4857 1 vs 2 13.54583 2e-04
pvals <- c(pvals, a$`p-value`)
Adding hierarchy type as a fixed effect improved the variance model (no fixed effects) for anxiety.
Step 2: control variables
Checking which potential control variables are related to anxiety:
for(x in cvs){
data <- na.omit(study1.long[c(vars,x)])
print(x)
dv <- data$anxiety
iv <- data$hierarchy
id <- data$Respondent.ID
cv <- data[[x]]
m1 <- lme(dv ~ 1, random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~ 1 | iv),
method="ML")
mcv <- lme(dv ~ cv, random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~ 1 | iv),
method="ML")
a <- anova(m1,mcv)
print(a)
pvals <- c(pvals, a$`p-value`)
}
## [1] "rel.length.c"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 4 1889.623 1906.921 -940.8116
## mcv 2 5 1853.985 1875.606 -921.9923 1 vs 2 37.63856 <.0001
## [1] "age.c"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 4 1912.517 1929.858 -952.2586
## mcv 2 5 1913.950 1935.626 -951.9751 1 vs 2 0.5670111 0.4514
## [1] "incomeperPPL.c"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 4 1488.719 1505.030 -740.3597
## mcv 2 5 1490.355 1510.744 -740.1777 1 vs 2 0.364047 0.5463
## [1] "N.Partners.3levels"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 4 1912.517 1929.858 -952.2586
## mcv 2 6 1914.753 1940.763 -951.3764 1 vs 2 1.764562 0.4138
## [1] "gender"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 4 1581.768 1598.345 -786.8839
## mcv 2 6 1583.705 1608.570 -785.8526 1 vs 2 2.062734 0.3565
## [1] "orientation.binary"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 4 1558.268 1574.802 -775.1341
## mcv 2 5 1556.183 1576.850 -773.0913 1 vs 2 4.085704 0.0432
## [1] "race.binary"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 4 1569.976 1586.536 -780.9881
## mcv 2 5 1571.976 1592.675 -780.9880 1 vs 2 0.0002802241 0.9866
## [1] "children.binary"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 4 1581.768 1598.345 -786.8839
## mcv 2 5 1583.352 1604.073 -786.6760 1 vs 2 0.4157617 0.5191
## [1] "education.3levels"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 4 1560.868 1577.384 -776.4339
## mcv 2 6 1563.118 1587.892 -775.5587 1 vs 2 1.750394 0.4168
## [1] "marital.status"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 4 1912.517 1929.858 -952.2586
## mcv 2 5 1884.621 1906.296 -937.3104 1 vs 2 29.89654 <.0001
## [1] "cohab"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 4 1912.517 1929.858 -952.2586
## mcv 2 6 1873.617 1899.628 -930.8086 1 vs 2 42.90002 <.0001
## [1] "coparent"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 4 1912.517 1929.858 -952.2586
## mcv 2 6 1905.634 1931.644 -946.8169 1 vs 2 10.88358 0.0043
The following variables were related to anxiety: relationship length, sexual orientation, marital status, cohabitation, co-parenting.
Model selection strategy: start with a model that includes hierarchy (iv of interest), relationship length, sexual orientation, marital status, cohabitation, co-parenting, all as fixed effects variables. Then perform model selection through a backwards elimination process.
cvs.anx.fix <- c("rel.length.c", "orientation", "marital.status", "cohab", "coparent")
data <- na.omit(study1.long[c(vars,cvs.anx.fix)])
dv <- data$anxiety
iv <- data$hierarchy
id <- data$Respondent.ID
cv1 <- data$rel.length.c
cv6 <- data$orientation
cv10 <- data$marital.status
cv11 <- data$cohab
cv12 <- data$coparent
aic0 <- numeric()
m0 <- lme(dv ~ iv + cv1 + cv6 + cv10 + cv11 + cv12,
random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~ 1 | iv),
method="ML")
print(cvs.anx.fix[1])
## [1] "rel.length.c"
m <- update(m0, .~. -cv1)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[1] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 17 1521.022 1591.253 -743.5111
## m0 2 18 1517.847 1592.209 -740.9235 1 vs 2 5.175343 0.0229
pvals <- c(pvals, a$`p-value`)
print(cvs.anx.fix[2])
## [1] "orientation"
m <- update(m0, .~. -cv6)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[2] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 11 1513.181 1558.624 -745.5903
## m0 2 18 1517.847 1592.209 -740.9235 1 vs 2 9.333689 0.2296
pvals <- c(pvals, a$`p-value`)
print(cvs.anx.fix[3])
## [1] "marital.status"
m <- update(m0, .~. -cv10)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[3] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 17 1516.114 1586.345 -741.0571
## m0 2 18 1517.847 1592.209 -740.9235 1 vs 2 0.2672973 0.6052
pvals <- c(pvals, a$`p-value`)
print(cvs.anx.fix[4])
## [1] "cohab"
m <- update(m0, .~. -cv11)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[4] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 16 1526.345 1592.444 -747.1723
## m0 2 18 1517.847 1592.209 -740.9235 1 vs 2 12.49774 0.0019
pvals <- c(pvals, a$`p-value`)
print(cvs.anx.fix[5])
## [1] "coparent"
m <- update(m0, .~. -cv12)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[5] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 16 1515.306 1581.406 -741.6530
## m0 2 18 1517.847 1592.209 -740.9235 1 vs 2 1.459097 0.4821
pvals <- c(pvals, a$`p-value`)
print(aic0)
## [1] -3.175343 4.666311 1.732703 -8.497738 2.540903
Remove coparenting - it does not improve model with 4 control variables and increases AIC the most.
cvs.anx.fix <- c("rel.length.c", "orientation", "marital.status", "cohab")
data <- na.omit(study1.long[c(vars,cvs.anx.fix)])
dv <- data$anxiety
iv <- data$hierarchy
id <- data$Respondent.ID
cv1 <- data$rel.length.c
cv6 <- data$orientation
cv10 <- data$marital.status
cv11 <- data$cohab
aic0 <- numeric()
m0 <- lme(dv ~ iv + cv1 + cv6 + cv10 + cv11,
random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~ 1 | iv),
method="ML")
print(cvs.anx.fix[1])
## [1] "rel.length.c"
m <- update(m0, .~. -cv1)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[1] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 15 1517.321 1579.289 -743.6602
## m0 2 16 1515.306 1581.406 -741.6530 1 vs 2 4.014495 0.0451
pvals <- c(pvals, a$`p-value`)
print(cvs.anx.fix[2])
## [1] "orientation"
m <- update(m0, .~. -cv6)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[2] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 9 1510.106 1547.287 -746.053
## m0 2 16 1515.306 1581.406 -741.653 1 vs 2 8.800016 0.2673
pvals <- c(pvals, a$`p-value`)
print(cvs.anx.fix[3])
## [1] "marital.status"
m <- update(m0, .~. -cv10)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[3] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 15 1513.449 1575.418 -741.7248
## m0 2 16 1515.306 1581.406 -741.6530 1 vs 2 0.1435067 0.7048
pvals <- c(pvals, a$`p-value`)
print(cvs.anx.fix[4])
## [1] "cohab"
m <- update(m0, .~. -cv11)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[4] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 14 1523.208 1581.045 -747.6041
## m0 2 16 1515.306 1581.406 -741.6530 1 vs 2 11.90224 0.0026
pvals <- c(pvals, a$`p-value`)
print(aic0)
## [1] -2.014495 5.199984 1.856493 -7.902243
Remove orientation - it does not improve a model with 3 control variables and it increases AIC the most.
cvs.anx.fix <- c("rel.length.c", "marital.status", "cohab")
data <- na.omit(study1.long[c(vars,cvs.anx.fix)])
dv <- data$anxiety
iv <- data$hierarchy
id <- data$Respondent.ID
cv1 <- data$rel.length.c
cv10 <- data$marital.status
cv11 <- data$cohab
aic0 <- numeric()
m0 <- lme(dv ~ iv + cv1 + cv10 + cv11,
random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~ 1 | iv),
method="ML")
print(cvs.anx.fix[1])
## [1] "rel.length.c"
m <- update(m0, .~. -cv1)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[1] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 8 1835.759 1870.354 -909.8798
## m0 2 9 1831.461 1870.380 -906.7305 1 vs 2 6.298463 0.0121
pvals <- c(pvals, a$`p-value`)
print(cvs.anx.fix[2])
## [1] "marital.status"
m <- update(m0, .~. -cv10)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[2] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 8 1829.702 1864.297 -906.8510
## m0 2 9 1831.461 1870.380 -906.7305 1 vs 2 0.2410078 0.6235
pvals <- c(pvals, a$`p-value`)
print(cvs.anx.fix[3])
## [1] "cohab"
m <- update(m0, .~. -cv11)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[3] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 7 1838.741 1869.012 -912.3707
## m0 2 9 1831.461 1870.380 -906.7305 1 vs 2 11.28028 0.0036
pvals <- c(pvals, a$`p-value`)
print(aic0)
## [1] -4.298463 1.758992 -7.280277
Remove marital status.
cvs.anx.fix <- c("rel.length.c", "cohab")
data <- na.omit(study1.long[c(vars,cvs.anx.fix)])
dv <- data$anxiety
iv <- data$hierarchy
id <- data$Respondent.ID
cv1 <- data$rel.length.c
cv11 <- data$cohab
aic0 <- numeric()
m0 <- lme(dv ~ iv + cv1 + cv11,
random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~ 1 | iv),
method="ML")
print(cvs.anx.fix[1])
## [1] "rel.length.c"
m <- update(m0, .~. -cv1)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[1] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 7 1840.053 1870.323 -913.0264
## m0 2 8 1829.702 1864.297 -906.8510 1 vs 2 12.35068 4e-04
pvals <- c(pvals, a$`p-value`)
print(cvs.anx.fix[2])
## [1] "cohab"
m <- update(m0, .~. -cv11)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[2] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 6 1840.161 1866.108 -914.0808
## m0 2 8 1829.702 1864.297 -906.8510 1 vs 2 14.45945 7e-04
pvals <- c(pvals, a$`p-value`)
print(aic0)
## [1] -10.35068 -10.45945
Keep relationship length and cohabitation.
summary(m0)
## Linear mixed-effects model fit by maximum likelihood
## Data: NULL
## AIC BIC logLik
## 1829.702 1864.297 -906.851
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 0.7597411 0.9573307
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | iv
## Parameter estimates:
## NonHierarchical Hierarchical
## 1.000000 1.202774
## Fixed effects: dv ~ iv + cv1 + cv11
## Value Std.Error DF t-value
## (Intercept) 2.7506465 0.12221267 330 22.507049
## ivNonHierarchical -0.5434799 0.13959761 223 -3.893189
## cv1 -0.2122172 0.06029459 330 -3.519673
## cv11Living Separately, Long Distance -0.0406754 0.15228300 330 -0.267104
## cv11Living Together -0.4729560 0.13040235 330 -3.626898
## p-value
## (Intercept) 0.0000
## ivNonHierarchical 0.0001
## cv1 0.0005
## cv11Living Separately, Long Distance 0.7896
## cv11Living Together 0.0003
## Correlation:
## (Intr) ivNnHr cv1 c11SLD
## ivNonHierarchical -0.661
## cv1 0.198 0.075
## cv11Living Separately, Long Distance -0.316 -0.037 -0.182
## cv11Living Together -0.437 -0.050 -0.548 0.381
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.9022303 -0.5992908 -0.2149004 0.4898987 3.4076844
##
## Number of Observations: 558
## Number of Groups: 225
Residual plots:
res <- resid(m0, type = "pearson")
fit <- fitted(m0, level=0)
plot(ranef(m0, level = 1))
plot(res ~ iv)
plot(res ~ cv1) # rel. length
plot(res ~ cv11) # cohab
qqnorm(res)
qqline(res)
Differences in anxiety between non-hierarchical relationships and 3 levels of hierarchy (primary, secondary, and tertiary partners)
data <- na.omit(study1.long[c(vars,"hierarchy2","rel.length.c","cohab")])
dv <- data$anxiety
iv <- data$hierarchy
iv2 <- data$hierarchy2
id <- data$Respondent.ID
cv1 <- data$rel.length.c
cv11 <- data$cohab
m.anx <- lme(dv ~ 1, random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~ 1 | iv),
method="ML")
m.anx.h <- lme(dv ~ iv2, random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~ 1 | iv),
method="ML")
m.anx.hcv <- lme(dv ~ iv2 + cv1 + cv11, random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~ 1 | iv),
method="ML")
print(anova(m.anx, m.anx.h, m.anx.hcv))
## Model df AIC BIC logLik Test L.Ratio p-value
## m.anx 1 4 1889.623 1906.921 -940.8116
## m.anx.h 2 8 1851.256 1885.851 -917.6280 1 vs 2 46.36725 <.0001
## m.anx.hcv 3 11 1829.198 1876.766 -903.5989 2 vs 3 28.05813 <.0001
summary(m.anx.h)
## Linear mixed-effects model fit by maximum likelihood
## Data: NULL
## AIC BIC logLik
## 1851.256 1885.851 -917.628
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 0.7747196 1.003007
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | iv
## Parameter estimates:
## NonHierarchical Hierarchical
## 1.000000 1.132424
## Fixed effects: dv ~ iv2
## Value Std.Error DF t-value p-value
## (Intercept) 2.0335325 0.0904697 329 22.477502 0.0000
## iv2Other 0.6532175 0.3273944 329 1.995201 0.0468
## iv2Primary -0.0167979 0.1660676 329 -0.101151 0.9195
## iv2Secondary 0.9280205 0.1686191 329 5.503650 0.0000
## iv2Tertiary 0.9570468 0.2864278 329 3.341319 0.0009
## Correlation:
## (Intr) iv2Oth iv2Prm iv2Scn
## iv2Other -0.276
## iv2Primary -0.545 0.268
## iv2Secondary -0.537 0.232 0.513
## iv2Tertiary -0.316 0.123 0.294 0.272
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.4489913 -0.5054803 -0.2632330 0.4339564 3.5330230
##
## Number of Observations: 558
## Number of Groups: 225
summary(m.anx.hcv)
## Linear mixed-effects model fit by maximum likelihood
## Data: NULL
## AIC BIC logLik
## 1829.198 1876.766 -903.5989
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 0.7692255 0.9550746
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | iv
## Parameter estimates:
## NonHierarchical Hierarchical
## 1.000000 1.179662
## Fixed effects: dv ~ iv2 + cv1 + cv11
## Value Std.Error DF t-value
## (Intercept) 2.1721841 0.1103088 326 19.691853
## iv2Other 0.5374539 0.3260970 326 1.648141
## iv2Primary 0.2853026 0.1740617 326 1.639089
## iv2Secondary 0.7382674 0.1711560 326 4.313418
## iv2Tertiary 0.7899550 0.2867508 326 2.754849
## cv1 -0.1752501 0.0618196 326 -2.834862
## cv11Living Separately, Long Distance -0.0371184 0.1519449 326 -0.244289
## cv11Living Together -0.3763038 0.1356738 326 -2.773593
## p-value
## (Intercept) 0.0000
## iv2Other 0.1003
## iv2Primary 0.1022
## iv2Secondary 0.0000
## iv2Tertiary 0.0062
## cv1 0.0049
## cv11Living Separately, Long Distance 0.8072
## cv11Living Together 0.0059
## Correlation:
## (Intr) iv2Oth iv2Prm iv2Scn iv2Trt
## iv2Other -0.243
## iv2Primary -0.360 0.229
## iv2Secondary -0.496 0.238 0.405
## iv2Tertiary -0.307 0.128 0.240 0.286
## cv1 0.276 0.017 -0.202 0.048 0.015
## cv11Living Separately, Long Distance -0.395 0.006 0.030 0.019 0.052
## cv11Living Together -0.559 0.052 -0.132 0.167 0.104
## cv1 c11SLD
## iv2Other
## iv2Primary
## iv2Secondary
## iv2Tertiary
## cv1
## cv11Living Separately, Long Distance -0.178
## cv11Living Together -0.442 0.364
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.8626435 -0.5918645 -0.2003911 0.4431720 3.4584526
##
## Number of Observations: 558
## Number of Groups: 225
After controling for relationship length and cohabitation, anxiety in non-hierarchical relationships is lower, on average, than anxiety in secondary and tertiary relationships, but is at a similar level as primary ones.
Satisfaction
Step 1: add hierarchy type to model from previous analysis
data <- na.omit(study1.long[vars])
dv <- data$satisfaction
iv <- data$hierarchy
id <- data$Respondent.ID
m.sat <- lme(dv ~ 1, random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~ 1 | iv),
method="ML")
m.sat.h <- lme(dv ~ iv, random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~ 1 | iv),
method="ML")
a <- anova(m.sat, m.sat.h); print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m.sat 1 4 2969.064 2986.405 -1480.532
## m.sat.h 2 5 2956.329 2978.004 -1473.165 1 vs 2 14.7354 1e-04
summary(m.sat.h)
## Linear mixed-effects model fit by maximum likelihood
## Data: NULL
## AIC BIC logLik
## 2956.329 2978.004 -1473.164
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 1.592001 2.778264
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | iv
## Parameter estimates:
## NonHierarchical Hierarchical
## 1.000000 1.153484
## Fixed effects: dv ~ iv
## Value Std.Error DF t-value p-value
## (Intercept) 15.252474 0.2618853 339 58.24106 0e+00
## ivNonHierarchical 1.317279 0.3380101 223 3.89716 1e-04
## Correlation:
## (Intr)
## ivNonHierarchical -0.775
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.72377209 -0.62639807 0.09877309 0.69605258 1.86197608
##
## Number of Observations: 564
## Number of Groups: 225
pvals <- c(pvals, a$`p-value`)
Adding hierarchy type as a fixed effect improved the variance model (no fixed effects) for satisfaction.
Step 2: control variables
Checking which potential control variables are related to satisfaction:
for(x in cvs){
data <- na.omit(study1.long[c(vars,x)])
print(x)
dv <- data$satisfaction
iv <- data$hierarchy
id <- data$Respondent.ID
cv <- data[[x]]
m1 <- lme(dv ~ 1, random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~ 1 | iv),
method="ML")
mcv <- lme(dv ~ cv, random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~ 1 | iv),
method="ML")
a <- anova(m1,mcv)
print(a)
pvals <- c(pvals, a$`p-value`)
}
## [1] "rel.length.c"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 4 2941.305 2958.602 -1466.652
## mcv 2 5 2930.017 2951.638 -1460.008 1 vs 2 13.28801 3e-04
## [1] "age.c"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 4 2969.064 2986.405 -1480.532
## mcv 2 5 2970.799 2992.474 -1480.399 1 vs 2 0.2653941 0.6064
## [1] "incomeperPPL.c"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 4 2306.772 2323.083 -1149.386
## mcv 2 5 2308.769 2329.158 -1149.385 1 vs 2 0.002634695 0.9591
## [1] "N.Partners.3levels"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 4 2969.064 2986.405 -1480.532
## mcv 2 6 2965.087 2991.097 -1476.543 1 vs 2 7.977476 0.0185
## [1] "gender"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 4 2461.340 2477.917 -1226.670
## mcv 2 6 2460.173 2485.039 -1224.087 1 vs 2 5.166593 0.0755
## [1] "orientation.binary"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 4 2430.531 2447.065 -1211.266
## mcv 2 5 2430.610 2451.276 -1210.305 1 vs 2 1.921985 0.1656
## [1] "race.binary"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 4 2448.238 2464.797 -1220.119
## mcv 2 5 2450.140 2470.840 -1220.070 1 vs 2 0.09763496 0.7547
## [1] "children.binary"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 4 2461.340 2477.917 -1226.670
## mcv 2 5 2460.307 2481.028 -1225.153 1 vs 2 3.033148 0.0816
## [1] "education.3levels"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 4 2419.185 2435.701 -1205.592
## mcv 2 6 2421.665 2446.439 -1204.833 1 vs 2 1.519473 0.4678
## [1] "marital.status"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 4 2969.064 2986.405 -1480.532
## mcv 2 5 2960.798 2982.473 -1475.399 1 vs 2 10.26638 0.0014
## [1] "cohab"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 4 2969.064 2986.405 -1480.532
## mcv 2 6 2956.537 2982.548 -1472.269 1 vs 2 16.52699 3e-04
## [1] "coparent"
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 4 2969.064 2986.405 -1480.532
## mcv 2 6 2962.505 2988.515 -1475.252 1 vs 2 10.55977 0.0051
The following variables were related to satisfaction: relationship length, marital status, cohabitation, co-parenting.
Model selection strategy: start with a model that includes hierarchy (iv of interest), relationship length, marital status, cohabitation, co-parenting, all as fixed effects variables. Then perform model selection through a backwards elimination process.
cvs.sat.fix <- c("rel.length.c", "marital.status", "cohab", "coparent")
data <- na.omit(study1.long[c(vars,cvs.sat.fix)])
dv <- data$satisfaction
iv <- data$hierarchy
id <- data$Respondent.ID
cv1 <- data$rel.length.c
cv10 <- data$marital.status
cv11 <- data$cohab
cv12 <- data$coparent
aic0 <- numeric()
m0 <- lme(dv ~ iv + cv1 + cv10 + cv11 + cv12,
random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~ 1 | iv),
method="ML")
print(cvs.sat.fix[1])
## [1] "rel.length.c"
m <- update(m0, .~. -cv1)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[1] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 10 2918.595 2961.839 -1449.298
## m0 2 11 2918.783 2966.351 -1448.391 1 vs 2 1.8124 0.1782
pvals <- c(pvals, a$`p-value`)
print(cvs.sat.fix[2])
## [1] "marital.status"
m <- update(m0, .~. -cv10)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[2] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 10 2916.783 2960.027 -1448.392
## m0 2 11 2918.783 2966.351 -1448.391 1 vs 2 0.0003488189 0.9851
pvals <- c(pvals, a$`p-value`)
print(cvs.sat.fix[3])
## [1] "cohab"
m <- update(m0, .~. -cv11)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[3] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 9 2918.967 2957.886 -1450.483
## m0 2 11 2918.783 2966.351 -1448.391 1 vs 2 4.184037 0.1234
pvals <- c(pvals, a$`p-value`)
print(cvs.sat.fix[4])
## [1] "coparent"
m <- update(m0, .~. -cv12)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[4] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 9 2915.791 2954.711 -1448.896
## m0 2 11 2918.783 2966.351 -1448.391 1 vs 2 1.008725 0.6039
pvals <- c(pvals, a$`p-value`)
print(aic0)
## [1] 0.1876002 1.9996512 -0.1840373 2.9912755
Remove coparenting - it does not improve model with 4 control variables and increases AIC the most.
cvs.sat.fix <- c("rel.length.c", "marital.status", "cohab")
data <- na.omit(study1.long[c(vars,cvs.sat.fix)])
dv <- data$satisfaction
iv <- data$hierarchy
id <- data$Respondent.ID
cv1 <- data$rel.length.c
cv10 <- data$marital.status
cv11 <- data$cohab
aic0 <- numeric()
m0 <- lme(dv ~ iv + cv1 + cv10 + cv11,
random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~ 1 | iv),
method="ML")
print(cvs.sat.fix[1])
## [1] "rel.length.c"
m <- update(m0, .~. -cv1)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[1] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 8 2916.591 2951.186 -1450.295
## m0 2 9 2915.791 2954.711 -1448.896 1 vs 2 2.799484 0.0943
pvals <- c(pvals, a$`p-value`)
print(cvs.sat.fix[2])
## [1] "marital.status"
m <- update(m0, .~. -cv10)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[2] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 8 2913.801 2948.396 -1448.900
## m0 2 9 2915.791 2954.711 -1448.896 1 vs 2 0.009129234 0.9239
pvals <- c(pvals, a$`p-value`)
print(cvs.sat.fix[3])
## [1] "cohab"
m <- update(m0, .~. -cv11)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[3] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 7 2916.701 2946.972 -1451.351
## m0 2 9 2915.791 2954.711 -1448.896 1 vs 2 4.909633 0.0859
pvals <- c(pvals, a$`p-value`)
print(aic0)
## [1] -0.7994839 1.9908708 -0.9096331
Remove marital status.
cvs.sat.fix <- c("rel.length.c", "cohab")
data <- na.omit(study1.long[c(vars,cvs.sat.fix)])
dv <- data$satisfaction
iv <- data$hierarchy
id <- data$Respondent.ID
cv1 <- data$rel.length.c
cv11 <- data$cohab
aic0 <- numeric()
m0 <- lme(dv ~ iv + cv1 + cv11,
random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~ 1 | iv),
method="ML")
print(cvs.sat.fix[1])
## [1] "rel.length.c"
m <- update(m0, .~. -cv1)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[1] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 7 2916.617 2946.888 -1451.309
## m0 2 8 2913.801 2948.396 -1448.900 1 vs 2 4.816656 0.0282
pvals <- c(pvals, a$`p-value`)
print(cvs.sat.fix[2])
## [1] "cohab"
m <- update(m0, .~. -cv11)
a <- anova(m,m0)
diff <- a$AIC[2] - a$AIC[1]
aic0[2] <- diff
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 6 2915.594 2941.541 -1451.797
## m0 2 8 2913.801 2948.396 -1448.900 1 vs 2 5.793781 0.0552
pvals <- c(pvals, a$`p-value`)
print(aic0)
## [1] -2.816656 -1.793781
Revome cohabitation. It improves the model, but not significantly.
cvs.sat.fix <- c("rel.length.c")
data <- na.omit(study1.long[c(vars,cvs.sat.fix)])
dv <- data$satisfaction
iv <- data$hierarchy
id <- data$Respondent.ID
cv1 <- data$rel.length.c
m0 <- lme(dv ~ iv + cv1,
random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~ 1 | iv),
method="ML")
print(cvs.sat.fix[1])
## [1] "rel.length.c"
m <- update(m0, .~. -cv1)
a <- anova(m,m0)
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 5 2929.017 2950.639 -1459.508
## m0 2 6 2915.594 2941.541 -1451.797 1 vs 2 15.42234 1e-04
pvals <- c(pvals, a$`p-value`)
m <- update(m0, .~. -iv)
a <- anova(m,m0)
print(a)
## Model df AIC BIC logLik Test L.Ratio p-value
## m 1 5 2930.017 2951.638 -1460.008
## m0 2 6 2915.594 2941.541 -1451.797 1 vs 2 16.42221 1e-04
pvals <- c(pvals, a$`p-value`)
summary(m0)
## Linear mixed-effects model fit by maximum likelihood
## Data: NULL
## AIC BIC logLik
## 2915.594 2941.541 -1451.797
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 1.596808 2.751665
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | iv
## Parameter estimates:
## NonHierarchical Hierarchical
## 1.000000 1.145703
## Fixed effects: dv ~ iv + cv1
## Value Std.Error DF t-value p-value
## (Intercept) 15.210669 0.2616928 332 58.12414 0e+00
## ivNonHierarchical 1.392232 0.3382513 223 4.11597 1e-04
## cv1 0.539746 0.1366961 332 3.94851 1e-04
## Correlation:
## (Intr) ivNnHr
## ivNonHierarchical -0.775
## cv1 -0.054 0.066
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.6063417 -0.6053108 0.1151261 0.6845231 1.9090858
##
## Number of Observations: 558
## Number of Groups: 225
Residual plots:
res <- resid(m0, type = "pearson")
fit <- fitted(m0, level=0)
plot(ranef(m0, level = 1))
plot(res ~ iv)
plot(res ~ cv1) # rel. length
plot(res ~ cv11) # cohab
qqnorm(res)
qqline(res)
Differences in satisfaction between non-hierarchical relationships and 3 levels of hierarchy (primary, secondary, and tertiary partners)
data <- na.omit(study1.long[c(vars,"hierarchy2","rel.length.c","cohab")])
dv <- data$satisfaction
iv <- data$hierarchy
iv2 <- data$hierarchy2
id <- data$Respondent.ID
cv1 <- data$rel.length.c
m.sat <- lme(dv ~ 1, random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~ 1 | iv),
method="ML")
m.sat.h <- lme(dv ~ iv2, random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~ 1 | iv),
method="ML")
m.sat.hcv <- lme(dv ~ iv2 + cv1, random = list(id = pdSymm(form = ~ 1)),
weights = varIdent(form = ~ 1 | iv),
method="ML")
print(anova(m.sat, m.sat.h, m.sat.hcv))
## Model df AIC BIC logLik Test L.Ratio p-value
## m.sat 1 4 2941.305 2958.602 -1466.652
## m.sat.h 2 8 2878.146 2912.741 -1431.073 1 vs 2 71.15881 <.0001
## m.sat.hcv 3 9 2878.906 2917.825 -1430.453 2 vs 3 1.24020 0.2654
summary(m.sat.h)
## Linear mixed-effects model fit by maximum likelihood
## Data: NULL
## AIC BIC logLik
## 2878.146 2912.741 -1431.073
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 1.649997 2.776589
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | iv
## Parameter estimates:
## NonHierarchical Hierarchical
## 1.0000000 0.9995233
## Fixed effects: dv ~ iv2
## Value Std.Error DF t-value p-value
## (Intercept) 16.571627 0.2182967 329 75.91332 0.0000
## iv2Other -2.287679 0.7784189 329 -2.93888 0.0035
## iv2Primary 0.260075 0.3927875 329 0.66213 0.5084
## iv2Secondary -1.982843 0.3984492 329 -4.97640 0.0000
## iv2Tertiary -4.446815 0.6813313 329 -6.52666 0.0000
## Correlation:
## (Intr) iv2Oth iv2Prm iv2Scn
## iv2Other -0.280
## iv2Primary -0.556 0.251
## iv2Secondary -0.548 0.219 0.483
## iv2Tertiary -0.320 0.116 0.276 0.256
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.71423031 -0.57983191 0.09572108 0.66164997 1.90656048
##
## Number of Observations: 558
## Number of Groups: 225
summary(m.sat.hcv)
## Linear mixed-effects model fit by maximum likelihood
## Data: NULL
## AIC BIC logLik
## 2878.906 2917.825 -1430.453
##
## Random effects:
## Formula: ~1 | id
## (Intercept) Residual
## StdDev: 1.638991 2.759479
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | iv
## Parameter estimates:
## NonHierarchical Hierarchical
## 1.000000 1.013855
## Fixed effects: dv ~ iv2 + cv1
## Value Std.Error DF t-value p-value
## (Intercept) 16.581345 0.2172686 328 76.31725 0.0000
## iv2Other -2.249935 0.7833320 328 -2.87226 0.0043
## iv2Primary 0.110002 0.4151829 328 0.26495 0.7912
## iv2Secondary -1.913839 0.4035860 328 -4.74208 0.0000
## iv2Tertiary -4.392910 0.6865104 328 -6.39890 0.0000
## cv1 0.163908 0.1454126 328 1.12719 0.2605
## Correlation:
## (Intr) iv2Oth iv2Prm iv2Scn iv2Trt
## iv2Other -0.275
## iv2Primary -0.535 0.217
## iv2Secondary -0.532 0.219 0.399
## iv2Tertiary -0.313 0.116 0.232 0.259
## cv1 0.040 0.049 -0.319 0.148 0.073
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.6920637 -0.5945191 0.1124667 0.6491956 1.9112982
##
## Number of Observations: 558
## Number of Groups: 225
After controling for relationship length, satisfaction in non-hierarchical relationships is higher, on average, than satisfaction in secondary and tertiary relationships, but is at a similar level as primary ones.