fsolt / dotwhisker

Dot-and-Whisker Plots of Regression Results
https://fsolt.org/dotwhisker/
Other
60 stars 10 forks source link

Dodge and order vars not working dwplot #72

Closed aaanders05 closed 6 years ago

aaanders05 commented 7 years ago

ABDIII_EXLAPA.xlsx I have the attached data and am trying to create a dot and whisker plot using the following lines of code:

Creating the initial dataset

table3 <- read.csv("ABDIII_EXLAPA.csv")

Subsetting the data to include the columns necessary for analysis

newdata <- table3[c(1:4,6,8,9,13:31,50:53,61,63:66,70:75,77:79)]

Removing non-responses

mydata <- na.omit(newdata)

Model 1d: Protest with survey glm package and binary internet use, org.membership measures, motivations, grievance measures and all other variables

mydata$INCOME_QUINT_2 <- as.factor(mydata$INCOME_QUINT_1) mydata$EDUC_BAS_2 <- as.factor(mydata$EDUC_BAS)

svydata1d <- svydesign(id= ~PSU, strata= ~STRATA_2, weights=~WT, data=mydata)

summary(svydata1d)

options(survey.lonely.psu="adjust")

mysvylogit1d <- svyglm(PROTEST ~ INTERNET_1_BIN + FORM_ORG_BIN + POLIT_KNOW + POLIT_INT + ABS_GRIEVE + REL_GRIEVE + INCOME_QUINT_2 + QURAN + GENDER + SINGLE + EDUC_BAS_2 + STUDENT + EMPLOYED + AGE, design = svydata1d, family=quasibinomial(link=logit))

summary(mysvylogit1d)

Model 1f: Protest with survey glm package with all variables and organizational type

mydata$ORG_TYPE_21 <- as.factor(mydata$ORG_TYPE_2)

svydata1f <- svydesign(id= ~PSU, strata= ~STRATA_2, weights=~WT, data=mydata)

summary(svydata1f)

options(survey.lonely.psu="adjust")

mysvylogit1f <- svyglm(PROTEST ~ INTERNET_1_BIN + ORG_TYPE_21+ POLIT_KNOW + POLIT_INT + ABS_GRIEVE + REL_GRIEVE + INCOME_QUINT_2 + QURAN + GENDER + SINGLE + EDUC_BAS_2 + STUDENT + EMPLOYED + AGE, design = svydata1f, family=quasibinomial(link=logit))

summary(mysvylogit1f)

coefplot1 <- dwplot(list(mysvylogit1d, mysvylogit1f), dodge_size = 0.8, order_vars = c("INTERNET_1_BIN", "FORM_ORG_BIN", "ORG_TYPE_211", "ORG_TYPE_212", "POLIT_KNOW", "POLIT_INT", "ABS_GRIEVE", "REL_GRIEVE", "INCOME_QUINT_22", "INCOME_QUINT_23", "INCOME_QUINT_24", "INCOME_QUINT_25", "QURAN", "GENDER", "SINGLE", "EDUC_BAS_22", "EDUC_BAS_23", "STUDENT", "EMPLOYED", "AGE")) %>% relabel_predictors(c(INTERNET_1_BIN = "Internet use", FORM_ORG_BIN = "Organizational membership", POLIT_KNOW = "Political knowledge", POLIT_INT = "Political interest", ABS_GRIEVE = "Absolute grievance", REL_GRIEVE = "Relative grievance", INCOME_QUINT_22 = "Income quintile (21-40)", INCOME_QUINT_23 = "Income quintile (41-60)", INCOME_QUINT_24 = "Income quintile (61-80)", INCOME_QUINT_25 = "Income quintile (81-100)", QURAN = "Religiosity", GENDER = "Gender", SINGLE = "Marital status", EDUC_BAS_22 = "High School Educated", EDUC_BAS_23 = "College Educated", STUDENT = "Student", EMPLOYED = "Employment status", AGE = "Age", ORG_TYPE_211 = "Expressive org. membership", ORG_TYPE_212 = "Instrumental org. membership")) coefplot1 <- coefplot1 + theme_bw() + theme(panel.grid.minor.y = element_blank(), panel.grid.major.y=element_blank(), panel.grid.minor.x = element_blank(), panel.grid.major.x=element_blank()) + xlab("Coefficient Estimate") + ylab("") + xlim(-1,2) + geom_vline(xintercept = 0, colour="grey60", linetype=2) + theme(legend.position = "none") + annotate("text", x = 1.5, y = 10, size = 3, hjust = 0, label = "R^2 == 0.169", parse = TRUE) + annotate("text", x = 1.5, y = 9, size = 3, hjust = 0, label = "N = 6663") + theme(axis.text.y = element_text(family="CMUSansSerif", size = 11), axis.text.x = element_text(family="CMUSansSerif", size = 12), axis.title=element_text(family = "CMUSansSerif", size=14)) + ggtitle("Determinants of individual protest participation") + theme(plot.title = element_text(family = "CMUSansSerif", face = "bold", size = 14, hjust = 0.5))

However, the final plot does not display all estimates for model 2. Is there a way to fix this? Also even if I remove the order_vars specification, the dodge command doesn't create the spacing I need. Is there a way to rectify this issue as well?

stefan-mueller commented 7 years ago

Some general remarks first.

No offense taken, but if you require help, make sure that we don't have to wrangle your code for 15 minutes.

Why didn't your code work? Basically, because ordering and relabeling needs to be applied to a data frame. Please see my solution below (and some comments).

library(survey)
library(dotwhisker)
library(dplyr)
library(tidyverse)

##Creating the initial dataset (please upload correct data!)
table3 <- read.csv("ABDIII_EXLAPA.csv")

##Subsetting the data to include the columns necessary for analysis
# NOTE: not a good way to subset data -- look at select and filter from dplyr

newdata <- table3[c(1:4,6,8,9,13:31,50:53,61,63:66,70:75,77:79)]

##Removing non-responses
mydata <- na.omit(newdata)

##Model 1d: Protest with survey glm package and binary internet use, org.membership measures, motivations, grievance measures and all other variables

mydata$INCOME_QUINT_2 <- as.factor(mydata$INCOME_QUINT_1)
mydata$EDUC_BAS_2 <- as.factor(mydata$EDUC_BAS)

svydata1d <- svydesign(id = ~PSU, strata = ~STRATA_2, 
                       weights = ~WT, data = mydata)

summary(svydata1d)

options(survey.lonely.psu="adjust")

## Run models

mysvylogit1d <- svyglm(PROTEST ~ INTERNET_1_BIN + 
                         FORM_ORG_BIN + # not in model 1f
                         POLIT_KNOW + 
                         POLIT_INT + 
                         ABS_GRIEVE + 
                         REL_GRIEVE + 
                         INCOME_QUINT_2 + 
                         QURAN + 
                         GENDER + 
                         SINGLE + 
                         EDUC_BAS_2 + 
                         STUDENT + 
                         EMPLOYED + 
                         AGE,
                         design = svydata1d, 
                       family=quasibinomial(link=logit))

summary(mysvylogit1d)

mydata$ORG_TYPE_21 <- as.factor(mydata$ORG_TYPE_2)

svydata1f <- svydesign(id= ~PSU, strata= ~STRATA_2, weights=~WT, data=mydata)

options(survey.lonely.psu="adjust")

mysvylogit1f <- svyglm(PROTEST ~ INTERNET_1_BIN + 
                         POLIT_KNOW + 
                         POLIT_INT + 
                         ABS_GRIEVE + 
                         REL_GRIEVE + 
                         INCOME_QUINT_2 + 
                         QURAN + 
                         GENDER + 
                         SINGLE + 
                         EDUC_BAS_2 + 
                         STUDENT + 
                         EMPLOYED + 
                         AGE +
                         ORG_TYPE_21, # not in model 1e
                          design = svydata1f, 
                       family = quasibinomial(link=logit))

summary(mysvylogit1f)

# Disply a basic dwplot
dwplot(list(mysvylogit1d, mysvylogit1f))

# works correctly

order_vars_original <- c("INTERNET_1_BIN", 
                         "FORM_ORG_BIN", 
                         "ORG_TYPE_211", 
                         "ORG_TYPE_212",  
                         "POLIT_KNOW", 
                         "POLIT_INT", 
                         "ABS_GRIEVE", 
                         "REL_GRIEVE", 
                         "INCOME_QUINT_22",
                         "INCOME_QUINT_23", 
                         "INCOME_QUINT_24", 
                         "INCOME_QUINT_25", 
                         "QURAN", 
                         "GENDER", 
                         "SINGLE", 
                         "EDUC_BAS_22", 
                         "EDUC_BAS_23", 
                         "STUDENT", 
                         "EMPLOYED", 
                         "AGE")

data_plot <- rbind(broom::tidy(mysvylogit1d) %>%  mutate(model = "Model 1"),
              broom::tidy(mysvylogit1f) %>%  mutate(model = "Model 2")) %>% 
  filter(term != "(Intercept)") %>%  
  mutate(term = factor(term, levels = order_vars_original)) %>% 
  relabel_predictors(c(INTERNET_1_BIN = "Internet use", 
                       FORM_ORG_BIN = "Organizational membership", 
                       POLIT_KNOW = "Political knowledge", 
                       POLIT_INT = "Political interest", 
                       ABS_GRIEVE = "Absolute grievance", 
                       REL_GRIEVE = "Relative grievance", 
                       INCOME_QUINT_22 = "Income quintile (21-40)", 
                       INCOME_QUINT_23 = "Income quintile (41-60)", 
                       INCOME_QUINT_24 = "Income quintile (61-80)", 
                       INCOME_QUINT_25 = "Income quintile (81-100)", 
                       QURAN = "Religiosity", 
                       GENDER = "Gender", 
                       SINGLE = "Marital status", 
                       EDUC_BAS_22 = "High School Educated", 
                       EDUC_BAS_23 = "College Educated", 
                       STUDENT = "Student", 
                       EMPLOYED = "Employment status", 
                       AGE = "Age", 
                       ORG_TYPE_211 = "Expressive org. membership", 
                       ORG_TYPE_212 = "Instrumental org. membership"))

## Only use each option once (e.g. theme) and improve structure
dwplot(data_plot, dodge_size = 0.8) +
  theme_bw() + 
  ggtitle("Determinants of individual protest participation") + 
  xlab("Coefficient Estimate") + 
  ylab("") + 
  xlim(-1,2) + 
  geom_vline(xintercept = 0, colour="grey60", linetype=2) + 
  theme(legend.position = "none") +
  annotate("text", x = 1.5, y = 10, size = 3, hjust = 0, label = "R^2 == 0.169", parse = TRUE) + annotate("text", x = 1.5, y = 9, size = 3, hjust = 0, label = "N = 6663") + 
  theme(axis.text.y = element_text(size = 11), 
        axis.text.x = element_text(size = 12), 
        axis.title=element_text(size=14),
        plot.title = element_text(face = "bold", size = 14, hjust = 0.5),
        legend.position = "none",
        panel.grid.minor.y = element_blank(), 
        panel.grid.major.y=element_blank(), 
        panel.grid.minor.x = element_blank(), 
        panel.grid.major.x=element_blank()) 

Edit: I see that there is still something wrong with the ordering/renaming, but this code seems to work. So go from there and check the parts from order_vars_original onwards. If there are serious issues you cannot solve, please let me know.

dwplot(list(mysvylogit1d, mysvylogit1f))
fsolt commented 6 years ago

It looks like @stefan-mueller solved your issue, @aaanders05, so I'm going to close this. If not, please feel welcome to reopen and let us know what else needs fixed.