Install packages
install.packages("tidyverse", repos = "https://cloud.r-project.org")
##
## The downloaded binary packages are in
## /var/folders/5_/389qrkvs1sd7nkp792bslx5r0000gn/T//RtmpxHQjjJ/downloaded_packages
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.1 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
install.packages("psych", repos = "https://cloud.r-project.org")
##
## The downloaded binary packages are in
## /var/folders/5_/389qrkvs1sd7nkp792bslx5r0000gn/T//RtmpxHQjjJ/downloaded_packages
library(psych)
##
## Attaching package: 'psych'
##
## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha
install.packages("knitr",repos = "https://cloud.r-project.org")
##
## The downloaded binary packages are in
## /var/folders/5_/389qrkvs1sd7nkp792bslx5r0000gn/T//RtmpxHQjjJ/downloaded_packages
library("knitr")
install.packages("kableExtra",repos = "https://cloud.r-project.org")
##
## The downloaded binary packages are in
## /var/folders/5_/389qrkvs1sd7nkp792bslx5r0000gn/T//RtmpxHQjjJ/downloaded_packages
library("kableExtra")
##
## Attaching package: 'kableExtra'
##
## The following object is masked from 'package:dplyr':
##
## group_rows
install.packages("apaTables", repos = "https://cloud.r-project.org")
##
## The downloaded binary packages are in
## /var/folders/5_/389qrkvs1sd7nkp792bslx5r0000gn/T//RtmpxHQjjJ/downloaded_packages
library(apaTables)
install.packages("modelsummary", repos = "https://cloud.r-project.org")
##
## The downloaded binary packages are in
## /var/folders/5_/389qrkvs1sd7nkp792bslx5r0000gn/T//RtmpxHQjjJ/downloaded_packages
library("modelsummary")
##
## Attaching package: 'modelsummary'
##
## The following object is masked from 'package:psych':
##
## SD
install.packages("margins", repos = "https://cloud.r-project.org")
##
## The downloaded binary packages are in
## /var/folders/5_/389qrkvs1sd7nkp792bslx5r0000gn/T//RtmpxHQjjJ/downloaded_packages
library(margins)
Download data
dt <- read_csv("affair.csv")
## Rows: 6366 Columns: 9
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (9): affair, rate_marriage, age, yrs_married, children, religious, educ,...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Creating summary of statistics
summary_table1 <- describe(dt)
table1 <- knitr::kable(summary_table1,
"html",
caption = "Summary Statistics for Extramarital Affair Data",
digits = 2) %>%
footnote(general = "N = 6366") %>%
kable_styling(font_size = 12)
table1
Summary Statistics for Extramarital Affair Data
|
|
vars
|
n
|
mean
|
sd
|
median
|
trimmed
|
mad
|
min
|
max
|
range
|
skew
|
kurtosis
|
se
|
|
affair
|
1
|
6366
|
0.32
|
0.47
|
0
|
0.28
|
0.00
|
0.0
|
1.0
|
1.0
|
0.76
|
-1.42
|
0.01
|
|
rate_marriage
|
2
|
6366
|
4.11
|
0.96
|
4
|
4.24
|
1.48
|
1.0
|
5.0
|
4.0
|
-1.01
|
0.56
|
0.01
|
|
age
|
3
|
6366
|
29.08
|
6.85
|
27
|
28.48
|
7.41
|
17.5
|
42.0
|
24.5
|
0.58
|
-0.73
|
0.09
|
|
yrs_married
|
4
|
6366
|
9.01
|
7.28
|
6
|
8.22
|
5.19
|
0.5
|
23.0
|
22.5
|
0.73
|
-0.79
|
0.09
|
|
children
|
5
|
6366
|
1.40
|
1.43
|
1
|
1.21
|
1.48
|
0.0
|
5.5
|
5.5
|
0.89
|
0.24
|
0.02
|
|
religious
|
6
|
6366
|
2.43
|
0.88
|
2
|
2.41
|
1.48
|
1.0
|
4.0
|
3.0
|
-0.03
|
-0.73
|
0.01
|
|
educ
|
7
|
6366
|
14.21
|
2.18
|
14
|
13.97
|
2.97
|
9.0
|
20.0
|
11.0
|
0.79
|
0.48
|
0.03
|
|
occupation
|
8
|
6366
|
3.42
|
0.94
|
3
|
3.39
|
1.48
|
1.0
|
6.0
|
5.0
|
0.34
|
-0.09
|
0.01
|
|
occupation_husb
|
9
|
6366
|
3.85
|
1.35
|
4
|
3.88
|
1.48
|
1.0
|
6.0
|
5.0
|
-0.38
|
-0.84
|
0.02
|
|
Note:
|
|
N = 6366
|
LOGIT MODEL switching variables to factor
dt$affair <- factor(dt$affair)
dt$rate_marriage <- factor(dt$rate_marriage)
dt$religious <- factor(dt$religious)
BASE MODEl
Logit_BASE <- glm(dt$affair ~ dt$rate_marriage, data = dt, family = "binomial")
Creating Model 1 with all variables
Logit_model1 <- glm(affair ~ rate_marriage + yrs_married + age + religious + occupation + occupation_husb + children + educ, data = dt, family = "binomial")
summary(Logit_model1)
##
## Call:
## glm(formula = affair ~ rate_marriage + yrs_married + age + religious +
## occupation + occupation_husb + children + educ, family = "binomial",
## data = dt)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0823 -0.8136 -0.5751 1.0030 2.4359
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.051001 0.367275 5.584 2.35e-08 ***
## rate_marriage2 -0.419786 0.264683 -1.586 0.112741
## rate_marriage3 -0.790169 0.247205 -3.196 0.001391 **
## rate_marriage4 -1.677971 0.242983 -6.906 4.99e-12 ***
## rate_marriage5 -2.410610 0.244000 -9.880 < 2e-16 ***
## yrs_married 0.110158 0.010950 10.060 < 2e-16 ***
## age -0.060078 0.010284 -5.842 5.16e-09 ***
## religious2 -0.327460 0.084888 -3.858 0.000115 ***
## religious3 -0.632303 0.085767 -7.372 1.68e-13 ***
## religious4 -1.289919 0.131682 -9.796 < 2e-16 ***
## occupation 0.157295 0.034054 4.619 3.86e-06 ***
## occupation_husb 0.013778 0.022972 0.600 0.548658
## children -0.002952 0.031637 -0.093 0.925655
## educ -0.036186 0.015563 -2.325 0.020068 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 8005.1 on 6365 degrees of freedom
## Residual deviance: 6923.8 on 6352 degrees of freedom
## AIC: 6951.8
##
## Number of Fisher Scoring iterations: 4
Model 2 Only for significant control variables
Logit_model2 <- glm(affair ~ rate_marriage + yrs_married + age + religious + occupation, data = dt, family = "binomial")
summary(Logit_model2)
##
## Call:
## glm(formula = affair ~ rate_marriage + yrs_married + age + religious +
## occupation, family = "binomial", data = dt)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0344 -0.8071 -0.5735 0.9980 2.4380
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.766379 0.339161 5.208 1.91e-07 ***
## rate_marriage2 -0.415685 0.264345 -1.573 0.115832
## rate_marriage3 -0.787952 0.246829 -3.192 0.001411 **
## rate_marriage4 -1.676600 0.242556 -6.912 4.77e-12 ***
## rate_marriage5 -2.412604 0.243529 -9.907 < 2e-16 ***
## yrs_married 0.115758 0.009420 12.288 < 2e-16 ***
## age -0.065333 0.009996 -6.536 6.31e-11 ***
## religious2 -0.316010 0.084589 -3.736 0.000187 ***
## religious3 -0.626016 0.085453 -7.326 2.37e-13 ***
## religious4 -1.299991 0.131218 -9.907 < 2e-16 ***
## occupation 0.133050 0.031539 4.219 2.46e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 8005.1 on 6365 degrees of freedom
## Residual deviance: 6929.4 on 6355 degrees of freedom
## AIC: 6951.4
##
## Number of Fisher Scoring iterations: 4
Testing to for the better model
anova(Logit_model2, Logit_model1, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: affair ~ rate_marriage + yrs_married + age + religious + occupation
## Model 2: affair ~ rate_marriage + yrs_married + age + religious + occupation +
## occupation_husb + children + educ
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 6355 6929.4
## 2 6352 6923.8 3 5.5314 0.1368
Creating Side-By-Side Table
SBST4 <- list("Logit Base Model" = Logit_BASE,
"Logit full Model" = Logit_model1,
"Logit Reduced Model" = Logit_model2)
modelsummary(SBST4,
fmt = 5,
stars = TRUE,
output = "table4.2.docx")
## Warning: Using `$` in model formulas can produce unexpected results. Specify your
## model using the `data` argument instead.
## Try: affair ~ rate_marriage,
## data = dt
AMEmodel2 <- summary(margins(Logit_model2))
AME_model1 <- summary(margins(Logit_model1))
SBST4.1 <- list("Logit full Model AME" = AME_model1,
"Logit Reduced Model AME" = AMEmodel2)
SBST4.1
## $`Logit full Model AME`
## factor AME SE z p lower upper
## age -0.0109 0.0019 -5.8933 0.0000 -0.0146 -0.0073
## children -0.0005 0.0058 -0.0933 0.9257 -0.0118 0.0108
## educ -0.0066 0.0028 -2.3283 0.0199 -0.0121 -0.0010
## occupation 0.0286 0.0062 4.6447 0.0000 0.0166 0.0407
## occupation_husb 0.0025 0.0042 0.5998 0.5486 -0.0057 0.0107
## rate_marriage2 -0.0889 0.0539 -1.6478 0.0994 -0.1946 0.0168
## rate_marriage3 -0.1733 0.0498 -3.4767 0.0005 -0.2710 -0.0756
## rate_marriage4 -0.3713 0.0485 -7.6546 0.0000 -0.4664 -0.2762
## rate_marriage5 -0.5008 0.0482 -10.3885 0.0000 -0.5953 -0.4063
## religious2 -0.0652 0.0170 -3.8303 0.0001 -0.0985 -0.0318
## religious3 -0.1214 0.0167 -7.2515 0.0000 -0.1542 -0.0886
## religious4 -0.2238 0.0207 -10.8062 0.0000 -0.2644 -0.1832
## yrs_married 0.0201 0.0019 10.3344 0.0000 0.0163 0.0239
##
## $`Logit Reduced Model AME`
## factor AME SE z p lower upper
## age -0.0119 0.0018 -6.6085 0.0000 -0.0154 -0.0084
## occupation 0.0243 0.0057 4.2381 0.0000 0.0130 0.0355
## rate_marriage2 -0.0880 0.0539 -1.6334 0.1024 -0.1936 0.0176
## rate_marriage3 -0.1729 0.0498 -3.4728 0.0005 -0.2705 -0.0753
## rate_marriage4 -0.3714 0.0484 -7.6668 0.0000 -0.4663 -0.2764
## rate_marriage5 -0.5016 0.0481 -10.4210 0.0000 -0.5959 -0.4072
## religious2 -0.0630 0.0170 -3.7095 0.0002 -0.0962 -0.0297
## religious3 -0.1203 0.0167 -7.2052 0.0000 -0.1530 -0.0876
## religious4 -0.2250 0.0206 -10.9423 0.0000 -0.2653 -0.1847
## yrs_married 0.0211 0.0016 12.7946 0.0000 0.0179 0.0243
table1.2 <- knitr::kable(AME_model1,
"html",
caption = "AME FULL model for Extramarital Affair Data",
digits = 5) %>%
footnote(general = "N = 6366") %>%
kable_styling(font_size = 12)
table1.2
AME FULL model for Extramarital Affair Data
|
factor
|
AME
|
SE
|
z
|
p
|
lower
|
upper
|
|
age
|
-0.01094
|
0.00186
|
-5.89325
|
0.00000
|
-0.01458
|
-0.00730
|
|
children
|
-0.00054
|
0.00576
|
-0.09331
|
0.92565
|
-0.01183
|
0.01075
|
|
educ
|
-0.00659
|
0.00283
|
-2.32834
|
0.01989
|
-0.01214
|
-0.00104
|
|
occupation
|
0.02864
|
0.00617
|
4.64468
|
0.00000
|
0.01656
|
0.04073
|
|
occupation_husb
|
0.00251
|
0.00418
|
0.59982
|
0.54862
|
-0.00569
|
0.01071
|
|
rate_marriage2
|
-0.08888
|
0.05394
|
-1.64784
|
0.09938
|
-0.19460
|
0.01684
|
|
rate_marriage3
|
-0.17330
|
0.04985
|
-3.47672
|
0.00051
|
-0.27100
|
-0.07561
|
|
rate_marriage4
|
-0.37130
|
0.04851
|
-7.65463
|
0.00000
|
-0.46637
|
-0.27623
|
|
rate_marriage5
|
-0.50082
|
0.04821
|
-10.38851
|
0.00000
|
-0.59530
|
-0.40633
|
|
religious2
|
-0.06517
|
0.01701
|
-3.83027
|
0.00013
|
-0.09852
|
-0.03182
|
|
religious3
|
-0.12140
|
0.01674
|
-7.25146
|
0.00000
|
-0.15422
|
-0.08859
|
|
religious4
|
-0.22382
|
0.02071
|
-10.80622
|
0.00000
|
-0.26441
|
-0.18322
|
|
yrs_married
|
0.02006
|
0.00194
|
10.33445
|
0.00000
|
0.01626
|
0.02386
|
|
Note:
|
|
N = 6366
|
table1.2 <- knitr::kable(AMEmodel2,
"html",
caption = "AME Reduced model for Extramarital Affair Data",
digits = 5) %>%
footnote(general = "N = 6366") %>%
kable_styling(font_size = 12)
table1.2
AME Reduced model for Extramarital Affair Data
|
factor
|
AME
|
SE
|
z
|
p
|
lower
|
upper
|
|
age
|
-0.01191
|
0.00180
|
-6.60853
|
0.00000
|
-0.01544
|
-0.00838
|
|
occupation
|
0.02425
|
0.00572
|
4.23806
|
0.00002
|
0.01304
|
0.03547
|
|
rate_marriage2
|
-0.08802
|
0.05389
|
-1.63344
|
0.10238
|
-0.19364
|
0.01760
|
|
rate_marriage3
|
-0.17293
|
0.04979
|
-3.47282
|
0.00052
|
-0.27052
|
-0.07533
|
|
rate_marriage4
|
-0.37138
|
0.04844
|
-7.66685
|
0.00000
|
-0.46632
|
-0.27644
|
|
rate_marriage5
|
-0.50157
|
0.04813
|
-10.42103
|
0.00000
|
-0.59590
|
-0.40723
|
|
religious2
|
-0.06297
|
0.01698
|
-3.70953
|
0.00021
|
-0.09625
|
-0.02970
|
|
religious3
|
-0.12027
|
0.01669
|
-7.20518
|
0.00000
|
-0.15299
|
-0.08755
|
|
religious4
|
-0.22499
|
0.02056
|
-10.94234
|
0.00000
|
-0.26528
|
-0.18469
|
|
yrs_married
|
0.02110
|
0.00165
|
12.79456
|
0.00000
|
0.01787
|
0.02433
|
|
Note:
|
|
N = 6366
|