Week III - Logit/Probit Models II
R. H. Lowe Jr., M.S. | February 18, 2019 | Dem 7283 | University of Texas at San Antonio
Research Purpose
Before running the models, all required packages are installed:
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(questionr)
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, unionAll variables are recoded accordingly, starting with the DV (pregnancy intentions). Weights (wmweight) and PSU (HH1) are included prior to telling R the survey design:
library(haven)
wm <- read_sav("Desktop/Data/Ghana MICS 2011 SPSS Datasets/wm.sav")
#Pregnancy Intentions: "When you got pregnant, did you want to get pregnant at that time?"
wm$pregintnt<-Recode(wm$DB1,recodes="1=0;2=1;9=NA")
#Age: "How old are you?"
wm$agec<-Recode(wm$WB2CA,recodes="15:19='15-19';20:29='20-29';30:39='30-39';40:49='40-49';else=NA",as.factor=T)
wm$agec<-relevel(wm$agec,ref='40-49')
wm$teen<-Recode(wm$WB2CA,recodes="15:19=1;9=NA;else=0")
wm$yadult<-Recode(wm$WB2CA,recodes="20:29=1;9=NA;else=0")
wm$adult<-Recode(wm$WB2CA,recodes="30:39=1;9=NA;else=0")
wm$midage<-Recode(wm$WB2CA,recodes="40:49=1;9=NA;else=0")
#Educational Attainment:"What is the highest level of school you attended?"
wm$educ<-Recode(wm$WB4,recodes="0='1Pre';1='2Prim';2='3Sec';3='4High';4:6=NA", as.factor=T)
wm$educ<-relevel(wm$educ,ref='4High')
#Contraceptive Use:"Are you currently doing something or using any method to delay or avoid getting pregnant?"
wm$cntrcpt<-Recode(wm$CP2,recodes="1='1Cntrcpt';2='2NoCntrcpt';8:9=NA",as.factor=T)
#Marriage:"Are you currently married or living together with a man as if married?"
wm$marrg<-Recode(wm$MA1,recodes="1='1Marrd';2='2Cohab';3='3Sngl';9=NA",as.factor=T)
sub<-wm%>%
select(pregintnt,agec,educ,cntrcpt,marrg,wmweight,teen,yadult,adult,midage,HH1) %>%
filter( complete.cases(.))
options(survey.lonely.psu = "adjust")
des<-svydesign(ids=~HH1,weights=~wmweight, data=sub)The percentage of unintended preganancies by age is plotted and a survey-corrected chi-square test for independence is performed. Results indicate a signficant relationship between unintentional pregnancy and age, F (2.92, 1802.2) = 7.08), p < .05
library(ggplot2)
cat<-svyby(formula = ~pregintnt, by = ~agec, design = des, FUN = svymean, na.rm=T)
svychisq(~pregintnt+agec, design = des)##
## Pearson's X^2: Rao & Scott adjustment
##
## data: svychisq(~pregintnt + agec, design = des)
## F = 7.0848, ndf = 2.9208, ddf = 1802.2000, p-value = 0.0001169qplot(x=cat$agec,y=cat$pregintnt, data=cat ,xlab="Age", ylab="% Pregnancy Intentions" )+
geom_errorbar(aes(x=agec, ymin=pregintnt-1.96*se,ymax= pregintnt+1.96*se), width=.25)+
ggtitle(label = "% of Ghanaian Women with Unintended Pregnancies by Age")
The percentage of unintended preganancies by educational attainment is plotted and another survey-corrected chi-square test for independence performed. Results indicate that there is a significant relationship between unintended pregnancies and educational attainment, F (2.55, 1575) = 3.0584), p < .05
library(ggplot2)
cat<-svyby(formula = ~pregintnt, by = ~educ, design = des, FUN = svymean, na.rm=T)
svychisq(~pregintnt+educ, design = des)##
## Pearson's X^2: Rao & Scott adjustment
##
## data: svychisq(~pregintnt + educ, design = des)
## F = 3.0584, ndf = 2.5528, ddf = 1575.0000, p-value = 0.03489qplot(x=cat$educ,y=cat$pregintnt, data=cat ,xlab="Education", ylab="% Pregnancy Intentions" )+
geom_errorbar(aes(x=educ, ymin=pregintnt-1.96*se,ymax= pregintnt+1.96*se), width=.25)+
ggtitle(label = "% of Ghanaian Women with Unintended Pregnancies by Education")
Another cross tabulation is performed with both covariates:
catdog<-svyby(formula = ~pregintnt, by = ~agec+educ, design = des, FUN = svymean, na.rm=T)
catdog## agec educ pregintnt se
## 40-49.4High 40-49 4High 0.0000000 0.00000000
## 15-19.4High 15-19 4High 1.0000000 0.00000000
## 20-29.4High 20-29 4High 0.4296271 0.06602270
## 30-39.4High 30-39 4High 0.3164448 0.07577250
## 15-19.1Pre 15-19 1Pre 1.0000000 0.00000000
## 20-29.1Pre 20-29 1Pre 0.6519622 0.28341094
## 30-39.1Pre 30-39 1Pre 0.6443783 0.32433651
## 40-49.2Prim 40-49 2Prim 0.6148152 0.09270836
## 15-19.2Prim 15-19 2Prim 0.7408586 0.08687212
## 20-29.2Prim 20-29 2Prim 0.5062683 0.04490183
## 30-39.2Prim 30-39 2Prim 0.4371690 0.04363962
## 40-49.3Sec 40-49 3Sec 0.5159563 0.10379197
## 15-19.3Sec 15-19 3Sec 0.7577802 0.06695439
## 20-29.3Sec 20-29 3Sec 0.5394866 0.03658701
## 30-39.3Sec 30-39 3Sec 0.4763527 0.04716227Here, I run Logit & Probit Models, both of which include all covariates in the intial survey design. Odds Ratios and Confidence intervals are performed as well:
fit.logit<-svyglm(pregintnt~agec+educ+cntrcpt+marrg,
design= des,
family=binomial)## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!exp(coef(fit.logit))#odds ratios## (Intercept) agec15-19 agec20-29 agec30-39
## 0.7557847 1.6217066 0.7188254 0.6664974
## educ1Pre educ2Prim educ3Sec cntrcpt2NoCntrcpt
## 4.7259128 1.3710198 1.5127027 0.7160433
## marrg2Cohab marrg3Sngl
## 2.3729842 2.7594222exp(confint(fit.logit)) #confidence intervals for odds ratios## 2.5 % 97.5 %
## (Intercept) 0.3604891 1.5845432
## agec15-19 0.7077647 3.7158284
## agec20-29 0.3781836 1.3662939
## agec30-39 0.3530257 1.2583188
## educ1Pre 0.9304392 24.0039886
## educ2Prim 0.8193033 2.2942605
## educ3Sec 0.9293711 2.4621698
## cntrcpt2NoCntrcpt 0.5165317 0.9926168
## marrg2Cohab 1.6979681 3.3163485
## marrg3Sngl 1.7064306 4.4621860#probit model
fit.probit<-svyglm(pregintnt~agec+educ+cntrcpt+marrg,
design=des,
family=binomial(link= "probit"))## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!exp(confint(fit.probit)) #confidence intervals for odds ratios## 2.5 % 97.5 %
## (Intercept) 0.5291179 1.3232040
## agec15-19 0.8127579 2.2292783
## agec20-29 0.5496716 1.2194245
## agec30-39 0.5265041 1.1583438
## educ1Pre 1.0459086 6.7854585
## educ2Prim 0.8853588 1.6638537
## educ3Sec 0.9562822 1.7383989
## cntrcpt2NoCntrcpt 0.6662039 0.9956586
## marrg2Cohab 1.3921459 2.1050337
## marrg3Sngl 1.3965949 2.5231537myexp<-function(x){exp(x)}
stargazer(fit.logit,
fit.probit,
type = "text",
style="demography",
covariate.labels=c("Age 15-19" ,"Age 20-29","Age 30-39","Preschool", "Primary School", "Secondary Education","Contraceptive Use", "Living with Male, Unmarried", "Single"),
ci=T, column.sep.width = "10pt",apply.coef = myexp)##
## -----------------------------------------------------------
## pregintnt
## survey-weighted survey-weighted
## logistic probit
## Model 1 Model 2
## -----------------------------------------------------------
## Age 15-19 1.622*** 1.346***
## (0.793, 2.451) (0.842, 1.851)
## Age 20-29 0.719* 0.819***
## (0.077, 1.361) (0.420, 1.217)
## Age 30-39 0.666* 0.781***
## (0.031, 1.302) (0.387, 1.175)
## Preschool 4.726*** 2.664***
## (3.101, 6.351) (1.729, 3.599)
## Primary School 1.371*** 1.214***
## (0.856, 1.886) (0.898, 1.529)
## Secondary Education 1.513*** 1.289***
## (1.026, 2.000) (0.991, 1.588)
## Contraceptive Use 0.716*** 0.814***
## (0.389, 1.043) (0.614, 1.015)
## Living with Male, Unmarried 2.373*** 1.712***
## (2.038, 2.708) (1.505, 1.919)
## Single 2.759*** 1.877***
## (2.279, 3.240) (1.581, 2.173)
## Constant 0.756* 0.837***
## (0.015, 1.496) (0.378, 1.295)
## N 1,291 1,291
## Log Likelihood -786.262 -786.233
## AIC 1,592.525 1,592.466
## -----------------------------------------------------------
## *p < .05; **p < .01; ***p < .001The models mutually suggest that women aged 15-19 had an increased chance of reporting an unintended pregnancy compared to their middle aged counterparts (reference group = Age 40-49). Similarly, Ghanaian women with just a primary school education had an increased chance in reporting an unintended pregnancy as well. Interestingly, the probit model has stronger significance levels than the logit model across age cagtegories and the constant. Next, I retrieve marginal effects (and plot them) to determine the rate of change between the dependent and independent variables.
#Logit marginal effects
log.marg<-coef(fit.logit)*mean(dlogis(predict(fit.logit)), na.rm=T)
#Progit marginal effects
prob.marg<-coef(fit.probit)*mean(dnorm(predict(fit.probit)), na.rm=T)
effects<-data.frame(effect=c(log.marg, prob.marg),
term=rep(names(log.marg),2),
model=c(rep("logit", length(log.marg)),rep("probit", length(log.marg)))
)
effects%>%
ggplot(aes(x=term, y=effect))+geom_point(aes(color=model, group=model, shape=model),
position=position_jitterdodge(jitter.width = 1),
size=2)+
ylab("Marginal Effect")+
xlab("Model Term")+
geom_abline(intercept = 0, slope=0)+
theme(axis.text.x = element_text(angle = 45, hjust = 1))+
ggtitle(label = "Comparison of marginal effects in Logit and Probit models")
We see that marginal effects are comparable between the logit and probit model, though large distinctions exist with the first two age categories of 15-19 and 20-29.
#get a series of predicted probabilites for different "types" of people for each model
#expand.grid will generate all possible combinations of values you specify
dat<-expand.grid(educ=levels(wm$educ), agec=levels(wm$agec),marrg=levels(wm$marrg),cntrcpt=levels(wm$cntrcpt))
#You MAY need to get rid of impossible cases here
#generate the fitted values
fit<-predict(fit.logit, newdat=dat,type="response")
fitp<-predict(fit.probit, newdat=dat,type="response")
#add the values to the fake data
dat$fitted.prob.lrm<-round(fit, 3)
dat$fitted.prob.pro<-round(fitp, 3)
#Print the fitted probabilities for the first 20 cases
knitr::kable(head(dat, n=20))| educ | agec | marrg | cntrcpt | fitted.prob.lrm | fitted.prob.pro |
|---|---|---|---|---|---|
| 4High | 40-49 | 1Marrd | 1Cntrcpt | 0.430 | 0.429 |
| 1Pre | 40-49 | 1Marrd | 1Cntrcpt | 0.781 | 0.789 |
| 2Prim | 40-49 | 1Marrd | 1Cntrcpt | 0.509 | 0.506 |
| 3Sec | 40-49 | 1Marrd | 1Cntrcpt | 0.533 | 0.530 |
| 4High | 15-19 | 1Marrd | 1Cntrcpt | 0.551 | 0.547 |
| 1Pre | 15-19 | 1Marrd | 1Cntrcpt | 0.853 | 0.864 |
| 2Prim | 15-19 | 1Marrd | 1Cntrcpt | 0.627 | 0.623 |
| 3Sec | 15-19 | 1Marrd | 1Cntrcpt | 0.650 | 0.645 |
| 4High | 20-29 | 1Marrd | 1Cntrcpt | 0.352 | 0.353 |
| 1Pre | 20-29 | 1Marrd | 1Cntrcpt | 0.720 | 0.726 |
| 2Prim | 20-29 | 1Marrd | 1Cntrcpt | 0.427 | 0.427 |
| 3Sec | 20-29 | 1Marrd | 1Cntrcpt | 0.451 | 0.451 |
| 4High | 30-39 | 1Marrd | 1Cntrcpt | 0.335 | 0.335 |
| 1Pre | 30-39 | 1Marrd | 1Cntrcpt | 0.704 | 0.710 |
| 2Prim | 30-39 | 1Marrd | 1Cntrcpt | 0.409 | 0.408 |
| 3Sec | 30-39 | 1Marrd | 1Cntrcpt | 0.432 | 0.432 |
| 4High | 40-49 | 2Cohab | 1Cntrcpt | 0.642 | 0.640 |
| 1Pre | 40-49 | 2Cohab | 1Cntrcpt | 0.894 | 0.910 |
| 2Prim | 40-49 | 2Cohab | 1Cntrcpt | 0.711 | 0.710 |
| 3Sec | 40-49 | 2Cohab | 1Cntrcpt | 0.731 | 0.730 |
There are several alarming findings in this analysis. Specifically, I specify the probability of reporting an unintended pregnancy for pre-school educated, cohabitiating woman aged between 40-49 years, compared to a highly educated woman aged between 15-19.
dat[which(dat$educ=="1Pre"&dat$agec=="40-49"&dat$marrg=="2Cohab"),]## educ agec marrg cntrcpt fitted.prob.lrm fitted.prob.pro
## 18 1Pre 40-49 2Cohab 1Cntrcpt 0.894 0.910
## 66 1Pre 40-49 2Cohab 2NoCntrcpt 0.859 0.872dat[which(dat$educ=="4High"&dat$agec=="30-39"&dat$marrg=="1Marrd"),]## educ agec marrg cntrcpt fitted.prob.lrm fitted.prob.pro
## 13 4High 30-39 1Marrd 1Cntrcpt 0.335 0.335
## 61 4High 30-39 1Marrd 2NoCntrcpt 0.265 0.264Here, we see that highly educated Ghanaian women in union and aged between 30-39 have only a 33.5 percent probability of reporting an unintentional pregnancy; contrast to pre-school educated Ghanaian women, unmarried but cohabiting with male, and aged between 40-49 - who have an estimated probability of reporting unintentional pregancy of about 89 percent. Women with higher probabilities however are mostly using some sort of contraceptive, while those who have less likelihood in reporting unintential pregnancy are not currently using a contraceptive, which could potentially put them at risk.
Homework 2, Part II (for week of 2.18.19)
Nested Model Comparison
fit.logit1<-svyglm(pregintnt~agec, design = des, family=binomial) #age only## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!fit.logit2<-svyglm(pregintnt~agec+educ, design = des, family=binomial) #age+education## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!fit.logit3<-svyglm(pregintnt~agec+educ+marrg, design = des, family=binomial) #age+education+marriage+marriage## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!summary(fit.logit1)##
## Call:
## svyglm(formula = pregintnt ~ agec, design = des, family = binomial)
##
## Survey design:
## svydesign(ids = ~HH1, weights = ~wmweight, data = sub)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.1957 0.2743 0.714 0.4757
## agec15-19 0.9523 0.3881 2.454 0.0144 *
## agec20-29 -0.1422 0.3002 -0.474 0.6359
## agec30-39 -0.4238 0.3056 -1.387 0.1661
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1.000775)
##
## Number of Fisher Scoring iterations: 4summary(fit.logit2)##
## Call:
## svyglm(formula = pregintnt ~ agec + educ, design = des, family = binomial)
##
## Survey design:
## svydesign(ids = ~HH1, weights = ~wmweight, data = sub)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.26751 0.35394 -0.756 0.4501
## agec15-19 0.93690 0.38919 2.407 0.0164 *
## agec20-29 -0.09316 0.30058 -0.310 0.7567
## agec30-39 -0.38590 0.30502 -1.265 0.2063
## educ1Pre 1.42373 0.88086 1.616 0.1065
## educ2Prim 0.42295 0.24968 1.694 0.0908 .
## educ3Sec 0.51803 0.23488 2.206 0.0278 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1.000577)
##
## Number of Fisher Scoring iterations: 4summary(fit.logit3)##
## Call:
## svyglm(formula = pregintnt ~ agec + educ + marrg, design = des,
## family = binomial)
##
## Survey design:
## svydesign(ids = ~HH1, weights = ~wmweight, data = sub)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.4806 0.3718 -1.293 0.1966
## agec15-19 0.4533 0.4149 1.093 0.2750
## agec20-29 -0.2988 0.3250 -0.919 0.3582
## agec30-39 -0.3745 0.3192 -1.173 0.2412
## educ1Pre 1.5597 0.8848 1.763 0.0784 .
## educ2Prim 0.2999 0.2602 1.153 0.2496
## educ3Sec 0.4074 0.2465 1.653 0.0989 .
## marrg2Cohab 0.8513 0.1704 4.996 7.65e-07 ***
## marrg3Sngl 0.9521 0.2383 3.995 7.25e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 0.999572)
##
## Number of Fisher Scoring iterations: 4The reduction of differences across the models are apparent when additonal covariates are added. For instance, the only significant predictor of unintended pregnancies in model 1 pertained to women aged between 15 and 19. In Model 2, we see that secondary education becomes significant. The last model completely diminishes the significance of the age 15-19 variable, as marriage becomes the stongest significant predictors.
regTermTest(fit.logit3,test.terms=~agec+educ, method="Wald", df=NULL)## Wald test for agec educ
## in svyglm(formula = pregintnt ~ agec + educ + marrg, design = des,
## family = binomial)
## F = 2.263655 on 6 and 609 df: p= 0.036061stargazer(fit.logit1, fit.logit2, fit.logit3,type = "text", style="demography", covariate.labels =c("age 15-19","age 20-29" ,"age 30-39","Preschool", "Primary School", "Secondary School", "Cohabiting","Single"), ci = T )##
## ----------------------------------------------------------------
## pregintnt
## Model 1 Model 2 Model 3
## ----------------------------------------------------------------
## age 15-19 0.952* 0.937* 0.453
## (0.192, 1.713) (0.174, 1.700) (-0.360, 1.266)
## age 20-29 -0.142 -0.093 -0.299
## (-0.731, 0.446) (-0.682, 0.496) (-0.936, 0.338)
## age 30-39 -0.424 -0.386 -0.375
## (-1.023, 0.175) (-0.984, 0.212) (-1.000, 0.251)
## Preschool 1.424 1.560
## (-0.303, 3.150) (-0.175, 3.294)
## Primary School 0.423 0.300
## (-0.066, 0.912) (-0.210, 0.810)
## Secondary School 0.518* 0.407
## (0.058, 0.978) (-0.076, 0.891)
## Cohabiting 0.851***
## (0.517, 1.185)
## Single 0.952***
## (0.485, 1.419)
## Constant 0.196 -0.268 -0.481
## (-0.342, 0.733) (-0.961, 0.426) (-1.209, 0.248)
## N 1,291 1,291 1,291
## Log Likelihood -817.327 -813.303 -789.118
## AIC 1,642.653 1,640.605 1,596.235
## ----------------------------------------------------------------
## *p < .05; **p < .01; ***p < .001Stratified Models
fit.logitint<-svyglm(pregintnt~agec*educ+marrg, design = des, family=binomial)#age*education interaction+marriage## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!regTermTest(fit.logitint, test.terms = ~agec:educ, method = "Wald", df=NULL)## Wald test for agec:educ
## in svyglm(formula = pregintnt ~ agec * educ + marrg, design = des,
## family = binomial)
## F = 97.48589 on 8 and 601 df: p= < 2.22e-16The wald test proves there is a signficant interaction between age and education. Now, we can stratify:
fit.unrestricted<-svyglm(pregintnt~agec+educ+marrg,design= des, family=binomial)## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!fit.logit.teen<-svyglm(pregintnt~(educ+marrg),design= subset(des, teen==1), family=binomial)## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!fit.logit.yadult<-svyglm(pregintnt~(educ+marrg),design= subset(des, yadult==1), family=binomial)## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!fit.logit.adult<-svyglm(pregintnt~(educ+marrg),design= subset(des, adult==1), family=binomial)## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!fit.logit.midage<-svyglm(pregintnt~(educ+marrg),design= subset(des, midage==1), family=binomial)## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!stargazer(fit.logit.teen, fit.logit.yadult, fit.logit.adult, fit.logit.midage,
type = "text", style="demography",
covariate.labels =c( "Preschool","Primary School", "Secondary School", "Cohabiting", "Single", "Intercept"),
column.labels = c("Teen 15-19", "Young Adult 20-29", "Adult 30-39", "Middle Age 40-49"),
ci=T )##
## -----------------------------------------------------------------------------------------
## pregintnt
## Teen 15-19 Young Adult 20-29 Adult 30-39 Middle Age 40-49
## Model 1 Model 2 Model 3 Model 4
## -----------------------------------------------------------------------------------------
## Preschool 0.137 1.306 1.448
## (-2.087, 2.361) (-1.231, 3.843) (-1.414, 4.311)
## Primary School -15.371*** 0.319 0.319 16.355***
## (-16.867, -13.875) (-0.388, 1.025) (-0.469, 1.107) (14.215, 18.495)
## Secondary School -15.178*** 0.347 0.529 15.959***
## (-16.382, -13.975) (-0.276, 0.970) (-0.249, 1.307) (13.831, 18.088)
## Cohabiting 1.794* 1.048*** 0.701* 0.032
## (0.197, 3.391) (0.585, 1.512) (0.094, 1.308) (-1.195, 1.258)
## Single 1.943* 1.169*** 0.406 -0.074
## (0.425, 3.460) (0.561, 1.776) (-0.581, 1.393) (-2.508, 2.361)
## Intercept 14.714*** -0.871** -0.854* -15.896***
## (13.300, 16.129) (-1.516, -0.225) (-1.556, -0.152) (-17.856, -13.936)
## N 124 620 460 87
## Log Likelihood -54.920 -384.970 -288.434 -54.216
## AIC 121.839 781.940 588.867 118.432
## -----------------------------------------------------------------------------------------
## *p < .05; **p < .01; ***p < .001beta.test<-function(model1, model2, betaname){
s1<-summary(model1)$coef
s2<-summary(model2)$coef
db <- ((s2[rownames(s2)==betaname,1]-s1[rownames(s1)==betaname,1]))^2
sd <-s2[rownames(s2)==betaname,2]^2+s1[rownames(s1)==betaname,2]^2
td <- db/sd
beta1=s1[rownames(s1)==betaname,1]
beta2=s2[rownames(s2)==betaname,1]
pv<-1- pchisq(td, df = 1)
print(list(beta=betaname,beta1=beta1, beta2=beta2, x2=td, pvalue=pv))
}
beta.test(fit.logit.teen, fit.logit.midage, "educ3Sec")## $beta
## [1] "educ3Sec"
##
## $beta1
## [1] -15.17834
##
## $beta2
## [1] 15.95941
##
## $x2
## [1] 622.7513
##
## $pvalue
## [1] 0beta.test(fit.logit.teen, fit.logit.midage, "marrg3Sngl")## $beta
## [1] "marrg3Sngl"
##
## $beta1
## [1] 1.942613
##
## $beta2
## [1] -0.07352654
##
## $x2
## [1] 1.897979
##
## $pvalue
## [1] 0.1683047We see that the effect of secondary education is the same between the two specified groups. There is no mutal effect in reference to singlehood however.