library(car)
## Loading required package: carData
library(mice)
## Warning: package 'mice' was built under R version 3.5.2
## Loading required package: lattice
##
## Attaching package: 'mice'
## The following objects are masked from 'package:base':
##
## cbind, rbind
library(ggplot2)
library(dplyr)
## Warning: package 'dplyr' was built under R version 3.5.2
##
## 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, union
library(foreign)
brfss2017=read.xport("~/Desktop/DataFolder/brfss.XPT ")
samp<-sample(1:dim(brfss2017)[1],size=25000)
brfss2017<-brfss2017[samp,]
#Recode variables of interest
brfss2017$menth<-Recode(brfss2017$MENTHLTH,recodes="88=0;77=NA;99=NA",as.factor=T)
brfss2017$hlthinsur<-Recode(brfss2017$HLTHPLN1,recodes="1='Y';2='N';else=NA",as.factor=T)
brfss2017$medcost<-Recode(brfss2017$MEDCOST,recodes="1='Y';2='N';else=NA",as.factor=T)
brfss2017$sxorient<-Recode(brfss2017$SXORIENT,recodes="1='S';2='G';else=NA",as.factor=T)
brfss2017$marital<-Recode(brfss2017$MARITAL, recodes="1='M';5='S';else=NA",as.factor=T)
#Make a smaller sample size
summary(brfss2017[,c("menth","hlthinsur","medcost","sxorient","marital")])
## menth hlthinsur medcost sxorient marital
## 0 :16656 N : 1974 N :22268 G : 159 M :12982
## 30 : 1359 Y :22914 Y : 2657 S :10609 S : 4062
## 2 : 1196 NA's: 112 NA's: 75 NA's:14232 NA's: 7956
## 5 : 921
## 1 : 784
## (Other): 3708
## NA's : 376
#Look at patterns of missingness
md.pattern(brfss2017[,c("menth","hlthinsur","medcost","sxorient","marital")])
## medcost hlthinsur menth marital sxorient
## 7169 1 1 1 1 1 0
## 9577 1 1 1 1 0 1
## 3403 1 1 1 0 1 1
## 4306 1 1 1 0 0 2
## 82 1 1 0 1 1 1
## 108 1 1 0 1 0 2
## 56 1 1 0 0 1 2
## 118 1 1 0 0 0 3
## 22 1 0 1 1 1 1
## 44 1 0 1 1 0 2
## 11 1 0 1 0 1 2
## 23 1 0 1 0 0 3
## 1 1 0 0 1 1 2
## 4 1 0 0 1 0 3
## 1 1 0 0 0 1 3
## 10 0 1 1 1 1 1
## 22 0 1 1 1 0 2
## 12 0 1 1 0 1 2
## 22 0 1 1 0 0 3
## 3 0 1 0 1 0 3
## 1 0 0 1 1 1 2
## 2 0 0 1 0 0 4
## 1 0 0 0 1 0 4
## 2 0 0 0 0 0 5
## 75 112 376 7956 14232 22751
#Let's see how pairs of variables are missing together...
md.pairs(brfss2017[,c("menth","hlthinsur","medcost","sxorient","marital")])
## $rr
## menth hlthinsur medcost sxorient marital
## menth 24624 24521 24555 10628 16845
## hlthinsur 24521 24888 24819 10732 16971
## medcost 24555 24819 24925 10745 17007
## sxorient 10628 10732 10745 10768 7285
## marital 16845 16971 17007 7285 17044
##
## $rm
## menth hlthinsur medcost sxorient marital
## menth 0 103 69 13996 7779
## hlthinsur 367 0 69 14156 7917
## medcost 370 106 0 14180 7918
## sxorient 140 36 23 0 3483
## marital 199 73 37 9759 0
##
## $mr
## menth hlthinsur medcost sxorient marital
## menth 0 367 370 140 199
## hlthinsur 103 0 106 36 73
## medcost 69 69 0 23 37
## sxorient 13996 14156 14180 0 9759
## marital 7779 7917 7918 3483 0
##
## $mm
## menth hlthinsur medcost sxorient marital
## menth 376 9 6 236 177
## hlthinsur 9 112 6 76 39
## medcost 6 6 75 52 38
## sxorient 236 76 52 14232 4473
## marital 177 39 38 4473 7956
#find the most common value for sexual orientation
mcv.sxorient<-factor(names(which.max(table(brfss2017$sxorient))),levels=levels(brfss2017$sxorient))
mcv.sxorient
## [1] S
## Levels: G S
#find the most common value for marital status
mcv.marital<-factor(names(which.max(table(brfss2017$marital))),levels=levels(brfss2017$marital))
mcv.marital
## [1] M
## Levels: M S
#find the most common value for outstanding medical bills
mcv.medcost<-factor(names(which.max(table(brfss2017$medcost))),levels=levels(brfss2017$medcost))
mcv.medcost
## [1] N
## Levels: N Y
#find the most common value for having health insurance
mcv.hlthinsur<-factor(names(which.max(table(brfss2017$hlthinsur))),levels=levels(brfss2017$hlthinsur))
mcv.hlthinsur
## [1] Y
## Levels: N Y
#find the most common value for self reported mental health
mcv.menth<-factor(names(which.max(table(brfss2017$menth))),levels=levels(brfss2017$menth))
mcv.menth
## [1] 0
## 31 Levels: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 ... 30
#impute the cases for sexual orientation
brfss2017$sxorient.imp<-as.factor(ifelse(is.na(brfss2017$sxorient)==T, mcv.sxorient, brfss2017$sxorient))
levels(brfss2017$sxorient.imp)<-levels(brfss2017$sxorient)
prop.table(table(brfss2017$sxorient))
##
## G S
## 0.01476597 0.98523403
prop.table(table(brfss2017$sxorient.imp))
##
## G S
## 0.00636 0.99364
barplot(prop.table(table(brfss2017$sxorient)), main="Original Data", ylim=c(0, .6))
barplot(prop.table(table(brfss2017$sxorient.imp)), main="Imputed Data",ylim=c(0, .6))
#impute the cases for marital status
brfss2017$marital.imp<-as.factor(ifelse(is.na(brfss2017$marital)==T, mcv.marital, brfss2017$marital))
levels(brfss2017$marital.imp)<-levels(brfss2017$marital)
prop.table(table(brfss2017$marital))
##
## M S
## 0.7616757 0.2383243
prop.table(table(brfss2017$marital.imp))
##
## M S
## 0.83752 0.16248
barplot(prop.table(table(brfss2017$marital)), main="Original Data", ylim=c(0, .6))
barplot(prop.table(table(brfss2017$marital.imp)), main="Imputed Data",ylim=c(0, .6))
#impute the cases for outstanding medical bills
brfss2017$medcost.imp<-as.factor(ifelse(is.na(brfss2017$medcost)==T, mcv.medcost, brfss2017$medcost))
levels(brfss2017$medcost.imp)<-levels(brfss2017$medcost)
prop.table(table(brfss2017$medcost))
##
## N Y
## 0.8934002 0.1065998
prop.table(table(brfss2017$medcost.imp))
##
## N Y
## 0.89372 0.10628
barplot(prop.table(table(brfss2017$medcost)), main="Original Data", ylim=c(0, .6))
barplot(prop.table(table(brfss2017$medcost.imp)), main="Imputed Data",ylim=c(0, .6))
#impute the cases for having health insurance
brfss2017$hlthinsur.imp<-as.factor(ifelse(is.na(brfss2017$hlthinsur)==T, mcv.hlthinsur, brfss2017$hlthinsur))
levels(brfss2017$hlthinsur.imp)<-levels(brfss2017$hlthinsur)
prop.table(table(brfss2017$hlthinsur))
##
## N Y
## 0.07931533 0.92068467
prop.table(table(brfss2017$hlthinsur.imp))
##
## N Y
## 0.07896 0.92104
barplot(prop.table(table(brfss2017$hlthinsur)), main="Original Data", ylim=c(0, .6))
barplot(prop.table(table(brfss2017$hlthinsur.imp)), main="Imputed Data",ylim=c(0, .6))
#impute the cases for self reported mental health
brfss2017$menth.imp<-as.factor(ifelse(is.na(brfss2017$menth)==T, mcv.menth, brfss2017$menth))
levels(brfss2017$menth.imp)<-levels(brfss2017$menth)
prop.table(table(brfss2017$menth))
##
## 0 1 2 3 4
## 6.764133e-01 3.183886e-02 4.857050e-02 3.086420e-02 1.555393e-02
## 5 6 7 8 9
## 3.740253e-02 3.492528e-03 1.376706e-02 2.558480e-03 5.279402e-04
## 10 11 12 13 14
## 2.534113e-02 4.061079e-05 1.543210e-03 6.091618e-04 4.548408e-03
## 15 16 17 18 19
## 2.481319e-02 4.061079e-04 6.903834e-04 5.685510e-04 4.061079e-05
## 20 21 22 23 24
## 1.441683e-02 6.091618e-04 1.624431e-04 2.030539e-04 2.842755e-04
## 25 26 27 28 29
## 6.457115e-03 2.030539e-04 4.873294e-04 1.624431e-03 7.716049e-04
## 30
## 5.519006e-02
prop.table(table(brfss2017$menth.imp))
##
## 0 1 2 3 4 5 6 7 8
## 0.68128 0.03136 0.04784 0.03040 0.01532 0.03684 0.00344 0.01356 0.00252
## 9 10 11 12 13 14 15 16 17
## 0.00052 0.02496 0.00004 0.00152 0.00060 0.00448 0.02444 0.00040 0.00068
## 18 19 20 21 22 23 24 25 26
## 0.00056 0.00004 0.01420 0.00060 0.00016 0.00020 0.00028 0.00636 0.00020
## 27 28 29 30
## 0.00048 0.00160 0.00076 0.05436
barplot(prop.table(table(brfss2017$menth)), main="Original Data", ylim=c(0, .6))
barplot(prop.table(table(brfss2017$menth.imp)), main="Imputed Data",ylim=c(0, .6))
#Multiple imputation
imp<-mice(data=brfss2017[,c("menth","hlthinsur","medcost","sxorient","marital")],seed=22,m=5)
##
## iter imp variable
## 1 1 menth hlthinsur medcost sxorient marital
## 1 2 menth hlthinsur medcost sxorient marital
## 1 3 menth hlthinsur medcost sxorient marital
## 1 4 menth hlthinsur medcost sxorient marital
## 1 5 menth hlthinsur medcost sxorient marital
## 2 1 menth hlthinsur medcost sxorient marital
## 2 2 menth hlthinsur medcost sxorient marital
## 2 3 menth hlthinsur medcost sxorient marital
## 2 4 menth hlthinsur medcost sxorient marital
## 2 5 menth hlthinsur medcost sxorient marital
## 3 1 menth hlthinsur medcost sxorient marital
## 3 2 menth hlthinsur medcost sxorient marital
## 3 3 menth hlthinsur medcost sxorient marital
## 3 4 menth hlthinsur medcost sxorient marital
## 3 5 menth hlthinsur medcost sxorient marital
## 4 1 menth hlthinsur medcost sxorient marital
## 4 2 menth hlthinsur medcost sxorient marital
## 4 3 menth hlthinsur medcost sxorient marital
## 4 4 menth hlthinsur medcost sxorient marital
## 4 5 menth hlthinsur medcost sxorient marital
## 5 1 menth hlthinsur medcost sxorient marital
## 5 2 menth hlthinsur medcost sxorient marital
## 5 3 menth hlthinsur medcost sxorient marital
## 5 4 menth hlthinsur medcost sxorient marital
## 5 5 menth hlthinsur medcost sxorient marital
## Warning: Number of logged events: 100
print(imp)
## Class: mids
## Number of multiple imputations: 5
## Imputation methods:
## menth hlthinsur medcost sxorient marital
## "polyreg" "logreg" "logreg" "logreg" "logreg"
## PredictorMatrix:
## menth hlthinsur medcost sxorient marital
## menth 0 1 1 1 1
## hlthinsur 1 0 1 1 1
## medcost 1 1 0 1 1
## sxorient 1 1 1 0 1
## marital 1 1 1 1 0
## Number of logged events: 100
## it im dep meth out
## 1 1 1 hlthinsur logreg menth11, menth19
## 2 1 1 medcost logreg menth11, menth19
## 3 1 1 sxorient logreg menth11, menth19, menth22, menth23, menth26
## 4 1 1 marital logreg menth11, menth19
## 5 1 2 hlthinsur logreg menth11, menth19
## 6 1 2 medcost logreg menth11, menth19
#observe the new, imputed data
#we do not see any NA's
dat.imp<-complete(imp,action=1)
head(dat.imp,n=20)
## menth hlthinsur medcost sxorient marital
## 1 2 Y N S M
## 2 15 N N S S
## 3 0 Y N S M
## 4 0 Y N S M
## 5 0 Y N S M
## 6 0 Y Y S M
## 7 0 Y N S M
## 8 0 Y N S M
## 9 0 Y N S M
## 10 0 Y N S M
## 11 0 Y N S M
## 12 0 Y N S S
## 13 0 Y N S M
## 14 0 Y N S M
## 15 3 Y N S S
## 16 15 Y N S M
## 17 0 Y N S M
## 18 0 Y N S M
## 19 30 Y N S M
## 20 20 Y N S M
#compare to the original data
#notice the NA's
head(brfss2017[,c("menth","hlthinsur","medcost","sxorient","marital")],n=20)
## menth hlthinsur medcost sxorient marital
## 298797 2 Y N S M
## 246803 15 N N <NA> <NA>
## 360948 0 Y N S M
## 433713 0 Y N <NA> M
## 310183 0 Y N <NA> <NA>
## 315088 0 Y Y <NA> <NA>
## 256052 0 Y N <NA> M
## 317562 0 Y N <NA> <NA>
## 178508 0 Y N <NA> M
## 117801 0 Y N <NA> <NA>
## 64447 0 Y N S M
## 7229 0 Y N <NA> <NA>
## 419294 0 Y N <NA> M
## 292236 0 Y N S M
## 411807 3 Y N S <NA>
## 375166 15 Y N <NA> <NA>
## 211469 0 Y N <NA> <NA>
## 281572 0 Y N <NA> M
## 445953 30 Y N <NA> M
## 120685 20 Y N S M
#Let's look at the variablity in the 5 different imputations
#the 5 different imputations follow a similar pattern
fit.menth<-with(data=imp,expr=lm(menth~sxorient+marital+hlthinsur+medcost))
## Warning in model.response(mf, "numeric"): using type = "numeric" with a
## factor response will be ignored
## Warning in Ops.factor(y, z$residuals): '-' not meaningful for factors
## Warning in model.response(mf, "numeric"): using type = "numeric" with a
## factor response will be ignored
## Warning in Ops.factor(y, z$residuals): '-' not meaningful for factors
## Warning in model.response(mf, "numeric"): using type = "numeric" with a
## factor response will be ignored
## Warning in Ops.factor(y, z$residuals): '-' not meaningful for factors
## Warning in model.response(mf, "numeric"): using type = "numeric" with a
## factor response will be ignored
## Warning in Ops.factor(y, z$residuals): '-' not meaningful for factors
## Warning in model.response(mf, "numeric"): using type = "numeric" with a
## factor response will be ignored
## Warning in Ops.factor(y, z$residuals): '-' not meaningful for factors
fit.menth
## call :
## with.mids(data = imp, expr = lm(menth ~ sxorient + marital +
## hlthinsur + medcost))
##
## call1 :
## mice(data = brfss2017[, c("menth", "hlthinsur", "medcost", "sxorient",
## "marital")], m = 5, seed = 22)
##
## nmis :
## menth hlthinsur medcost sxorient marital
## 376 112 75 14232 7956
##
## analyses :
## [[1]]
##
## Call:
## lm(formula = menth ~ sxorient + marital + hlthinsur + medcost)
##
## Coefficients:
## (Intercept) sxorientS maritalS hlthinsurY medcostY
## 5.3161 -1.9596 1.9078 0.2718 4.5950
##
##
## [[2]]
##
## Call:
## lm(formula = menth ~ sxorient + marital + hlthinsur + medcost)
##
## Coefficients:
## (Intercept) sxorientS maritalS hlthinsurY medcostY
## 5.1878 -1.9503 2.0282 0.3634 4.5608
##
##
## [[3]]
##
## Call:
## lm(formula = menth ~ sxorient + marital + hlthinsur + medcost)
##
## Coefficients:
## (Intercept) sxorientS maritalS hlthinsurY medcostY
## 5.2895 -1.9688 1.8999 0.3093 4.5382
##
##
## [[4]]
##
## Call:
## lm(formula = menth ~ sxorient + marital + hlthinsur + medcost)
##
## Coefficients:
## (Intercept) sxorientS maritalS hlthinsurY medcostY
## 4.8801 -1.5894 1.9064 0.3398 4.6158
##
##
## [[5]]
##
## Call:
## lm(formula = menth ~ sxorient + marital + hlthinsur + medcost)
##
## Coefficients:
## (Intercept) sxorientS maritalS hlthinsurY medcostY
## 4.8865 -1.6205 1.8623 0.3665 4.6522
#variation in mental health
#variation in mental health is consistant across 5 imputations
with(data=imp,exp=sd(menth))
## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm = na.rm): Calling var(x) on a factor x is deprecated and will become an error.
## Use something like 'all(duplicated(x)[-1L])' to test for a constant vector.
## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm = na.rm): Calling var(x) on a factor x is deprecated and will become an error.
## Use something like 'all(duplicated(x)[-1L])' to test for a constant vector.
## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm = na.rm): Calling var(x) on a factor x is deprecated and will become an error.
## Use something like 'all(duplicated(x)[-1L])' to test for a constant vector.
## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm = na.rm): Calling var(x) on a factor x is deprecated and will become an error.
## Use something like 'all(duplicated(x)[-1L])' to test for a constant vector.
## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm = na.rm): Calling var(x) on a factor x is deprecated and will become an error.
## Use something like 'all(duplicated(x)[-1L])' to test for a constant vector.
## call :
## with.mids(data = imp, expr = sd(menth))
##
## call1 :
## mice(data = brfss2017[, c("menth", "hlthinsur", "medcost", "sxorient",
## "marital")], m = 5, seed = 22)
##
## nmis :
## menth hlthinsur medcost sxorient marital
## 376 112 75 14232 7956
##
## analyses :
## [[1]]
## [1] 7.901641
##
## [[2]]
## [1] 7.893235
##
## [[3]]
## [1] 7.888821
##
## [[4]]
## [1] 7.892257
##
## [[5]]
## [1] 7.88554
#frequency data for sxorient
#Frequency is consistant across 5 imputations
with(data=imp,exp=(prop.table(table(sxorient))))
## call :
## with.mids(data = imp, expr = (prop.table(table(sxorient))))
##
## call1 :
## mice(data = brfss2017[, c("menth", "hlthinsur", "medcost", "sxorient",
## "marital")], m = 5, seed = 22)
##
## nmis :
## menth hlthinsur medcost sxorient marital
## 376 112 75 14232 7956
##
## analyses :
## [[1]]
## sxorient
## G S
## 0.01612 0.98388
##
## [[2]]
## sxorient
## G S
## 0.01732 0.98268
##
## [[3]]
## sxorient
## G S
## 0.0174 0.9826
##
## [[4]]
## sxorient
## G S
## 0.01768 0.98232
##
## [[5]]
## sxorient
## G S
## 0.01764 0.98236
#frequency data for marital
#Frequency is consistant across 5 imputations
with(data=imp,exp=(prop.table(table(marital))))
## call :
## with.mids(data = imp, expr = (prop.table(table(marital))))
##
## call1 :
## mice(data = brfss2017[, c("menth", "hlthinsur", "medcost", "sxorient",
## "marital")], m = 5, seed = 22)
##
## nmis :
## menth hlthinsur medcost sxorient marital
## 376 112 75 14232 7956
##
## analyses :
## [[1]]
## marital
## M S
## 0.75932 0.24068
##
## [[2]]
## marital
## M S
## 0.75756 0.24244
##
## [[3]]
## marital
## M S
## 0.75916 0.24084
##
## [[4]]
## marital
## M S
## 0.75912 0.24088
##
## [[5]]
## marital
## M S
## 0.7588 0.2412
#frequency data for health insurance
#Frequency is consistant across 5 imputations
with(data=imp,exp=(prop.table(table(hlthinsur))))
## call :
## with.mids(data = imp, expr = (prop.table(table(hlthinsur))))
##
## call1 :
## mice(data = brfss2017[, c("menth", "hlthinsur", "medcost", "sxorient",
## "marital")], m = 5, seed = 22)
##
## nmis :
## menth hlthinsur medcost sxorient marital
## 376 112 75 14232 7956
##
## analyses :
## [[1]]
## hlthinsur
## N Y
## 0.07936 0.92064
##
## [[2]]
## hlthinsur
## N Y
## 0.07948 0.92052
##
## [[3]]
## hlthinsur
## N Y
## 0.07928 0.92072
##
## [[4]]
## hlthinsur
## N Y
## 0.07932 0.92068
##
## [[5]]
## hlthinsur
## N Y
## 0.07952 0.92048
#frequency data for medcost
#Frequency is consistant across 5 imputations
with(data=imp,exp=(prop.table(table(medcost))))
## call :
## with.mids(data = imp, expr = (prop.table(table(medcost))))
##
## call1 :
## mice(data = brfss2017[, c("menth", "hlthinsur", "medcost", "sxorient",
## "marital")], m = 5, seed = 22)
##
## nmis :
## menth hlthinsur medcost sxorient marital
## 376 112 75 14232 7956
##
## analyses :
## [[1]]
## medcost
## N Y
## 0.89332 0.10668
##
## [[2]]
## medcost
## N Y
## 0.89328 0.10672
##
## [[3]]
## medcost
## N Y
## 0.89324 0.10676
##
## [[4]]
## medcost
## N Y
## 0.8934 0.1066
##
## [[5]]
## medcost
## N Y
## 0.89328 0.10672
esp.p<-pool(fit.menth)
## Warning in Ops.factor(r, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(weighted.residuals(object), 2): '^' not meaningful
## for factors
## Warning in Ops.factor(r, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(weighted.residuals(object), 2): '^' not meaningful
## for factors
## Warning in Ops.factor(r, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(weighted.residuals(object), 2): '^' not meaningful
## for factors
## Warning in Ops.factor(r, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(weighted.residuals(object), 2): '^' not meaningful
## for factors
## Warning in Ops.factor(r, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(res, 2): '^' not meaningful for factors
## Warning in Ops.factor(weighted.residuals(object), 2): '^' not meaningful
## for factors
## Warning in Ops.factor(r, 2): '^' not meaningful for factors
## Warning in Ops.factor(r, 2): '^' not meaningful for factors
## Warning in Ops.factor(r, 2): '^' not meaningful for factors
## Warning in Ops.factor(r, 2): '^' not meaningful for factors
## Warning in Ops.factor(r, 2): '^' not meaningful for factors
print(esp.p)
## Class: mipo m = 5
## estimate ubar b t dfcom df riv lambda fmi
## (Intercept) 5.1119784 NA 0.045883302 NA 24995 NA NA NA NA
## sxorientS -1.8177208 NA 0.037892597 NA 24995 NA NA NA NA
## maritalS 1.9209147 NA 0.003942219 NA 24995 NA NA NA NA
## hlthinsurY 0.3301485 NA 0.001589063 NA 24995 NA NA NA NA
## medcostY 4.5924095 NA 0.002015001 NA 24995 NA NA NA NA
#We see significance across all variables except for sexual orientation. When an intersection is created with sexual orientation and marital status we do see significance for married same sex