Read in the individual data (or a pairwise dataset)
library(tidyr)
library(dplyr)
#install.packages("lme4")
library(lme4)
acitelli_ind <- read.csv(file.choose(), header=TRUE)Convert individual data to pairwise. I also create a simhobs variable that will be our binary response and two dummary variables that will be useful for estimating separate random intercepts for men and women.
tempA <- acitelli_ind %>%
mutate(genderE = gender, partnum = 1) %>%
mutate(gender = ifelse(gender == 1, "A", "P")) %>%
gather(variable, value, self_pos:genderE) %>%
unite(var_gender, variable, gender) %>%
spread(var_gender, value)
tempB <- acitelli_ind %>%
mutate(genderE = gender, partnum = 2) %>%
mutate(gender = ifelse(gender == 1, "P", "A")) %>%
gather(variable, value, self_pos:genderE)%>%
unite(var_gender, variable, gender) %>%
spread(var_gender, value)
acitelli_pair <- bind_rows(tempA, tempB) %>%
arrange(cuplid) %>%
mutate(gender_A = ifelse(genderE_A == 1, "hus", "wife"),
gender_A = as.factor(gender_A),
simhob_bin_A = ifelse(simhob_A == 1, 1, 0),
man = ifelse(genderE_A == 1, 1, 0),
woman = ifelse(genderE_A == 1, 0, 1))
rm(tempA, tempB)apim_bin <- glmer(simhob_bin_A ~ other_pos_A + other_pos_P
+ (1|cuplid),
data = acitelli_pair,
family = binomial,
na.action = na.omit)
summary(apim_bin)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: simhob_bin_A ~ other_pos_A + other_pos_P + (1 | cuplid)
## Data: acitelli_pair
##
## AIC BIC logLik deviance df.resid
## 322.5 337.3 -157.3 314.5 292
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.90316 -0.44848 -0.33245 0.02569 2.03565
##
## Random effects:
## Groups Name Variance Std.Dev.
## cuplid (Intercept) 1.659 1.288
## Number of obs: 296, groups: cuplid, 148
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.3310 2.4903 -3.345 0.000821 ***
## other_pos_A 0.9195 0.3894 2.361 0.018217 *
## other_pos_P 0.6723 0.3756 1.790 0.073481 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) oth__A
## other_pos_A -0.738
## other_pos_P -0.708 0.058Interaction approach.
apim_bin_di <- glmer(simhob_bin_A ~ other_pos_A + other_pos_P + genderE_A
+ other_pos_A*genderE_A + other_pos_P*genderE_A
+ (man + woman - 1|cuplid),
data = acitelli_pair,
family = binomial,
na.action = na.omit)
summary(apim_bin_di)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## simhob_bin_A ~ other_pos_A + other_pos_P + genderE_A + other_pos_A *
## genderE_A + other_pos_P * genderE_A + (man + woman - 1 | cuplid)
## Data: acitelli_pair
##
## AIC BIC logLik deviance df.resid
## 314.7 347.9 -148.4 296.7 287
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.09002 -0.48703 -0.00949 0.03292 2.52080
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## cuplid man 183.878 13.560
## woman 0.191 0.437 1.00
## Number of obs: 296, groups: cuplid, 148
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.842e+01 6.918e-04 -55536 <2e-16 ***
## other_pos_A 6.104e+00 6.917e-04 8824 <2e-16 ***
## other_pos_P 1.435e+00 6.917e-04 2074 <2e-16 ***
## genderE_A -3.062e+01 6.918e-04 -44265 <2e-16 ***
## other_pos_A:genderE_A 5.015e+00 6.917e-04 7250 <2e-16 ***
## other_pos_P:genderE_A 9.391e-01 6.917e-04 1358 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) oth__A oth__P gndE_A o__A:E
## other_pos_A 0.000
## other_pos_P 0.000 0.000
## genderE_A 0.000 0.000 0.000
## othr__A:E_A 0.000 0.000 0.000 0.000
## othr__P:E_A 0.000 0.000 0.000 0.000 0.000
## convergence code: 0
## Model failed to converge with max|grad| = 0.103721 (tol = 0.001, component 1)
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?Two-intercept model.
apim_bin_di_two <- glmer(simhob_bin_A ~ gender_A + other_pos_A:gender_A + other_pos_P:gender_A - 1
+ (man + woman - 1|cuplid),
data = acitelli_pair,
family = binomial,
na.action = na.omit)
summary(apim_bin_di_two)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## simhob_bin_A ~ gender_A + other_pos_A:gender_A + other_pos_P:gender_A -
## 1 + (man + woman - 1 | cuplid)
## Data: acitelli_pair
##
## AIC BIC logLik deviance df.resid
## 252.6 285.8 -117.3 234.6 287
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.013522 -0.004989 -0.000136 0.002841 0.134981
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## cuplid man 172856.6 415.76
## woman 860.6 29.34 0.02
## Number of obs: 296, groups: cuplid, 148
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## gender_Ahus -27.8603 59.6381 -0.467 0.640
## gender_Awife -19.3113 18.2034 -1.061 0.289
## gender_Ahus:other_pos_A 1.8910 7.4050 0.255 0.798
## gender_Awife:other_pos_A 1.5086 1.8566 0.813 0.416
## gender_Ahus:other_pos_P 0.2470 2.8218 0.088 0.930
## gender_Awife:other_pos_P 0.6637 2.7059 0.245 0.806
##
## Correlation of Fixed Effects:
## gndr_Ah gndr_Aw gndr_Ah:__A gndr_Aw:__A gndr_Ah:__P
## gender_Awif 0.810
## gndr_Ah:__A -0.932 -0.726
## gndr_Aw:__A -0.028 -0.375 -0.003
## gndr_Ah:__P -0.332 -0.268 0.029 0.150
## gndr_Aw:__P -0.737 -0.847 0.681 -0.122 0.179
## convergence code: 0
## unable to evaluate scaled gradient
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues
## failure to converge in 10000 evaluations