** This is a course material for PSY6111-01-00. Any questions or comments regarding the material should be addressed to the course instructor, Jeong Eun Cheon (email: cheonje@yonsei.ac.kr) **
We will be utilizing data from Cheon et al. (2024). In our study, SEM was used to estimate the effects of both wives’ and husbands’ neuroticism on marital satisfaction. We aim to partially replicate these findings using MLM. Specifically, we will compare and contrast results using ordinary regression, SEM, and MLM. It is important to note that this exercise is not designed to perfectly replicate all findings from Cheon et al. (2024), but rather to serve as a teaching tool for exploring how dyadic data can be analyzed using MLM.
Reference:
Cheon, J. E., Yang, H., Kang, I., & Kim, Y. H. (2024). Totality model of neuroticism: The sum of spousal neuroticism and marital satisfaction. Personality and Individual Differences, 228, 112726.
library(foreign)
data_wide <-read.spss("Totality_Study1_wide.sav",use.value.label=FALSE, to.data.frame=TRUE)
names(data_wide)## [1] "DyadNo" "UID" "Wgender" "Wage" "Wmstate"
## [6] "Wchildnum" "Wn1" "Wn2_r" "Wn3" "Wn4_r"
## [11] "Wn5" "Wn6" "Wn7" "Wn8" "Wn9"
## [16] "Wn10" "WMSAT1" "WMSAT2" "WMSAT3" "WMSAT4"
## [21] "WMSAT5" "WMSAT6" "Hgender" "Hage" "Hchildnum"
## [26] "Hn1" "Hn2_r" "Hn3" "Hn4_r" "Hn5"
## [31] "Hn6" "Hn7" "Hn8" "Hn9" "Hn10"
## [36] "HMSAT1" "HMSAT2" "HMSAT3" "HMSAT4" "HMSAT5"
## [41] "HMSAT6" "Hn2r" "Hn4r" "Wn2r" "Wn4r"
## [46] "HN" "WN" "WMsat" "HMsat" "av_marlngth"names(data_long)## [1] "DyadNo" "UID" "av_marlngth" "role" "gender"
## [6] "age" "mstate" "childnum" "n1" "n2_r"
## [11] "n3" "n4_r" "n5" "n6" "n7"
## [16] "n8" "n9" "n10" "MSAT1" "MSAT2"
## [21] "MSAT3" "MSAT4" "MSAT5" "MSAT6" "n2r"
## [26] "n4r" "N" "Msat"summary(lm(Msat ~ N, data_long))##
## Call:
## lm(formula = Msat ~ N, data = data_long)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.7545 -0.4528 0.1595 0.6917 1.4714
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.65790 0.17042 39.067 < 0.0000000000000002 ***
## N -0.18822 0.03978 -4.731 0.00000309 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9948 on 406 degrees of freedom
## Multiple R-squared: 0.05225, Adjusted R-squared: 0.04991
## F-statistic: 22.38 on 1 and 406 DF, p-value: 0.000003089Research Question: Is there a gender difference in the strength of the association between neuroticism and marital satisfaction?
summary(lm(Msat ~ N*gender, data_long))##
## Call:
## lm(formula = Msat ~ N * gender, data = data_long)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.7680 -0.4535 0.1576 0.6875 1.4823
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.57653 0.55828 11.780 <0.0000000000000002 ***
## N -0.16487 0.14136 -1.166 0.244
## gender 0.05351 0.38227 0.140 0.889
## N:gender -0.01472 0.08929 -0.165 0.869
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9972 on 404 degrees of freedom
## Multiple R-squared: 0.05232, Adjusted R-squared: 0.04528
## F-statistic: 7.435 on 3 and 404 DF, p-value: 0.00007397No, the gender difference is not statistically significant.
First, let’s use Structural Equation Modeling (SEM) to estimate the effects of both actor and partner neuroticism on marital satisfaction through the Actor-Partner Interdependence Model (APIM).
library(lavaan)
APIM <- 'HMsat ~ ah*HN + ph*WN
WMsat ~ aw*WN + pw*HN
HMsat ~~ WMsat
HN ~~ WN'
APIMout <- sem(APIM, data = data_wide ,missing="ml")
summary(APIMout, fit.measures = TRUE, standardize=TRUE, rsquare=TRUE, ci = TRUE)## lavaan 0.6.17 ended normally after 28 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 14
##
## Number of observations 204
## Number of missing patterns 1
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
##
## Model Test Baseline Model:
##
## Test statistic 161.550
## Degrees of freedom 6
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 1.000
## Tucker-Lewis Index (TLI) 1.000
##
## Robust Comparative Fit Index (CFI) 1.000
## Robust Tucker-Lewis Index (TLI) 1.000
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1119.672
## Loglikelihood unrestricted model (H1) -1119.672
##
## Akaike (AIC) 2267.345
## Bayesian (BIC) 2313.798
## Sample-size adjusted Bayesian (SABIC) 2269.442
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.000
## P-value H_0: RMSEA <= 0.050 NA
## P-value H_0: RMSEA >= 0.080 NA
##
## Robust RMSEA 0.000
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.000
## P-value H_0: Robust RMSEA <= 0.050 NA
## P-value H_0: Robust RMSEA >= 0.080 NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.000
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Observed
## Observed information based on Hessian
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
## HMsat ~
## HN (ah) -0.181 0.056 -3.231 0.001 -0.290 -0.071
## WN (ph) -0.093 0.056 -1.671 0.095 -0.202 0.016
## WMsat ~
## WN (aw) -0.196 0.067 -2.912 0.004 -0.329 -0.064
## HN (pw) -0.190 0.068 -2.807 0.005 -0.323 -0.057
## Std.lv Std.all
##
## -0.181 -0.219
## -0.093 -0.113
##
## -0.196 -0.196
## -0.190 -0.189
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
## .HMsat ~~
## .WMsat 0.651 0.080 8.132 0.000 0.494 0.808
## HN ~~
## WN -0.013 0.086 -0.155 0.877 -0.181 0.155
## Std.lv Std.all
##
## 0.651 0.693
##
## -0.013 -0.011
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
## .HMsat 7.067 0.334 21.180 0.000 6.413 7.721
## .WMsat 7.367 0.404 18.222 0.000 6.575 8.159
## HN 3.545 0.077 45.881 0.000 3.393 3.696
## WN 4.657 0.078 60.018 0.000 4.505 4.809
## Std.lv Std.all
## 7.067 7.777
## 7.367 6.644
## 3.545 3.212
## 4.657 4.202
##
## Variances:
## Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
## .HMsat 0.776 0.077 10.100 0.000 0.625 0.926
## .WMsat 1.139 0.113 10.100 0.000 0.918 1.360
## HN 1.218 0.121 10.100 0.000 0.981 1.454
## WN 1.228 0.122 10.100 0.000 0.990 1.467
## Std.lv Std.all
## 0.776 0.940
## 1.139 0.926
## 1.218 1.000
## 1.228 1.000
##
## R-Square:
## Estimate
## HMsat 0.060
## WMsat 0.074#gender difference
APIM_G <- paste0(APIM, '
ah == aw
ph == pw')
#indistinguishable
APIM_IND <- paste0(APIM, '
ah == aw
ph == pw
#equal predictor variance
HN ~~ v2*HN
WN ~~ v2*WN
# equal outcome intercepts
HMsat ~ o1*1
WMsat ~ o1*1
# equal residual variances for outcome variables
HMsat ~~ v1*HMsat
WMsat ~~ v1*WMsat')
APIMout_G <- sem(APIM_G, data = data_wide ,missing="ml")
APIMout_IND <- sem(APIM_IND, data = data_wide ,missing="ml")
anova(APIMout, APIMout_G, APIMout_IND)##
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq)
## APIMout 0 2267.3 2313.8 0.0000
## APIMout_G 2 2266.2 2306.0 2.8073 2.8073 0.044483 2 0.2456972
## APIMout_IND 5 2280.0 2309.8 22.6412 19.8339 0.165850 3 0.0001837 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(APIMout_G, fit.measures = TRUE, standardize=TRUE, rsquare=TRUE, ci = TRUE)## lavaan 0.6.17 ended normally after 31 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 14
## Number of equality constraints 2
##
## Number of observations 204
## Number of missing patterns 1
##
## Model Test User Model:
##
## Test statistic 2.807
## Degrees of freedom 2
## P-value (Chi-square) 0.246
##
## Model Test Baseline Model:
##
## Test statistic 161.550
## Degrees of freedom 6
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.995
## Tucker-Lewis Index (TLI) 0.984
##
## Robust Comparative Fit Index (CFI) 0.995
## Robust Tucker-Lewis Index (TLI) 0.984
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1121.076
## Loglikelihood unrestricted model (H1) -1119.672
##
## Akaike (AIC) 2266.152
## Bayesian (BIC) 2305.969
## Sample-size adjusted Bayesian (SABIC) 2267.950
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.044
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.153
## P-value H_0: RMSEA <= 0.050 0.408
## P-value H_0: RMSEA >= 0.080 0.394
##
## Robust RMSEA 0.044
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.153
## P-value H_0: Robust RMSEA <= 0.050 0.408
## P-value H_0: Robust RMSEA >= 0.080 0.394
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.023
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Observed
## Observed information based on Hessian
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
## HMsat ~
## HN (ah) -0.174 0.043 -4.075 0.000 -0.258 -0.090
## WN (ph) -0.122 0.043 -2.857 0.004 -0.206 -0.038
## WMsat ~
## WN (aw) -0.174 0.043 -4.075 0.000 -0.258 -0.090
## HN (pw) -0.122 0.043 -2.857 0.004 -0.206 -0.038
## Std.lv Std.all
##
## -0.174 -0.211
## -0.122 -0.149
##
## -0.174 -0.176
## -0.122 -0.123
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
## .HMsat ~~
## .WMsat 0.651 0.080 8.104 0.000 0.493 0.808
## HN ~~
## WN -0.013 0.086 -0.155 0.877 -0.181 0.155
## Std.lv Std.all
##
## 0.651 0.690
##
## -0.013 -0.011
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
## .HMsat 7.180 0.328 21.860 0.000 6.536 7.824
## .WMsat 7.023 0.331 21.211 0.000 6.374 7.672
## HN 3.545 0.077 45.881 0.000 3.393 3.696
## WN 4.657 0.078 60.018 0.000 4.505 4.809
## Std.lv Std.all
## 7.180 7.873
## 7.023 6.411
## 3.545 3.212
## 4.657 4.202
##
## Variances:
## Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
## .HMsat 0.777 0.077 10.094 0.000 0.626 0.928
## .WMsat 1.145 0.114 10.070 0.000 0.922 1.368
## HN 1.218 0.121 10.100 0.000 0.981 1.454
## WN 1.228 0.122 10.100 0.000 0.990 1.467
## Std.lv Std.all
## 0.777 0.934
## 1.145 0.954
## 1.218 1.000
## 1.228 1.000
##
## R-Square:
## Estimate
## HMsat 0.066
## WMsat 0.046
##
## Constraints:
## |Slack|
## ah - (aw) 0.000
## ph - (pw) 0.000summary(APIMout_IND, fit.measures = TRUE, standardize=TRUE, rsquare=TRUE, ci = TRUE)## lavaan 0.6.17 ended normally after 22 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 14
## Number of equality constraints 5
##
## Number of observations 204
## Number of missing patterns 1
##
## Model Test User Model:
##
## Test statistic 22.641
## Degrees of freedom 5
## P-value (Chi-square) 0.000
##
## Model Test Baseline Model:
##
## Test statistic 161.550
## Degrees of freedom 6
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.887
## Tucker-Lewis Index (TLI) 0.864
##
## Robust Comparative Fit Index (CFI) 0.887
## Robust Tucker-Lewis Index (TLI) 0.864
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1130.993
## Loglikelihood unrestricted model (H1) -1119.672
##
## Akaike (AIC) 2279.986
## Bayesian (BIC) 2309.849
## Sample-size adjusted Bayesian (SABIC) 2281.334
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.132
## 90 Percent confidence interval - lower 0.080
## 90 Percent confidence interval - upper 0.189
## P-value H_0: RMSEA <= 0.050 0.007
## P-value H_0: RMSEA >= 0.080 0.949
##
## Robust RMSEA 0.132
## 90 Percent confidence interval - lower 0.080
## 90 Percent confidence interval - upper 0.189
## P-value H_0: Robust RMSEA <= 0.050 0.007
## P-value H_0: Robust RMSEA >= 0.080 0.949
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.086
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Observed
## Observed information based on Hessian
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
## HMsat ~
## HN (ah) -0.213 0.043 -4.959 0.000 -0.297 -0.129
## WN (ph) -0.117 0.043 -2.725 0.006 -0.201 -0.033
## WMsat ~
## WN (aw) -0.213 0.043 -4.959 0.000 -0.297 -0.129
## HN (pw) -0.117 0.043 -2.725 0.006 -0.201 -0.033
## Std.lv Std.all
##
## -0.213 -0.231
## -0.117 -0.127
##
## -0.213 -0.231
## -0.117 -0.127
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
## .HMsat ~~
## .WMsat 0.646 0.081 7.944 0.000 0.486 0.805
## HN ~~
## WN -0.013 0.086 -0.155 0.877 -0.181 0.155
## Std.lv Std.all
##
## 0.646 0.669
##
## -0.013 -0.011
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
## .HMsat (o1) 7.239 0.337 21.467 0.000 6.578 7.900
## .WMsat (o1) 7.239 0.337 21.467 0.000 6.578 7.900
## HN 3.545 0.077 45.779 0.000 3.393 3.696
## WN 4.657 0.077 60.151 0.000 4.506 4.809
## Std.lv Std.all
## 7.239 7.111
## 7.239 7.111
## 3.545 3.205
## 4.657 4.211
##
## Variances:
## Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
## HN (v2) 1.223 0.086 14.282 0.000 1.055 1.391
## WN (v2) 1.223 0.086 14.282 0.000 1.055 1.391
## .HMsat (v1) 0.965 0.081 11.870 0.000 0.805 1.124
## .WMsat (v1) 0.965 0.081 11.870 0.000 0.805 1.124
## Std.lv Std.all
## 1.223 1.000
## 1.223 1.000
## 0.965 0.931
## 0.965 0.931
##
## R-Square:
## Estimate
## HMsat 0.069
## WMsat 0.069
##
## Constraints:
## |Slack|
## ah - (aw) 0.000
## ph - (pw) 0.000Transition to MLM Analysis:
First, we need to adjust the data format for MLM analysis. Let’s convert it to a pairwise format.
# Create the female (Gender == 1) and male (Gender == 0) datasets
female_data <- data_long %>%
mutate(partnum = 1) %>%
gather(variable, value, -DyadNo, -gender) %>%
mutate(gender = ifelse(gender == 2, "A", "P")) %>%
unite(var_gender, variable, gender) %>%
spread(var_gender, value)
male_data <- data_long %>%
mutate(partnum = 2) %>%
gather(variable, value, -DyadNo, -gender) %>%
mutate(gender = ifelse(gender == 2, "P", "A")) %>%
unite(var_gender, variable, gender) %>%
spread(var_gender, value)
# Combine female and male data
dyad_pairwise <- bind_rows(female_data, male_data) %>%
arrange(DyadNo) %>%
dplyr::select(DyadNo, everything())
names(dyad_pairwise)## [1] "DyadNo" "age_A" "age_P" "av_marlngth_A"
## [5] "av_marlngth_P" "childnum_A" "childnum_P" "Msat_A"
## [9] "Msat_P" "MSAT1_A" "MSAT1_P" "MSAT2_A"
## [13] "MSAT2_P" "MSAT3_A" "MSAT3_P" "MSAT4_A"
## [17] "MSAT4_P" "MSAT5_A" "MSAT5_P" "MSAT6_A"
## [21] "MSAT6_P" "mstate_A" "mstate_P" "N_A"
## [25] "N_P" "n1_A" "n1_P" "n10_A"
## [29] "n10_P" "n2_r_A" "n2_r_P" "n2r_A"
## [33] "n2r_P" "n3_A" "n3_P" "n4_r_A"
## [37] "n4_r_P" "n4r_A" "n4r_P" "n5_A"
## [41] "n5_P" "n6_A" "n6_P" "n7_A"
## [45] "n7_P" "n8_A" "n8_P" "n9_A"
## [49] "n9_P" "partnum_A" "partnum_P" "role_A"
## [53] "role_P" "UID_A" "UID_P"# make gender variable
dyad_pairwise <- dyad_pairwise %>% mutate(gender_A = ifelse(partnum_A == 1, 1, 0))Before proceeding, let’s calculate the Intraclass Correlation Coefficient (ICC) to assess the importance of the dyadic nature of our data.
library(nlme)##
## Attaching package: 'nlme'## The following object is masked from 'package:usdm':
##
## Variogram## The following object is masked from 'package:dplyr':
##
## collapsedyad_pairwise <- dyad_pairwise %>%
mutate(across(Msat_A:N_P, as.numeric))
APIM_mlm <- lme(Msat_A ~ 1,
data = dyad_pairwise,
random = ~ 1 | DyadNo,
na.action = na.omit)
summary(APIM_mlm)## Linear mixed-effects model fit by REML
## Data: dyad_pairwise
## AIC BIC logLik
## 1058.288 1070.315 -526.1442
##
## Random effects:
## Formula: ~1 | DyadNo
## (Intercept) Residual
## StdDev: 0.8409876 0.5797626
##
## Fixed effects: Msat_A ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) 5.886029 0.06550417 204 89.85732 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.3165049 -0.3716188 0.0406183 0.3861620 1.9948664
##
## Number of Observations: 408
## Number of Groups: 2040.8409876/(0.8409876+0.5797626)## [1] 0.5919321#0.5919321Approximately 59% of the variance in marital satisfaction is attributable to differences between dyads.
Now, let’s proceed with the MLM analysis.
mlm <- gls(Msat_A ~ N_A + N_P,
correlation = corCompSymm(form = ~1|DyadNo),
data = dyad_pairwise)
summary(mlm)## Generalized least squares fit by REML
## Model: Msat_A ~ N_A + N_P
## Data: dyad_pairwise
## AIC BIC logLik
## 1045.529 1065.549 -517.7647
##
## Correlation Structure: Compound symmetry
## Formula: ~1 | DyadNo
## Parameter estimate(s):
## Rho
## 0.6706409
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 7.238668 0.3388598 21.361837 0.000
## N_A -0.212858 0.0431209 -4.936305 0.000
## N_P -0.116975 0.0431209 -2.712723 0.007
##
## Correlation:
## (Intr) N_A
## N_A -0.925
## N_P -0.925 0.772
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -4.5760872 -0.5124838 0.1753843 0.6926328 1.6064112
##
## Residual standard error: 0.9866725
## Degrees of freedom: 408 total; 405 residualAPIM_mlm <- lme(Msat_A ~ N_A + N_P,
data = dyad_pairwise,
random = ~ 1 | DyadNo,
na.action = na.omit)
summary(APIM_mlm)## Linear mixed-effects model fit by REML
## Data: dyad_pairwise
## AIC BIC logLik
## 1045.529 1065.549 -517.7647
##
## Random effects:
## Formula: ~1 | DyadNo
## (Intercept) Residual
## StdDev: 0.8080118 0.5662498
##
## Fixed effects: Msat_A ~ N_A + N_P
## Value Std.Error DF t-value p-value
## (Intercept) 7.238668 0.3388596 203 21.361849 0.0000
## N_A -0.212858 0.0431209 202 -4.936307 0.0000
## N_P -0.116975 0.0431209 202 -2.712724 0.0072
## Correlation:
## (Intr) N_A
## N_A -0.925
## N_P -0.925 0.772
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.35781824 -0.42960461 0.06726502 0.47476156 2.12284871
##
## Number of Observations: 408
## Number of Groups: 204#dyad_pairwise <- dyad_pairwise %>% rename(men = gender)
dyad_pairwise$women <- dyad_pairwise$gender_A #0 men, 1 women
dyad_pairwise <- dyad_pairwise %>% mutate(men = ifelse(women == 0, 1, 0)) #0 women, 1 men
APIM_mlm <- lme(Msat_A ~ men + women +
N_A:men + N_P:men +
N_A:women + N_P:women -1,
data = dyad_pairwise,
random = ~ -1 + men + women| DyadNo,
# correlation = corCompSymm(form=~1|DyadNo),# Further model within-dyad correlation
weights = varIdent(form=~1|women),
na.action = na.omit)
summary(APIM_mlm)## Linear mixed-effects model fit by REML
## Data: dyad_pairwise
## AIC BIC logLik
## 1045.359 1089.32 -511.6795
##
## Random effects:
## Formula: ~-1 + men + women | DyadNo
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## men 0.8319063 men
## women 1.0299387 0.771
## Residual 0.3086931
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | women
## Parameter estimates:
## 1 0
## 1 1
## Fixed effects: Msat_A ~ men + women + N_A:men + N_P:men + N_A:women + N_P:women - 1
## Value Std.Error DF t-value p-value
## men 7.066602 0.3361307 199 21.023377 0.0000
## women 7.366957 0.4072985 199 18.087366 0.0000
## men:N_A -0.180605 0.0563054 199 -3.207587 0.0016
## men:N_P -0.092964 0.0560561 199 -1.658409 0.0988
## women:N_A -0.196363 0.0679246 199 -2.890900 0.0043
## women:N_P -0.190100 0.0682268 199 -2.786291 0.0058
## Correlation:
## men women mn:N_A mn:N_P wm:N_A
## women 0.693
## men:N_A -0.602 -0.417
## men:N_P -0.783 -0.542 0.011
## women:N_A -0.542 -0.783 0.008 0.693
## women:N_P -0.417 -0.602 0.693 0.008 0.011
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.65584498 -0.25117227 0.05163055 0.24916563 1.30117986
##
## Number of Observations: 408
## Number of Groups: 204APIM_mlm <- lme(Msat_A ~ men + women +
N_A + N_P -1,
data = dyad_pairwise,
random = ~ -1 + men + women| DyadNo,
# correlation = corCompSymm(form=~1|DyadNo),
na.action = na.omit)
summary(APIM_mlm)## Linear mixed-effects model fit by REML
## Data: dyad_pairwise
## AIC BIC logLik
## 1035.48 1067.491 -509.7399
##
## Random effects:
## Formula: ~-1 + men + women | DyadNo
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## men 0.8304229 men
## women 1.0294364 0.769
## Residual 0.3091604
##
## Fixed effects: Msat_A ~ men + women + N_A + N_P - 1
## Value Std.Error DF t-value p-value
## men 7.180579 0.3280826 201 21.886496 0.0000
## women 7.023362 0.3308786 201 21.226403 0.0000
## N_A -0.174155 0.0427659 201 -4.072302 0.0001
## N_P -0.122345 0.0427299 201 -2.863206 0.0046
## Correlation:
## men women N_A
## women 0.979
## N_A -0.876 -0.915
## N_P -0.922 -0.869 0.683
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.58752287 -0.25037865 0.04519864 0.24634303 1.34869149
##
## Number of Observations: 408
## Number of Groups: 204