Read in the individual data (or a pairwise dataset)
library(tidyr)
library(dplyr)
#install.packages("lme4")
library(lme4)
acitelli_ind <- read.csv("/Users/randigarcia/Desktop/Data/acitelli.csv", header=TRUE)Convert individual data to pairwise. I also create a simhobs variable that will be our binary response, two dummy variables that will be useful for estimating separate random intercepts for men and women, and a count variable cigarettes.
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), #forced binary variable
man = ifelse(genderE_A == 1, 1, 0),
woman = ifelse(genderE_A == 1, 0, 1),
cigarettes_A = rpois(296, 0.7)) %>%
group_by(cuplid) %>%
mutate(cupcig = rpois(1, 0.7)) %>%
ungroup(cuplid) %>%
mutate(cigarettes_A = cigarettes_A + cupcig)
rm(tempA, tempB)To account for the nonindependence, we can make use of the glmer() function from the lme4 package. Note that we are asking for the variance of intercepts across dyads, that is the random intercept in traditional multilevel modeling. The gls() function in the nlme package does not have an option for specifying a link function (i.e., there is no family = option). The syntax of glmer() differs a bit from gls() in that the random effects are specified within the formula: + (1/cuplid).
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.09016 -0.48698 -0.00949 0.03293 2.52186
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## cuplid man 183.8755 13.5601
## woman 0.1908 0.4368 1.00
## Number of obs: 296, groups: cuplid, 148
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -38.413 15.095 -2.545 0.0109 *
## other_pos_A 6.104 2.453 2.489 0.0128 *
## other_pos_P 1.433 2.551 0.562 0.5744
## genderE_A -30.611 14.859 -2.060 0.0394 *
## other_pos_A:genderE_A 5.014 2.462 2.037 0.0417 *
## other_pos_P:genderE_A 0.937 2.559 0.366 0.7142
## ---
## 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.656
## other_pos_P -0.683 -0.098
## genderE_A 0.988 -0.647 -0.675
## othr__A:E_A -0.635 0.985 -0.113 -0.643
## othr__P:E_A -0.662 -0.114 0.986 -0.670 -0.134
## convergence code: 0
## unable to evaluate scaled gradient
## Model failed to converge: degenerate Hessian with 2 negative eigenvaluesTwo-intercept model.
#does not converge
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)
apim_bin_di_two <- glmer(simhob_bin_A ~ gender_A + other_pos_A:gender_A + other_pos_P:gender_A - 1
+ (1|cuplid),
data = acitelli_pair,
family = binomial,
na.action = na.omit)
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,
nAGQ = 0) #Adaptive Gauss-Hermite Quadrature
summary(apim_bin_di_two)## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 0) [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
## 330.4 363.6 -156.2 312.4 287
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.97710 -0.52321 -0.39856 0.00841 1.90671
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## cuplid man 1.257 1.121
## woman 1.038 1.019 1.00
## Number of obs: 296, groups: cuplid, 148
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## gender_Ahus -7.1637 2.7924 -2.566 0.0103 *
## gender_Awife -6.0917 2.4851 -2.451 0.0142 *
## gender_Ahus:other_pos_A 0.4259 0.5214 0.817 0.4141
## gender_Awife:other_pos_A 1.0173 0.4539 2.241 0.0250 *
## gender_Ahus:other_pos_P 0.9160 0.5003 1.831 0.0671 .
## gender_Awife:other_pos_P 0.1722 0.4655 0.370 0.7114
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## gndr_Ah gndr_Aw gndr_Ah:__A gndr_Aw:__A gndr_Ah:__P
## gender_Awif 0.153
## gndr_Ah:__A -0.641 -0.103
## gndr_Aw:__A -0.087 -0.611 -0.037
## gndr_Ah:__P -0.605 -0.090 -0.218 0.153
## gndr_Aw:__P -0.104 -0.629 0.164 -0.226 -0.036apim_poi <- glmer(cigarettes_A ~ other_pos_A + other_pos_P
+ (1|cuplid),
data = acitelli_pair,
family = poisson,
na.action = na.omit)
summary(apim_poi)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: cigarettes_A ~ other_pos_A + other_pos_P + (1 | cuplid)
## Data: acitelli_pair
##
## AIC BIC logLik deviance df.resid
## 859.4 874.2 -425.7 851.4 292
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.1350 -0.9103 -0.1137 0.5536 2.0678
##
## Random effects:
## Groups Name Variance Std.Dev.
## cuplid (Intercept) 0.1896 0.4354
## Number of obs: 296, groups: cuplid, 148
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.52660 0.68552 0.768 0.442
## other_pos_A -0.04827 0.11476 -0.421 0.674
## other_pos_P -0.03210 0.11493 -0.279 0.780
##
## Correlation of Fixed Effects:
## (Intr) oth__A
## other_pos_A -0.695
## other_pos_P -0.697 -0.022Interaction approach.
apim_poi_di <- glmer(cigarettes_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 = poisson,
na.action = na.omit)
summary(apim_poi_di)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## cigarettes_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
## 867.2 900.4 -424.6 849.2 287
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.1476 -0.8573 -0.0747 0.5109 2.0546
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## cuplid man 0.2391 0.4890
## woman 0.1516 0.3894 1.00
## Number of obs: 296, groups: cuplid, 148
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.55936 0.69255 0.808 0.419
## other_pos_A -0.04907 0.11664 -0.421 0.674
## other_pos_P -0.03971 0.11568 -0.343 0.731
## genderE_A -0.10278 0.55610 -0.185 0.853
## other_pos_A:genderE_A -0.05344 0.12225 -0.437 0.662
## other_pos_P:genderE_A 0.05875 0.12129 0.484 0.628
##
## Correlation of Fixed Effects:
## (Intr) oth__A oth__P gndE_A o__A:E
## other_pos_A -0.698
## other_pos_P -0.691 -0.025
## genderE_A 0.102 -0.134 -0.006
## othr__A:E_A -0.134 0.186 -0.004 -0.537
## othr__P:E_A 0.026 -0.044 0.013 -0.527 -0.428
## convergence code: 0
## Model failed to converge with max|grad| = 0.0114521 (tol = 0.001, component 1)Two-intercept model.
apim_poi_di_two <- glmer(cigarettes_A ~ gender_A + other_pos_A:gender_A + other_pos_P:gender_A - 1
+ (man + woman - 1|cuplid),
data = acitelli_pair,
family = poisson,
na.action = na.omit)
summary(apim_poi_di_two)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## cigarettes_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
## 867.2 900.4 -424.6 849.2 287
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.14761 -0.85705 -0.07483 0.51097 2.05495
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## cuplid man 0.2392 0.4891
## woman 0.1515 0.3893 1.00
## Number of obs: 296, groups: cuplid, 148
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## gender_Ahus 0.461750 0.931182 0.496 0.620
## gender_Awife 0.661191 0.842984 0.784 0.433
## gender_Ahus:other_pos_A -0.102970 0.183986 -0.560 0.576
## gender_Awife:other_pos_A 0.004639 0.152482 0.030 0.976
## gender_Ahus:other_pos_P 0.018262 0.168677 0.108 0.914
## gender_Awife:other_pos_P -0.098512 0.166535 -0.592 0.554
##
## Correlation of Fixed Effects:
## gndr_Ah gndr_Aw gndr_Ah:__A gndr_Aw:__A gndr_Ah:__P
## gender_Awif 0.217
## gndr_Ah:__A -0.659 -0.145
## gndr_Aw:__A -0.122 -0.567 -0.048
## gndr_Ah:__P -0.567 -0.121 -0.237 0.211
## gndr_Aw:__P -0.145 -0.660 0.215 -0.236 -0.047
## convergence code: 0
## Model failed to converge with max|grad| = 0.0121527 (tol = 0.001, component 1)#install.packages("gee")
library(gee)
apim_gee <- gee(simhob_bin_A ~ other_pos_A + other_pos_P,
id = cuplid,
data = acitelli_pair,
na.action = na.omit,
family = binomial,
corstr = "unstructured")## (Intercept) other_pos_A other_pos_P
## -6.2582893 0.6956734 0.5017679summary(apim_gee)##
## GEE: GENERALIZED LINEAR MODELS FOR DEPENDENT DATA
## gee S-function, version 4.13 modified 98/01/27 (1998)
##
## Model:
## Link: Logit
## Variance to Mean Relation: Binomial
## Correlation Structure: Unstructured
##
## Call:
## gee(formula = simhob_bin_A ~ other_pos_A + other_pos_P, id = cuplid,
## data = acitelli_pair, na.action = na.omit, family = binomial,
## corstr = "unstructured")
##
## Summary of Residuals:
## Min 1Q Median 3Q Max
## -0.43310561 -0.27624377 -0.19297150 0.07522507 0.86606888
##
##
## Coefficients:
## Estimate Naive S.E. Naive z Robust S.E. Robust z
## (Intercept) -6.2729871 1.8140043 -3.458088 1.7415077 -3.602044
## other_pos_A 0.6982148 0.2901570 2.406334 0.2851282 2.448775
## other_pos_P 0.5025443 0.2840654 1.769115 0.2605414 1.928846
##
## Estimated Scale Parameter: 0.9978716
## Number of Iterations: 2
##
## Working Correlation
## [,1] [,2]
## [1,] 1.0000000 0.2703261
## [2,] 0.2703261 1.0000000