Plot Map of Veterans (Sample Weighted)
plot_usmap(data=myprop, values="Veteran",
labels = TRUE, label_color = "white",
color="red")+
scale_fill_continuous(type = "viridis"
)+
theme(legend.position="right", plot.title =
element_text(hjust = 0.5, face="bold"))+
labs(fill="Proportion of Veterans")
## Warning: Use of `map_df$x` is discouraged. Use `x` instead.
## Warning: Use of `map_df$y` is discouraged. Use `y` instead.
## Warning: Use of `map_df$group` is discouraged. Use `group` instead.
## Warning: Use of `centroid_labels$x` is discouraged. Use `x` instead.
## Warning: Use of `centroid_labels$y` is discouraged. Use `y` instead.
## Warning: Use of `centroid_labels$abbr` is discouraged. Use `abbr` instead.
Apply Survey Weights and Design
require(survey)
## Loading required package: survey
## Warning: package 'survey' was built under R version 3.5.3
## Loading required package: grid
## Loading required package: Matrix
## Warning: package 'Matrix' was built under R version 3.5.3
##
## Attaching package: 'survey'
## The following object is masked from 'package:graphics':
##
## dotchart
# Set options for allowing a single observation per stratum
options(survey.lonely.psu = "adjust")
mysurvey <- svydesign(
id=~1,
strata = ~Stratum,
weights = ~Weights,
data = mydata)
Make Factors from Numeric Variables
library(psych)
## Warning: package 'psych' was built under R version 3.5.3
##
## Attaching package: 'psych'
## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha
newvars=mysurvey$variables
newvars$Age=as.factor(newvars$Age)
newvars$Race=as.factor(newvars$Race)
newvars$Married=as.factor(newvars$Married)
newvars$Income=as.factor(newvars$Income)
newvars$Education=as.factor(newvars$Education)
newvars$Employed=as.factor(newvars$Employed)
newvars$State=as.factor(newvars$state)
Descriptives
Boxplots
myplots=svyby(~BMI+Arthritis+Depression+Diabetes+SkinCancer+
Cancer+CHD+COPD+
Kidney+Stroke
,~Vet,mysurvey, svymean)
op <- par(cex = 0.7)
barplot(myplots, beside=TRUE, legend=TRUE, col=terrain.colors(10), ylab="Proportion", xlab="Veteran Status", ylim=c(0,.8))
Proportions by Vet Status
myplots=svyby(~BMI+Gender+Arthritis+Depression+Diabetes+SkinCancer+
Cancer+CHD+COPD+
Kidney+Stroke
,~Vet,mysurvey, svymean)
t(round(myplots,3))
## 0 1
## Vet 0.000 1.000
## BMI 0.688 0.758
## Gender 0.416 0.893
## Arthritis 0.287 0.380
## Depression 0.187 0.163
## Diabetes 0.122 0.183
## SkinCancer 0.073 0.147
## Cancer 0.078 0.128
## CHD 0.044 0.108
## COPD 0.073 0.107
## Kidney 0.034 0.049
## Stroke 0.036 0.065
## se.BMI 0.001 0.003
## se.Gender 0.001 0.002
## se.Arthritis 0.001 0.004
## se.Depression 0.001 0.003
## se.Diabetes 0.001 0.003
## se.SkinCancer 0.001 0.003
## se.Cancer 0.001 0.003
## se.CHD 0.001 0.002
## se.COPD 0.001 0.002
## se.Kidney 0.001 0.002
## se.Stroke 0.001 0.002
Categorical Distributions by Vet Status
round(svytable(~Age+Vet, mysurvey))
## Vet
## Age 0 1
## 1 18448002 1002954
## 2 31737049 2585366
## 3 33562529 2606825
## 4 39393961 4176664
## 5 45766172 5269496
## 6 57936867 15587087
round(svytable(~Race+Vet, mysurvey))
## Vet
## Race 0 1
## 1 150923402 23924151
## 2 23141574 3082954
## 3 8827609 450042
## 4 2835217 477955
## 5 34058788 2067875
## 6 7057989 1225416
round(svytable(~Married+Vet, mysurvey))
## Vet
## Married 0 1
## 1 117244545 19004330
## 2 27700300 4698798
## 3 19125947 2870047
## 4 5705795 625682
## 5 47306108 3365352
## 6 9761886 664183
round(svytable(~Income+Vet, mysurvey))
## Vet
## Income 0 1
## 1 9981703 590902
## 2 9446934 915126
## 3 13139276 1299729
## 4 15995550 1975529
## 5 17823584 2884678
## 6 22715736 3979065
## 7 27647254 5005299
## 8 69554937 10167270
## 9 40539605 4410794
round(svytable(~Education+Vet, mysurvey))
## Vet
## Education 0 1
## 1 565564 21420
## 2 7676759 243225
## 3 12502001 763200
## 4 56331664 7529521
## 5 58460447 10131974
## 6 91308144 12539052
round(svytable(~Employed+Vet, mysurvey))
## Vet
## Employed 0 1
## 1 107555934 11407575
## 2 21608269 2199635
## 3 5063241 480753
## 4 5309590 409354
## 5 13090478 225333
## 6 8944506 459432
## 7 48444756 14105300
## 8 16827806 1941011
round(svytable(~state+Vet,mysurvey))
## Vet
## state 0 1
## Alabama 3233678 576580
## Alaska 462167 94222
## Arizona 4617814 879501
## Arkansas 2004811 316407
## California 27912327 2893984
## Colorado 3886852 558015
## Connecticut 2560492 292154
## Delaware 664196 104928
## District of Columbia 528475 48932
## Florida 14589586 2616746
## Georgia 6980609 1093479
## Guam 89141 17905
## Hawaii 954759 167416
## Idaho 1119513 190960
## Illinois 9039076 871581
## Indiana 4547872 594000
## Iowa 2154207 280038
## Kansas 1924749 287840
## Kentucky 3107055 364340
## Louisiana 3173899 421327
## Maine 933395 155907
## Maryland 4137658 613030
## Massachusetts 4985265 554590
## Michigan 7038835 791335
## Minnesota 3881855 453944
## Mississippi 1995789 283290
## Missouri 4122989 639648
## Montana 705455 128720
## Nebraska 1275353 185260
## Nevada 2058387 317689
## New Hampshire 935457 159755
## New Jersey 6471842 606634
## New Mexico 1393960 211654
## New York 14444052 1283177
## North Carolina 6969404 1153448
## North Dakota 500037 80584
## Ohio 8055130 1055803
## Oklahoma 2582159 414473
## Oregon 2888845 447519
## Pennsylvania 8991610 1186068
## Puerto Rico 2613569 124048
## Rhode Island 746239 109610
## South Carolina 3416893 583880
## South Dakota 567254 96439
## Tennessee 4558582 730816
## Texas 18157415 3318213
## Utah 2024980 213214
## Vermont 449595 57449
## Virginia 5609508 1065707
## Washington 5067604 824309
## West Virginia 1264811 173838
## Wisconsin 4069188 476301
## Wyoming 380190 61683
Models
BMI
model1<-svyglm(na.omit(BMI~as.factor(Age)+as.factor(Race)+as.factor(Gender)+
as.factor(Married)+as.factor(Income)+
as.factor(Education)+as.factor(Employed)+as.factor(state)+
as.factor(Vet)),
design=mysurvey,family=quasibinomial)
model1$oddsratios=exp(model1$coefficients)
model1$lower=exp(log(model1$oddsratios)-qnorm(.975)*summary(model1)$coefficients[,2])
model1$upper=exp(log(model1$oddsratios)+qnorm(.975)*summary(model1)$coefficients[,2])
options(digits=3)
model1$lower[90]
## as.factor(Vet)1
## 1.07
model1$oddsratios[90]
## as.factor(Vet)1
## 1.12
model1$upper[90]
## as.factor(Vet)1
## 1.17
#summary(model1)
Arthritis
#models 2 thr
model2<-svyglm(na.omit(Arthritis~as.factor(Age)+as.factor(Race)+as.factor(Gender)+
as.factor(Married)+as.factor(Income)+
as.factor(Education)+as.factor(Employed)+
as.factor(state)+ as.factor(Vet)),
design=mysurvey,family=quasibinomial)
model2$oddsratios=exp(model2$coefficients)
model2$lower=exp(log(model2$oddsratios)-qnorm(.975)*summary(model1)$coefficients[,2])
model2$upper=exp(log(model2$oddsratios)+qnorm(.975)*summary(model1)$coefficients[,2])
options(digits=3)
model2$lower[90]
## as.factor(Vet)1
## 1.21
model2$oddsratios[90]
## as.factor(Vet)1
## 1.26
model2$upper[90]
## as.factor(Vet)1
## 1.31
#summary(model2)
Depression
model3<-svyglm(na.omit(Depression~as.factor(Age)+as.factor(Race)+as.factor(Gender)+
as.factor(Married)+as.factor(Income)+
as.factor(Education)+as.factor(Employed)+
as.factor(state)+ as.factor(Vet)),
design=mysurvey,family=quasibinomial)
model3$oddsratios=exp(model3$coefficients)
model3$lower=exp(log(model3$oddsratios)-qnorm(.975)*summary(model1)$coefficients[,2])
model3$upper=exp(log(model3$oddsratios)+qnorm(.975)*summary(model1)$coefficients[,2])
options(digits=3)
model3$lower[90]
## as.factor(Vet)1
## 1.24
model3$oddsratios[90]
## as.factor(Vet)1
## 1.29
model3$upper[90]
## as.factor(Vet)1
## 1.35
Diabetes
model4<-svyglm(na.omit(Diabetes~as.factor(Age)+as.factor(Race)+as.factor(Gender)+
as.factor(Married)+as.factor(Income)+
as.factor(Education)+as.factor(Employed)+
as.factor(state)+ as.factor(Vet)),
design=mysurvey,family=quasibinomial)
model4$oddsratios=exp(model4$coefficients)
model4$lower=exp(log(model4$oddsratios)-qnorm(.975)*summary(model4)$coefficients[,2])
model4$upper=exp(log(model4$oddsratios)+qnorm(.975)*summary(model4)$coefficients[,2])
options(digits=3)
model4$lower[90]
## as.factor(Vet)1
## 1.05
model4$oddsratios[90]
## as.factor(Vet)1
## 1.1
model4$upper[90]
## as.factor(Vet)1
## 1.16
Skin Cancer
model5<-svyglm(na.omit(SkinCancer~as.factor(Age)+as.factor(Race)+as.factor(Gender)+
as.factor(Married)+as.factor(Income)+
as.factor(Education)+as.factor(Employed)+
as.factor(state)+ as.factor(Vet)),
design=mysurvey,family=quasibinomial)
model5$oddsratios=exp(model5$coefficients)
model5$lower=exp(log(model5$oddsratios)-qnorm(.975)*summary(model5)$coefficients[,2])
model5$upper=exp(log(model5$oddsratios)+qnorm(.975)*summary(model5)$coefficients[,2])
options(digits=3)
model5$lower[90]
## as.factor(Vet)1
## 1.12
model5$oddsratios[90]
## as.factor(Vet)1
## 1.18
model5$upper[90]
## as.factor(Vet)1
## 1.25
Cancer
model6<-svyglm(na.omit(Cancer~as.factor(Age)+as.factor(Race)+as.factor(Gender)+
as.factor(Married)+as.factor(Income)+
as.factor(Education)+as.factor(Employed)+
as.factor(state)+ as.factor(Vet)),
design=mysurvey,family=quasibinomial)
model6$oddsratios=exp(model6$coefficients)
model6$lower=exp(log(model6$oddsratios)-qnorm(.975)*summary(model6)$coefficients[,2])
model6$upper=exp(log(model6$oddsratios)+qnorm(.975)*summary(model6)$coefficients[,2])
options(digits=3)
model6$lower[90]
## as.factor(Vet)1
## 1.3
model6$oddsratios[90]
## as.factor(Vet)1
## 1.38
model6$upper[90]
## as.factor(Vet)1
## 1.47
CHD
model7<-svyglm(na.omit(CHD~as.factor(Age)+as.factor(Race)+as.factor(Gender)+
as.factor(Married)+as.factor(Income)+
as.factor(Education)+as.factor(Employed)+
as.factor(state)+ as.factor(Vet)),
design=mysurvey,family=quasibinomial)
model7$oddsratios=exp(model7$coefficients)
model7$lower=exp(log(model7$oddsratios)-qnorm(.975)*summary(model7)$coefficients[,2])
model7$upper=exp(log(model7$oddsratios)+qnorm(.975)*summary(model7)$coefficients[,2])
options(digits=3)
model7$lower[90]
## as.factor(Vet)1
## 1.24
model7$oddsratios[90]
## as.factor(Vet)1
## 1.32
model7$upper[90]
## as.factor(Vet)1
## 1.41
COPD
model8<-svyglm(na.omit(COPD~as.factor(Age)+as.factor(Race)+as.factor(Gender)+
as.factor(Married)+as.factor(Income)+
as.factor(Education)+as.factor(Employed)+
as.factor(state)+ as.factor(Vet)),
design=mysurvey,family=quasibinomial)
model8$oddsratios=exp(model8$coefficients)
model8$lower=exp(log(model8$oddsratios)-qnorm(.975)*summary(model8)$coefficients[,2])
model8$upper=exp(log(model8$oddsratios)+qnorm(.975)*summary(model8)$coefficients[,2])
options(digits=3)
model8$lower[90]
## as.factor(Vet)1
## 1.37
model8$oddsratios[90]
## as.factor(Vet)1
## 1.46
model8$upper[90]
## as.factor(Vet)1
## 1.56
Kidney
model9<-svyglm(na.omit(Kidney~as.factor(Age)+as.factor(Race)+as.factor(Gender)+
as.factor(Married)+as.factor(Income)+
as.factor(Education)+as.factor(Employed)+
as.factor(state)+ as.factor(Vet)),
design=mysurvey,family=quasibinomial)
model9$oddsratios=exp(model9$coefficients)
model9$lower=exp(log(model9$oddsratios)-qnorm(.975)*summary(model9)$coefficients[,2])
model9$upper=exp(log(model9$oddsratios)+qnorm(.975)*summary(model9)$coefficients[,2])
options(digits=3)
model9$lower[90]
## as.factor(Vet)1
## 0.986
model9$oddsratios[90]
## as.factor(Vet)1
## 1.08
model9$upper[90]
## as.factor(Vet)1
## 1.19
Stroke
model10<-svyglm(na.omit(Stroke~as.factor(Age)+as.factor(Race)+as.factor(Gender)+
as.factor(Married)+as.factor(Income)+
as.factor(Education)+as.factor(Employed)+
as.factor(state)+ as.factor(Vet)),
design=mysurvey,family=quasibinomial)
model10$oddsratios=exp(model10$coefficients)
model10$lower=exp(log(model10$oddsratios)-qnorm(.975)*summary(model10)$coefficients[,2])
model10$upper=exp(log(model10$oddsratios)+qnorm(.975)*summary(model10)$coefficients[,2])
options(digits=3)
model10$lower[90]
## as.factor(Vet)1
## 1.17
model10$oddsratios[90]
## as.factor(Vet)1
## 1.28
model10$upper[90]
## as.factor(Vet)1
## 1.39