library(haven)
GSS <- read_sav("/Users/loweticjustice/Desktop/New/GSS7216_R4.sav")
library(car)## Loading required package: carDatalibrary(stargazer)##
## Please cite as:## Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.2. https://CRAN.R-project.org/package=stargazerlibrary(survey)## Loading required package: grid## Loading required package: Matrix## Loading required package: survival##
## Attaching package: 'survey'## The following object is masked from 'package:graphics':
##
## dotchartlibrary(ggplot2)
library(pander)
library(knitr)
library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:car':
##
## recode## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(psych)##
## Attaching package: 'psych'## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha## The following object is masked from 'package:car':
##
## logitlibrary(lme4)
library(lmerTest) ##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## steplibrary(arm)## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select##
## arm (Version 1.10-1, built: 2018-4-12)## Working directory is /Users/loweticjustice/Library/Mobile Documents/com~apple~CloudDocs/Desktop##
## Attaching package: 'arm'## The following objects are masked from 'package:psych':
##
## logit, rescale, sim## The following object is masked from 'package:car':
##
## logitlibrary(questionr)##
## Attaching package: 'questionr'## The following object is masked from 'package:psych':
##
## describeGSS$age<-cut(GSS$AGE, breaks=c(0,24,39,59,79,99))
#Church Attendence and Religion
GSS$attend<-Recode(GSS$ATTEND,recodes="6:8='Very Frequent Church Attendee';0:2='Infequent Church Attendee'; 3:5='Frequent Church Attendee'; else=NA",as.factor=T)
GSS$relig<-Recode(GSS$CONCLERG,recodes="1='Strong'; 2='Moderate';3='Low'; else=NA", as.factor=T)
#Race
GSS$race<-Recode(GSS$RACE,recodes=" 2='Blacks'; 1='Whites'; else=NA", as.factor=T)
GSS$race<-relevel(GSS$race,ref="Blacks")
GSS$black<-Recode(GSS$RACE,recodes="2='1'; else=0", as.factor=T)
GSS$black<-Recode(GSS$RACE,recodes="2='black'; else=NA", as.factor=T)
#Sex
GSS$sex<-Recode(GSS$SEX,recodes="1='Male'; 2='Female'; else=NA", as.factor=T)
GSS$sex<-relevel(GSS$sex,ref="Male")
#Year in clusters
GSS$year<-Recode(GSS$YEAR, recodes="1972:1979='1972:1979'; 1980:1989='1980:1989'; 1990:1999='1990:1999';2000:2009='2000:2009';2010:2016='2010:2016';else=NA", as.factor=T)
GSS$year<-relevel(GSS$year,ref="1972:1979")
#Region
GSS$region1<-Recode(GSS$REG16,recodes="1='New England';2='Middle Atlantic';3:4='North Central';5:7='South';8='Mountain';9='Pacific';else=NA", as.factor=T)
GSS$region1<-relevel(GSS$region1,ref="New England")
GSS$region2<-Recode(GSS$REG16,recodes="5:7='South';0=NA;else='Other Regions'", as.factor=T)
GSS$region2<-relevel(GSS$region2,ref="Other Regions")
GSS$south<-Recode(GSS$REG16,recodes="5:7='South';1:4='North';else=NA", as.factor=T)
GSS$south<-relevel(GSS$south,ref="South")
#Resdential
GSS$resi<-Recode(GSS$RES16,recodes="1:3='Rural';4:6='Urban';else=NA", as.factor=T)
GSS$resi<-relevel(GSS$resi,ref="Urban")
#Education
GSS$educn<-Recode(GSS$DEGREE,recodes="0:1='HS or Less';2:4='Higher Education';else=NA", as.factor=T)
GSS$educn<-relevel(GSS$educn,ref="Higher Education")
#Education
GSS$educ<-Recode(GSS$DEGREE,recodes="0:1=0;2:4=1;else=NA", as.factor=T)
#Familiy
GSS$nchld<-Recode(GSS$CHILDS,recodes="0 = 'No child';2='Less than 2'; 3:8='3 or More'; else=NA", as.factor=T)
GSS$nchld<-relevel(GSS$nchld,ref="Less than 2")
GSS$marst<-Recode(GSS$MARITAL,recodes="1='Married'; 2='Widowed';3='Divorced';4='Seperated';5='Never Married'; else=NA", as.factor=T)
GSS$marst<-relevel(GSS$marst,ref="Married")
GSS$trust<-Recode(GSS$TRUST,recodes="1 = 'Trust';2='No Trust'; else=NA", as.factor=T)
GSS$trust<-relevel(GSS$trust,ref="Trust")
GSS$cohort<-Recode(GSS$COHORT, recodes="1883:1927='Greatest Gen'; 1928:1945='Silent'; 1946:1964='Boomers';1965:1980='Generation X';1981:1996='Millennials';else=NA", as.factor=T)
GSS$cohort<-relevel(GSS$cohort,ref="Millennials")GSS$cohorts<-Recode(GSS$COHORT, recodes="1883:1927=1;1928:1945=2;1946:1964=3;1965:1980=4;1981:1996=5;else=NA")
GSS$bgens<-ifelse((GSS$RACE==2)&(GSS$cohorts==1), "1",
ifelse((GSS$RACE==2)&(GSS$cohorts==2),"2",
ifelse((GSS$RACE==2)&(GSS$cohorts==3),"3", ifelse((GSS$RACE==2)&(GSS$cohorts==4),"4",
ifelse((GSS$RACE==2)&(GSS$cohorts==5),"5",NA)))))
GSS$bgens1<-Recode(GSS$bgens,recodes="'1'='Greatest Gen';'2'='Silent Gen';'3'='Boomers';'4'='Generation X';'5'='Millennials';else=NA",as.factor=T)
GSS$bgen<-relevel(GSS$bgens1,ref="Greatest Gen")
summary(GSS$bgen)## Greatest Gen Boomers Generation X Millennials Silent Gen
## 1304 3411 1702 619 1704
## NA's
## 53726table(GSS$bgen,GSS$year)##
## 1972:1979 1980:1989 1990:1999 2000:2009 2010:2016
## Greatest Gen 505 517 201 72 9
## Boomers 346 1001 759 811 494
## Generation X 0 109 410 760 423
## Millennials 0 0 0 197 422
## Silent Gen 385 557 356 280 126GSS$abtdefect<-Recode(GSS$ABDEFECT,recodes="1=1; 2=0; else=NA")
GSS$abtnomore<-Recode(GSS$ABNOMORE,recodes="1=1; 2=0; else=NA")
GSS$abtpoor<-Recode(GSS$ABPOOR,recodes="1=1; 2=0; else=NA")
GSS$abtrape<-Recode(GSS$ABRAPE,recodes="1=1; 2=0; else=NA")
GSS$abtsing<-Recode(GSS$ABSINGLE,recodes="1=1; 2=0; else=NA")
GSS$abthlth<-Recode(GSS$ABHLTH,recodes="1=1; 2=0; else=NA")
GSS$abscald <- scale(GSS$abtdefect)+scale(GSS$abtnomore)+scale(GSS$abtpoor)+scale(GSS$abtrape)+scale(GSS$abtsing)+scale(GSS$abthlth)
summary(GSS$abscald)## V1
## Min. :-9.739
## 1st Qu.:-1.879
## Median : 0.609
## Mean : 0.151
## 3rd Qu.: 4.627
## Max. : 4.627
## NA's :23793fit.an<-lm(abscald~as.factor(YEAR),GSS)
anova(fit.an)## Analysis of Variance Table
##
## Response: abscald
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(YEAR) 29 4722 162.825 7.7149 < 2.2e-16 ***
## Residuals 38643 815567 21.105
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1library(lme4)
library(lmerTest)
library(arm)
fit<-lmer(abscald~(1|YEAR), data=GSS)
arm::display(fit, detail=T)## lmer(formula = abscald ~ (1 | YEAR), data = GSS)
## coef.est coef.se t value
## (Intercept) 0.15 0.07 2.34
##
## Error terms:
## Groups Name Std.Dev.
## YEAR (Intercept) 0.33
## Residual 4.59
## ---
## number of obs: 38673, groups: YEAR, 30
## AIC = 227753, DIC = 227739.5
## deviance = 227743.1yearlymeans<-aggregate(cbind(abscald)~YEAR,GSS,mean)
head(yearlymeans,n=40)## YEAR V1
## 1 1972 -0.07576717
## 2 1973 0.66063573
## 3 1974 0.80742543
## 4 1975 0.57686615
## 5 1976 0.58085070
## 6 1977 0.61979604
## 7 1978 0.10691514
## 8 1980 0.55143642
## 9 1982 0.37560855
## 10 1983 -0.11345876
## 11 1984 0.08055414
## 12 1985 -0.12367390
## 13 1987 -0.20262923
## 14 1988 -0.21745744
## 15 1989 0.29832077
## 16 1990 0.47716396
## 17 1991 0.49642675
## 18 1993 0.42645636
## 19 1994 0.52444695
## 20 1996 0.38358640
## 21 1998 -0.11685582
## 22 2000 -0.17850767
## 23 2002 -0.02298488
## 24 2004 -0.54928191
## 25 2006 -0.37904977
## 26 2008 -0.17243380
## 27 2010 -0.04434274
## 28 2012 -0.08974328
## 29 2014 -0.05702656
## 30 2016 -0.04499413Its better to see in plot:
plot(yearlymeans)
**What happened in the previous decade that incited abortion attitudes? Need to find in the literature.
Here is the actual multimodel:
multilevel<-lmer(abscald~bgen*sex*educ+(1|YEAR),
data=GSS, na.action=na.omit)
arm::display(multilevel, detail=T)## lmer(formula = abscald ~ bgen * sex * educ + (1 | YEAR), data = GSS,
## na.action = na.omit)
## coef.est coef.se t value
## (Intercept) -2.42 0.27 -8.92
## bgenBoomers 2.12 0.31 6.76
## bgenGeneration X 1.73 0.39 4.41
## bgenMillennials 2.57 0.54 4.72
## bgenSilent Gen 1.58 0.34 4.60
## sexFemale -0.97 0.32 -3.03
## educ1 4.33 0.98 4.41
## bgenBoomers:sexFemale 0.85 0.40 2.15
## bgenGeneration X:sexFemale 1.41 0.48 2.93
## bgenMillennials:sexFemale -0.18 0.67 -0.27
## bgenSilent Gen:sexFemale 0.84 0.44 1.91
## bgenBoomers:educ1 -3.80 1.06 -3.58
## bgenGeneration X:educ1 -2.13 1.20 -1.78
## bgenMillennials:educ1 -3.23 1.53 -2.11
## bgenSilent Gen:educ1 -2.39 1.15 -2.08
## sexFemale:educ1 -3.69 1.31 -2.81
## bgenBoomers:sexFemale:educ1 4.84 1.41 3.44
## bgenGeneration X:sexFemale:educ1 2.63 1.55 1.70
## bgenMillennials:sexFemale:educ1 5.51 1.96 2.81
## bgenSilent Gen:sexFemale:educ1 3.82 1.53 2.50
##
## Error terms:
## Groups Name Std.Dev.
## YEAR (Intercept) 0.47
## Residual 4.55
## ---
## number of obs: 5262, groups: YEAR, 30
## AIC = 30933.2, DIC = 30886.4
## deviance = 30887.8Variation is significantly represented across all predictors with exception of select marriage categories (seperated and widowed)
summary(multilevel)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: abscald ~ bgen * sex * educ + (1 | YEAR)
## Data: GSS
##
## REML criterion at convergence: 30889.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.56579 -0.71240 0.04303 0.87171 1.90000
##
## Random effects:
## Groups Name Variance Std.Dev.
## YEAR (Intercept) 0.2191 0.4681
## Residual 20.7005 4.5498
## Number of obs: 5262, groups: YEAR, 30
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -2.4198 0.2711 601.1184 -8.924
## bgenBoomers 2.1237 0.3144 4781.4751 6.756
## bgenGeneration X 1.7323 0.3927 3540.0643 4.412
## bgenMillennials 2.5695 0.5444 3536.9851 4.720
## bgenSilent Gen 1.5826 0.3439 5217.2411 4.602
## sexFemale -0.9745 0.3215 5240.5656 -3.032
## educ1 4.3325 0.9835 5231.9370 4.405
## bgenBoomers:sexFemale 0.8492 0.3954 5231.2617 2.148
## bgenGeneration X:sexFemale 1.4149 0.4831 5232.5020 2.929
## bgenMillennials:sexFemale -0.1813 0.6691 5232.5093 -0.271
## bgenSilent Gen:sexFemale 0.8413 0.4414 5230.4377 1.906
## bgenBoomers:educ1 -3.8005 1.0630 5231.4406 -3.575
## bgenGeneration X:educ1 -2.1348 1.1976 5233.5217 -1.783
## bgenMillennials:educ1 -3.2271 1.5263 5222.6015 -2.114
## bgenSilent Gen:educ1 -2.3896 1.1513 5226.9021 -2.076
## sexFemale:educ1 -3.6860 1.3131 5229.8002 -2.807
## bgenBoomers:sexFemale:educ1 4.8417 1.4074 5232.4984 3.440
## bgenGeneration X:sexFemale:educ1 2.6310 1.5492 5229.1208 1.698
## bgenMillennials:sexFemale:educ1 5.5084 1.9625 5227.6188 2.807
## bgenSilent Gen:sexFemale:educ1 3.8183 1.5260 5229.2263 2.502
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## bgenBoomers 1.59e-11 ***
## bgenGeneration X 1.06e-05 ***
## bgenMillennials 2.45e-06 ***
## bgenSilent Gen 4.28e-06 ***
## sexFemale 0.002445 **
## educ1 1.08e-05 ***
## bgenBoomers:sexFemale 0.031799 *
## bgenGeneration X:sexFemale 0.003418 **
## bgenMillennials:sexFemale 0.786423
## bgenSilent Gen:sexFemale 0.056713 .
## bgenBoomers:educ1 0.000353 ***
## bgenGeneration X:educ1 0.074708 .
## bgenMillennials:educ1 0.034539 *
## bgenSilent Gen:educ1 0.037986 *
## sexFemale:educ1 0.005016 **
## bgenBoomers:sexFemale:educ1 0.000586 ***
## bgenGeneration X:sexFemale:educ1 0.089513 .
## bgenMillennials:sexFemale:educ1 0.005022 **
## bgenSilent Gen:sexFemale:educ1 0.012372 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 20 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itarm::display(multilevel,detail=T)## lmer(formula = abscald ~ bgen * sex * educ + (1 | YEAR), data = GSS,
## na.action = na.omit)
## coef.est coef.se t value
## (Intercept) -2.42 0.27 -8.92
## bgenBoomers 2.12 0.31 6.76
## bgenGeneration X 1.73 0.39 4.41
## bgenMillennials 2.57 0.54 4.72
## bgenSilent Gen 1.58 0.34 4.60
## sexFemale -0.97 0.32 -3.03
## educ1 4.33 0.98 4.41
## bgenBoomers:sexFemale 0.85 0.40 2.15
## bgenGeneration X:sexFemale 1.41 0.48 2.93
## bgenMillennials:sexFemale -0.18 0.67 -0.27
## bgenSilent Gen:sexFemale 0.84 0.44 1.91
## bgenBoomers:educ1 -3.80 1.06 -3.58
## bgenGeneration X:educ1 -2.13 1.20 -1.78
## bgenMillennials:educ1 -3.23 1.53 -2.11
## bgenSilent Gen:educ1 -2.39 1.15 -2.08
## sexFemale:educ1 -3.69 1.31 -2.81
## bgenBoomers:sexFemale:educ1 4.84 1.41 3.44
## bgenGeneration X:sexFemale:educ1 2.63 1.55 1.70
## bgenMillennials:sexFemale:educ1 5.51 1.96 2.81
## bgenSilent Gen:sexFemale:educ1 3.82 1.53 2.50
##
## Error terms:
## Groups Name Std.Dev.
## YEAR (Intercept) 0.47
## Residual 4.55
## ---
## number of obs: 5262, groups: YEAR, 30
## AIC = 30933.2, DIC = 30886.4
## deviance = 30887.8#oddsratios
exp(fixef(multilevel)[-1])## bgenBoomers bgenGeneration X
## 8.36194723 5.65352246
## bgenMillennials bgenSilent Gen
## 13.05951522 4.86744531
## sexFemale educ1
## 0.37738163 76.13134828
## bgenBoomers:sexFemale bgenGeneration X:sexFemale
## 2.33775923 4.11624648
## bgenMillennials:sexFemale bgenSilent Gen:sexFemale
## 0.83418732 2.31943837
## bgenBoomers:educ1 bgenGeneration X:educ1
## 0.02236040 0.11826623
## bgenMillennials:educ1 bgenSilent Gen:educ1
## 0.03967324 0.09166487
## sexFemale:educ1 bgenBoomers:sexFemale:educ1
## 0.02507297 126.68298574
## bgenGeneration X:sexFemale:educ1 bgenMillennials:sexFemale:educ1
## 13.88765044 246.74772786
## bgenSilent Gen:sexFemale:educ1
## 45.52724962#CIs
exp(confint(multilevel,method="Wald"))## 2.5 % 97.5 %
## .sig01 NA NA
## .sigma NA NA
## (Intercept) 0.052275510 1.513208e-01
## bgenBoomers 4.515708699 1.548421e+01
## bgenGeneration X 2.618737000 1.220524e+01
## bgenMillennials 4.492635651 3.796233e+01
## bgenSilent Gen 2.480784301 9.550215e+00
## sexFemale 0.200981751 7.086061e-01
## educ1 11.077651772 5.232140e+02
## bgenBoomers:sexFemale 1.076979773 5.074485e+00
## bgenGeneration X:sexFemale 1.596851783 1.061056e+01
## bgenMillennials:sexFemale 0.224779415 3.095784e+00
## bgenSilent Gen:sexFemale 0.976435509 5.509626e+00
## bgenBoomers:educ1 0.002783671 1.796144e-01
## bgenGeneration X:educ1 0.011310386 1.236642e+00
## bgenMillennials:educ1 0.001991970 7.901554e-01
## bgenSilent Gen:educ1 0.009598197 8.754195e-01
## sexFemale:educ1 0.001912204 3.287587e-01
## bgenBoomers:sexFemale:educ1 8.030963399 1.998338e+03
## bgenGeneration X:sexFemale:educ1 0.666715125 2.892792e+02
## bgenMillennials:sexFemale:educ1 5.269334784 1.155448e+04
## bgenSilent Gen:sexFemale:educ1 2.287488158 9.061164e+02library(sjPlot)## Learn more about sjPlot with 'browseVignettes("sjPlot")'.plot_model(multilevel,facet.grid=F)