HRTC Stats
Hunter Ratliff, @HunterRatliff1
1/25/2019
If you’d like the source data, reach out to the BHAR office or DM me on Twitter
knitr::opts_chunk$set(warning=F, message=F, echo=T)
require(MASS)
require(tidyverse)
require(lubridate)
require(readr)
require(ggplot2)
require(knitr)
require(scales)
library(stargazer)
library(lme4)
require(survival)
library(survminer)
library(survMisc)
library(pander)
panderOptions("table.alignment.default", "left")
users <- readRDS("Data/users.RDS") %>%
mutate(SurvEnd2 = ifelse(SurvEnd > 135, 135, SurvEnd))
HRTC <- readRDS("Data/MergedDf2.RDS") %>%
mutate(Type = ifelse(campg.day %in% c(26), "Knowledge w/ EX", Type)) %>% # Recode this hotline as MORE
mutate(Type = ifelse(campg.day %in% c(63,82,86,124), "Knowledge w/ EX", Type)) %>% # Recode this fact as MORE
filter(campg.day < SurvEnd) %>%
filter(Type!="Open-ended") %>%
group_by(userID) %>%
mutate(userResponses = sum(Replied, na.rm = T)) %>%
ungroup() %>%
mutate(
Type = recode(Type, "Agree/disagree"="T/F"),
Type = recode(Type, "Knowledge w/ more info"="Knowledge w/ EX")
)
user.ids <- sample_n(count(HRTC, userID),10)$userID
HRTC[HRTC$mid=="m032",]$Category <- "Partner"
HRTC[HRTC$mid=="m038",]$Category <- "Partner"
HRTC[HRTC$mid=="m061",]$Category <- "Abuse/IPV"
HRTC <- HRTC %>%
mutate(Tag = Keywords) %>%
mutate(Tag = ifelse(stringr::str_detect(Tag, "Cyber"), "Cyber", Tag)) %>%
mutate(Tag = ifelse(stringr::str_detect(Tag, "Consent"), "Consent", Tag)) %>%
mutate(Tag = ifelse(stringr::str_detect(Tag, "Individuality"), "Individuality", Tag)) %>%
mutate(Tag = ifelse(stringr::str_detect(Tag, "Stalking"), "Stalking", Tag)) %>%
mutate(Tag = ifelse(stringr::str_detect(Tag, "Emotional"), "Emotional", Tag))df <- HRTC %>%
mutate(Replied = as.character(Replied)) %>%
mutate(
Replied = recode(Replied, "0"="No"),
Replied = recode(Replied, "1"="Yes")
) %>%
count(campg.day, Type, Replied) %>%
spread(Replied, n) %>%
group_by(Type) %>%
summarise(Yes = sum(Yes, na.rm=T), No=sum(No, na.rm=T)) %>%
mutate(
`Reply %` = Yes/(Yes+No),
n = Yes + No
)
conTable <- as.table(rbind(df$Yes, df$No))
dimnames(conTable) <- list(Replied = c("Yes", "No"),
qType = df$Type)Graph of response rate to each question type
Here’s just an overview of the crude percent of people who responded to each question type:
df %>% mutate(`Reply %` = percent(`Reply %`)) %>% kable()| Type | Yes | No | Reply % | n |
|---|---|---|---|---|
| Hotline | 19 | 1725 | 1.1% | 1744 |
| Knowledge | 78 | 4082 | 1.9% | 4160 |
| Knowledge w/ EX | 1318 | 6387 | 17.1% | 7705 |
| Quiz | 818 | 2290 | 26.3% | 3108 |
| T/F | 1482 | 3222 | 31.5% | 4704 |
df %>% ggplot(aes(x=Type, y=`Reply %`, fill=n)) +
geom_col() +
ggthemes::scale_fill_continuous_tableau(palette = "Blue") +
scale_y_continuous(labels=percent) + coord_flip()
TODO: Effect size..?
Contengency table
Contengency table of responded (yes, no) vs the type of question
kable(conTable)| Hotline | Knowledge | Knowledge w/ EX | Quiz | T/F | |
|---|---|---|---|---|---|
| Yes | 19 | 78 | 1318 | 818 | 1482 |
| No | 1725 | 4082 | 6387 | 2290 | 3222 |
Test Results
Here’s the test results
chisq.test(conTable)
##
## Pearson's Chi-squared test
##
## data: conTable
## X-squared = 1848.9, df = 4, p-value < 2.2e-16
chisq.test(conTable, simulate.p.value = TRUE, B = 10000)
##
## Pearson's Chi-squared test with simulated p-value (based on 10000
## replicates)
##
## data: conTable
## X-squared = 1848.9, df = NA, p-value = 9.999e-05fifer::chisq.plot(conTable)Residuals
chisq.test(conTable)$residuals %>%
round(2) %>% kable()| Hotline | Knowledge | Knowledge w/ EX | Quiz | T/F | |
|---|---|---|---|---|---|
| Yes | -16.30 | -23.96 | -0.50 | 12.02 | 23.32 |
| No | 7.47 | 10.97 | 0.23 | -5.50 | -10.68 |
Post-hoc analysis
Using Fisher’s Exact. Agrees nicely with the results from the logistic regression (next section)
fifer::chisq.post.hoc(t(conTable), control = "bonferroni", digits = 8) %>% pander()## Adjusted p-values used the bonferroni method.| comparison | raw.p | adj.p |
|---|---|---|
| Hotline vs. Knowledge | 0.03254 | 0.3254 |
| Hotline vs. Knowledge w/ EX | 0 | 0 |
| Hotline vs. Quiz | 0 | 0 |
| Hotline vs. T/F | 0 | 0 |
| Knowledge vs. Knowledge w/ EX | 0 | 0 |
| Knowledge vs. Quiz | 0 | 0 |
| Knowledge vs. T/F | 0 | 0 |
| Knowledge w/ EX vs. Quiz | 0 | 0 |
| Knowledge w/ EX vs. T/F | 0 | 0 |
| Quiz vs. T/F | 8.6e-07 | 8.57e-06 |
## Fisher's Exact
# fisher.test(conTable,simulate.p.value=TRUE,B=1e7)
fifer::chisq.post.hoc(t(conTable), control = "fdr", digits = 8) %>% pander()## Adjusted p-values used the fdr method.| comparison | raw.p | adj.p |
|---|---|---|
| Hotline vs. Knowledge | 0.03254 | 0.03254 |
| Hotline vs. Knowledge w/ EX | 0 | 0 |
| Hotline vs. Quiz | 0 | 0 |
| Hotline vs. T/F | 0 | 0 |
| Knowledge vs. Knowledge w/ EX | 0 | 0 |
| Knowledge vs. Quiz | 0 | 0 |
| Knowledge vs. T/F | 0 | 0 |
| Knowledge w/ EX vs. Quiz | 0 | 0 |
| Knowledge w/ EX vs. T/F | 0 | 0 |
| Quiz vs. T/F | 8.6e-07 | 9.5e-07 |
Note: These regressions are run on students who responded at least once
Some resources for binomal regression (for my reference):
Run a bimodal regression based on Question type + grade + race + gender + number of user's responses + Campaign number (e.g. is it the first message sent, second sent, last message sent in the campaign). In the future, can look at how category and/or keywords play into this.
The next table () shows the model in a more pretty table format. shows the same model, but with an output more similar to SPSS.
df <- HRTC %>%
filter(userResponses>0) %>%
# filter(Type!="Knowledge") %>%
mutate(Dummy = as.factor(zipcode)) %>%
mutate(userID = as.factor(userID)) %>%
mutate(Category = as.factor(Category)) %>%
mutate(Replied = ifelse(Replied==1, "Yes", "No")) %>%
mutate(Replied = factor(Replied, levels=c("No","Yes")))
mod <- glm(Replied ~ Type + grade + race + gender + userResponses + campg.num,
data=df, family=binomial(link = "logit"))
# round(exp(coef(mod)),4)pander(mod)| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | -6.354 | 0.3701 | -17.17 | 4.606e-66 |
| TypeKnowledge | 0.3154 | 0.3071 | 1.027 | 0.3044 |
| TypeKnowledge w/ EX | 3.931 | 0.2709 | 14.51 | 1.009e-47 |
| TypeQuiz | 4.912 | 0.2784 | 17.64 | 1.197e-69 |
| TypeT/F | 5.029 | 0.2737 | 18.38 | 2.05e-75 |
| grade11 | -0.09469 | 0.08896 | -1.064 | 0.2871 |
| grade12 | -0.0538 | 0.08536 | -0.6303 | 0.5285 |
| grade7 | 0.07815 | 0.1302 | 0.6004 | 0.5482 |
| grade8 | 0.03561 | 0.1804 | 0.1974 | 0.8435 |
| grade9 | -0.095 | 0.0932 | -1.019 | 0.3081 |
| raceAsian | 0.1379 | 0.2794 | 0.4937 | 0.6215 |
| raceBlack | 0.2564 | 0.2511 | 1.021 | 0.3072 |
| raceHispanic | 0.1618 | 0.2426 | 0.6667 | 0.505 |
| raceOther | 0.2364 | 0.2966 | 0.797 | 0.4254 |
| raceWhite | 0.1511 | 0.2419 | 0.6248 | 0.5321 |
| genderMale | -0.008795 | 0.06531 | -0.1347 | 0.8929 |
| genderOther | -0.1408 | 0.3921 | -0.3591 | 0.7195 |
| userResponses | 0.1489 | 0.003566 | 41.76 | 0 |
| campg.num | -0.03787 | 0.002095 | -18.08 | 4.631e-73 |
# summary(mod)
mod.table <- as.data.frame(summary(mod)$coefficients) %>%
mutate(`Exp(Estimate)` = exp(Estimate))
mod.exp.ci <- exp(confint(mod)) # Using MASS package
x <- cbind(mod.table, mod.exp.ci) %>%
rownames_to_column() %>%
mutate_if(is.double, round, 4) %>%
rename(`Exp(B) Lower95CI`=`2.5 %`) %>%
rename(`Exp(B) Upper95CI`=`97.5 %`) %>%
column_to_rownames() %>%
as.matrix()
# mod.summary <- summary(mod)
# mod.summary$coefficients <- x
# mod.summary
pander(x)| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | -6.354 | 0.3701 | -17.17 | 0 |
| TypeKnowledge | 0.3154 | 0.3071 | 1.027 | 0.3044 |
| TypeKnowledge w/ EX | 3.931 | 0.2709 | 14.51 | 0 |
| TypeQuiz | 4.912 | 0.2784 | 17.64 | 0 |
| TypeT/F | 5.029 | 0.2737 | 18.38 | 0 |
| grade11 | -0.0947 | 0.089 | -1.065 | 0.2871 |
| grade12 | -0.0538 | 0.0854 | -0.6303 | 0.5285 |
| grade7 | 0.0782 | 0.1302 | 0.6004 | 0.5482 |
| grade8 | 0.0356 | 0.1804 | 0.1974 | 0.8435 |
| grade9 | -0.095 | 0.0932 | -1.019 | 0.3081 |
| raceAsian | 0.1379 | 0.2794 | 0.4937 | 0.6215 |
| raceBlack | 0.2564 | 0.2511 | 1.021 | 0.3072 |
| raceHispanic | 0.1618 | 0.2426 | 0.6667 | 0.505 |
| raceOther | 0.2364 | 0.2966 | 0.797 | 0.4254 |
| raceWhite | 0.1511 | 0.2419 | 0.6248 | 0.5321 |
| genderMale | -0.0088 | 0.0653 | -0.1347 | 0.8929 |
| genderOther | -0.1408 | 0.3921 | -0.3591 | 0.7195 |
| userResponses | 0.1489 | 0.0036 | 41.76 | 0 |
| campg.num | -0.0379 | 0.0021 | -18.08 | 0 |
| Exp(Estimate) | Exp(B) Lower95CI | Exp(B) Upper95CI | |
|---|---|---|---|
| (Intercept) | 0.0017 | 8e-04 | 0.0035 |
| TypeKnowledge | 1.371 | 0.7659 | 2.572 |
| TypeKnowledge w/ EX | 50.95 | 30.93 | 90.02 |
| TypeQuiz | 135.9 | 81.14 | 243.2 |
| TypeT/F | 152.7 | 92.17 | 271.1 |
| grade11 | 0.9097 | 0.7641 | 1.083 |
| grade12 | 0.9476 | 0.8017 | 1.12 |
| grade7 | 1.081 | 0.8364 | 1.393 |
| grade8 | 1.036 | 0.7232 | 1.468 |
| grade9 | 0.9094 | 0.7574 | 1.091 |
| raceAsian | 1.148 | 0.6644 | 1.988 |
| raceBlack | 1.292 | 0.7918 | 2.121 |
| raceHispanic | 1.176 | 0.7322 | 1.897 |
| raceOther | 1.267 | 0.7072 | 2.264 |
| raceWhite | 1.163 | 0.7257 | 1.875 |
| genderMale | 0.9912 | 0.8718 | 1.126 |
| genderOther | 0.8687 | 0.3736 | 1.771 |
| userResponses | 1.161 | 1.153 | 1.169 |
| campg.num | 0.9628 | 0.9589 | 0.9668 |
rm(mod.table, mod.exp.ci, x)stargazer(mod, type='latex', title="Log-odds",
single.row=T, intercept.top=T, intercept.bottom=F)# stargazer(mod, coef=list(exp(coef(mod))), p.auto=FALSE, type = 'latex',
# single.row=T, intercept.top=T, intercept.bottom=F, title="Exponential Log odds")
stargazer(mod, apply.coef=exp, apply.se=exp, type="latex", p.auto=F,
single.row=T, intercept.top=T, intercept.bottom=F, title="Exponential of Log odds")summary(mod)anova(mod, test="Chisq") #see ?anova.glm
par(mfrow = c(2, 2))
plot(mod)Prediction
Based on modeling only the significant variables from above, which has a higher AIC (second model’s AIC=12,331) than the model above (AIC=8,090). Note that the figure shows three points in time (beginning, midpoint, and end), but assumes userResponses is the same across the board
mod2 <- glm(Replied ~ Type + userResponses + campg.num,
data=df, family=binomial(link = "logit"))
newDf <- data.frame(
userResponses = rep(mean(df$userResponses),15),
campg.num = c(rep(min(df$campg.num),5), rep(mean(df$campg.num),5), rep(max(df$campg.num),5)),
Time = c(rep("Start",5), rep("Midpoint",5), rep("End",5)),
Type = rep(unique(df$Type),3)
)
cbind(newDf, predict(mod2, newdata=newDf, type="response", se.fit=TRUE)) %>%
rename(prob=fit, se.prob=se.fit) %>%
mutate(
ll = prob - 1.96*se.prob,
ul = prob + 1.96*se.prob,
) %>%
ggplot(aes(x=Type, y = prob)) +
geom_errorbar(aes(ymin = ll, ymax = ul, color=Time), width = 0.2, lty=1, lwd=1) +
geom_point(shape=18, size=3, fill="black") +
scale_color_discrete(name="Time during campaign") +
labs(title= " Predicted probabilities", x="Question Type", y="Pr(respond to message)") 
rm(newDf, mod2) The next output below models the fixed effects of variables of interest (e.g. question type) while accounting for random effects (userID). Note: the model doesn’t converge if gender, race, or grade is included, presumably because those effects are accounted for in the random effects of the userID
df <- HRTC %>%
filter(userResponses>0) %>%
mutate(Dummy = as.factor(zipcode)) %>%
mutate(userID = as.factor(userID)) %>%
mutate(Replied = ifelse(Replied==1, "Yes", "No")) %>%
mutate(Replied = factor(Replied, levels=c("No","Yes")))
multiMod <- glmer(Replied ~ Type + (1 | userID),
data=df, family=binomial(link = "logit"))
# control=glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 10000)))
summary(multiMod)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Replied ~ Type + (1 | userID)
## Data: df
##
## AIC BIC logLik deviance df.resid
## 12952.6 12999.4 -6470.3 12940.6 18068
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1588 -0.3957 -0.1975 -0.0403 24.7967
##
## Random effects:
## Groups Name Variance Std.Dev.
## userID (Intercept) 2.119 1.456
## Number of obs: 18074, groups: userID, 414
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.3448 0.2489 -21.474 <2e-16 ***
## TypeKnowledge 0.5810 0.2626 2.212 0.0269 *
## TypeKnowledge w/ EX 3.4863 0.2392 14.572 <2e-16 ***
## TypeQuiz 4.2730 0.2426 17.613 <2e-16 ***
## TypeT/F 4.7069 0.2411 19.524 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TypKnw TKw/EX TypeQz
## TypeKnowldg -0.842
## TypKnwlw/EX -0.942 0.876
## TypeQuiz -0.934 0.864 0.967
## TypeT/F -0.943 0.870 0.974 0.966# Basic AIC = 11,174
# Gender AIC = 8,037
# Grade AIC = 7,861 (no convergence)
# Race AIC = 8,105 (no convergence)
# Gender grade AIC = 7,574 (no convergence)plot(multiMod)
lattice::qqmath(multiMod)# Consider using broom::augment.glm()
newDf <- data.frame(
userID = unique(df$userID),
Type = rep(unique(df$Type),414)
)
cbind(newDf, prob=predict(multiMod, newdata=newDf)) %>%
mutate(prob= exp(prob)) %>%
filter(prob<1) %>%
ggplot(aes(x=Type, y=prob)) + geom_boxplot(varwidth=T)Some resources for survival analysis (for my reference):
Gender
df <- users %>%
filter(gender!="Other")
survObj <- Surv(time=df$SurvEnd2, event=df$Censor)
survfit(survObj ~ gender, data = df) %>%
ggsurvplot(data = df, pval=T, #conf.int=T,
risk.table = TRUE,
fun = "pct",
size = 1,
legend = "bottom")
# ggsurv() +
# scale_y_continuous(labels = scales::percent) +
# scale_x_continuous(name="Days into campaign") +
# scale_color_brewer(palette="Set2", na.value = "grey50")Grade
df <- users %>%
filter(grade!="other")
survObj <- Surv(time=df$SurvEnd2, event=df$Censor)
survfit(survObj ~ grade, data = df) %>%
ggsurvplot(data = df, pval=T, #conf.int=T,
risk.table = TRUE,
fun = "pct",
size = 1,
legend = "bottom")
Race
df <- users
survObj <- Surv(time=df$SurvEnd2, event=df$Censor)
survfit(survObj ~ race, data = df) %>%
ggsurvplot(data = df, pval=T, #conf.int=T,
risk.table = TRUE,
fun = "pct",
size = 1,
legend = "bottom")
Combo
df <- users %>%
filter(grade!="other", gender!="Other")
survObj <- Surv(time=df$SurvEnd2, event=df$Censor)
# Fit a Cox proportional hazards model
coxph(survObj ~ grade + race, data=df) %>% summary()# Export to SPSS
# HRTC %>% filter(userResponses>1) %>% select(userID, message_id, Type, mid, Replied:userResponses) %>% write_csv("~/Downloads/data.csv")So there’s quite a few NA’s (people who didn’t have their gender, grade, or race recorded). A few people did respond to the gender & grade questions, but didn’t fit into a larger category, so I dropped them from the analysis (total of 7 students)
Total number of students: 525
Gender
users %>%
semi_join(count(HRTC, userID)) %>%
count(gender, sort=T) %>%
mutate(sansNA = ifelse(is.na(gender),0,n)) %>%
mutate(gender = ifelse(is.na(gender),"NA",as.character(gender))) %>%
filter(gender!="Other") %>%
mutate(sansNA = percent(sansNA/sum(sansNA))) %>%
mutate(Percent = percent(n/sum(n))) %>%
# select(gender, n, Percent, sansNA) %>%
kable()| gender | n | sansNA | Percent |
|---|---|---|---|
| Female | 238 | 72.1% | 45.9% |
| NA | 189 | 0.0% | 36.4% |
| Male | 92 | 27.9% | 17.7% |
Grade
users %>%
semi_join(count(HRTC, userID)) %>%
count(grade, sort=T) %>%
mutate(sansNA = ifelse(is.na(grade),0,n)) %>%
mutate(grade = ifelse(is.na(grade),"NA",as.character(grade))) %>%
filter(grade!="other") %>%
mutate(sansNA = percent(sansNA/sum(sansNA))) %>%
mutate(Percent = percent(n/sum(n))) %>%
kable()| grade | n | sansNA | Percent |
|---|---|---|---|
| NA | 207 | 0.0% | 39.5% |
| 12 | 99 | 31.2% | 18.9% |
| 10 | 64 | 20.2% | 12.2% |
| 11 | 58 | 18.3% | 11.1% |
| 9 | 58 | 18.3% | 11.1% |
| 7 | 25 | 7.9% | 4.8% |
| 8 | 13 | 4.1% | 2.5% |
Race
Note that unlike gender & grade, I didn’t drop the Other race because it’s representing Two or more races
users %>%
semi_join(count(HRTC, userID)) %>%
count(race, sort=T) %>%
mutate(sansNA = ifelse(is.na(race),0,n)) %>%
mutate(race = ifelse(is.na(race),"NA",as.character(race))) %>%
mutate(sansNA = percent(sansNA/sum(sansNA))) %>%
mutate(Percent = percent(n/sum(n))) %>%
kable()| race | n | sansNA | Percent |
|---|---|---|---|
| NA | 189 | 0.0% | 36.0% |
| White | 157 | 46.7% | 29.9% |
| Hispanic | 96 | 28.6% | 18.3% |
| Black | 50 | 14.9% | 9.5% |
| Asian | 17 | 5.1% | 3.2% |
| Other | 12 | 3.6% | 2.3% |
| AmInd/AkN | 4 | 1.2% | 0.8% |
Compare to Census data (2017 ACS 5-year estimates) from Galveston County:
data_frame(
Race = c("White (non-Hispanic)", "Hispanic or Latino",
"Black or African American", "Asian", "Two or more races",
"American Indian and Alaska Native"),
Percent = c(0.459, 0.287, 0.208, 0.032, 0.025, 0.005)
) %>%
mutate(Percent = percent(Percent)) %>%
# arrange(desc(Percent))
kable()| Race | Percent |
|---|---|
| White (non-Hispanic) | 45.9% |
| Hispanic or Latino | 28.7% |
| Black or African American | 20.8% |
| Asian | 3.2% |
| Two or more races | 2.5% |
| American Indian and Alaska Native | 0.5% |
Session info
pander(sessionInfo())R version 3.3.2 (2016-10-31)
**Platform:** x86_64-apple-darwin13.4.0 (64-bit)
locale: en_US.UTF-8||en_US.UTF-8||en_US.UTF-8||C||en_US.UTF-8||en_US.UTF-8
attached base packages: stats, graphics, grDevices, utils, datasets, methods and base
other attached packages: bindrcpp(v.0.2.2), pander(v.0.6.3), survMisc(v.0.5.5), survminer(v.0.4.3), ggpubr(v.0.2), magrittr(v.1.5), survival(v.2.40-1), lme4(v.1.1-19), Matrix(v.1.2-7.1), stargazer(v.5.2.2), scales(v.0.5.0), knitr(v.1.18), lubridate(v.1.7.1), forcats(v.0.2.0), stringr(v.1.2.0), dplyr(v.0.7.8), purrr(v.0.2.4), readr(v.1.1.1), tidyr(v.0.7.2), tibble(v.1.4.2), ggplot2(v.3.1.0.9000), tidyverse(v.1.2.1) and MASS(v.7.3-45)
loaded via a namespace (and not attached): TH.data(v.1.0-10), minqa(v.1.2.4), colorspace(v.1.3-2), modeltools(v.0.2-21), rprojroot(v.1.3-2), estimability(v.1.3), htmlTable(v.1.13.1), base64enc(v.0.1-3), rstudioapi(v.0.7), mice(v.2.46.0), fifer(v.1.1), mvtnorm(v.1.0-5), coin(v.1.2-2), xml2(v.1.1.1), codetools(v.0.2-15), splines(v.3.3.2), spam(v.2.1-2), Formula(v.1.2-3), jsonlite(v.1.5), nloptr(v.1.2.1), broom(v.0.5.1.9000), km.ci(v.0.5-2), cluster(v.2.0.5), httr(v.1.3.1), randomForestSRC(v.2.5.1), lsmeans(v.2.30-0), emmeans(v.1.3.2), backports(v.1.1.2), assertthat(v.0.2.0), lazyeval(v.0.2.1), strucchange(v.1.5-1), cli(v.1.0.1), acepack(v.1.4.1), htmltools(v.0.3.6), tools(v.3.3.2), dotCall64(v.0.9-5), gtable(v.0.2.0), glue(v.1.2.0), ggthemes(v.3.3.0), maps(v.3.1.1), Rcpp(v.1.0.0), cellranger(v.1.1.0), nlme(v.3.1-128), rvest(v.0.3.2), zoo(v.1.8-4), hms(v.0.4.0), parallel(v.3.3.2), sandwich(v.2.5-0), RColorBrewer(v.1.1-2), fields(v.9.0), yaml(v.2.1.16), gridExtra(v.2.2.1), KMsurv(v.0.1-5), rpart(v.4.1-10), latticeExtra(v.0.6-28), stringi(v.1.1.6), highr(v.0.6), plotrix(v.3.7-4), randomForest(v.4.6-12), checkmate(v.1.9.1), rlang(v.0.3.0.1), pkgconfig(v.2.0.2), evaluate(v.0.10), lattice(v.0.20-34), bindr(v.0.1.1), htmlwidgets(v.1.3), labeling(v.0.3), cmprsk(v.2.2-7), cowplot(v.0.9.3), tidyselect(v.0.2.3), plyr(v.1.8.4), R6(v.2.2.2), generics(v.0.0.2), Hmisc(v.4.1-1), multcomp(v.1.4-8), pillar(v.1.3.1), haven(v.1.1.0), foreign(v.0.8-67), withr(v.2.1.2), nnet(v.7.3-12), modelr(v.0.1.1), crayon(v.1.3.4), party(v.1.2-4), rmarkdown(v.1.8), grid(v.3.3.2), readxl(v.1.0.0), data.table(v.1.10.4-3), digest(v.0.6.13), xtable(v.1.8-2), stats4(v.3.3.2) and munsell(v.0.5.0)