Mental Health Patient Improvement - Ordinal Logistic Regression
Mackenzie Zendt
9/9/2020
need to put consult number into bins!!!
Categorize Depression Score (Outcome)
We will divide “last PHQ9 score” into three categories - people who have no depression (scores between 0 and 5), people who have some depression (scores between 6 and 10), and people who still have depression (scores 10 and above). Below, we print out a table to ensure our coding worked correctly.
communities$none_mild_some <- cut(communities$lastphq9, breaks = c(0, 5, 10, 28),
labels = c("None", "Mild", "Some"))
table(communities$lastphq9, communities$none_mild_some) %>% knitr::kable() #check to see it worked```| None | Mild | Some | |
|---|---|---|---|
| 0 | 0 | 0 | 0 |
| 1 | 10 | 0 | 0 |
| 2 | 10 | 0 | 0 |
| 3 | 14 | 0 | 0 |
| 4 | 18 | 0 | 0 |
| 5 | 10 | 0 | 0 |
| 6 | 0 | 20 | 0 |
| 7 | 0 | 14 | 0 |
| 8 | 0 | 14 | 0 |
| 9 | 0 | 18 | 0 |
| 10 | 0 | 23 | 0 |
| 11 | 0 | 0 | 28 |
| 12 | 0 | 0 | 29 |
| 13 | 0 | 0 | 23 |
| 14 | 0 | 0 | 20 |
| 15 | 0 | 0 | 24 |
| 16 | 0 | 0 | 17 |
| 17 | 0 | 0 | 19 |
| 18 | 0 | 0 | 9 |
| 19 | 0 | 0 | 12 |
| 20 | 0 | 0 | 13 |
| 21 | 0 | 0 | 12 |
| 22 | 0 | 0 | 12 |
| 23 | 0 | 0 | 9 |
| 24 | 0 | 0 | 3 |
| 25 | 0 | 0 | 9 |
| 27 | 0 | 0 | 1 |
Independent Variables: Age Category, Gender, and Number of Consults
# with categorical variables
with(communities, table(gender, none_mild_some)) %>% knitr::kable()| None | Mild | Some | |
|---|---|---|---|
| F | 51 | 75 | 204 |
| M | 6 | 9 | 29 |
with(communities, table(age_cat, none_mild_some)) %>% knitr::kable()| None | Mild | Some | |
|---|---|---|---|
| 0-18 | 1 | 3 | 15 |
| 18-34 | 19 | 29 | 102 |
| 35-49 | 23 | 34 | 78 |
| 40-64 | 16 | 17 | 28 |
| 65+ | 3 | 6 | 15 |
with(communities, table(consult_cat, none_mild_some)) %>% knitr::kable()| None | Mild | Some | |
|---|---|---|---|
| 1 | 7 | 25 | 78 |
| 10-20 | 14 | 11 | 20 |
| 2-4 | 18 | 28 | 83 |
| 20+ | 9 | 7 | 12 |
| 5-7 | 10 | 11 | 32 |
| 8-10 | 4 | 7 | 15 |
Ordinal Logistic Regression Model
Describe the model in words.
# Regression Model ####
#choose a baseline (one of the outcomes) and set it with the relevel funcion - choosing "still has depression"
communities$none_mild_some1 <- relevel(communities$none_mild_some, ref = "Some")
# then run the regression using multinom
test <- multinom(none_mild_some1 ~ age_cat + gender + consult.number, data = communities)## # weights: 24 (14 variable)
## initial value 408.683771
## iter 10 value 328.612607
## iter 20 value 327.735670
## final value 327.735490
## converged# and ask for a summary
summary(test) ## Call:
## multinom(formula = none_mild_some1 ~ age_cat + gender + consult.number,
## data = communities)
##
## Coefficients:
## (Intercept) age_cat18-34 age_cat35-49 age_cat40-64 age_cat65+ genderM
## None -2.915668 0.6727675 1.263048 1.740369 0.819362 0.10970110
## Mild -2.070389 0.6654353 1.092045 1.355439 1.103736 -0.04719274
## consult.number
## None 0.06132512
## Mild 0.02432714
##
## Std. Errors:
## (Intercept) age_cat18-34 age_cat35-49 age_cat40-64 age_cat65+ genderM
## None 1.0448178 1.0738989 1.0682127 1.0948216 1.2320797 0.4922643
## Mild 0.7631977 0.7863392 0.7867136 0.8222901 0.9074467 0.4138577
## consult.number
## None 0.01733581
## Mild 0.01767899
##
## Residual Deviance: 655.471
## AIC: 683.471Explanation:
Add p-values
z <- summary(test)$coefficients/summary(test)$standard.errors
# z - could print out z, but not necessary
# here are the pvalues
# 2-tailed z test
p <- (1 - pnorm(abs(z), 0, 1)) * 2
p## (Intercept) age_cat18-34 age_cat35-49 age_cat40-64 age_cat65+ genderM
## None 0.005261048 0.5310054 0.2370494 0.11191655 0.5060355 0.8236523
## Mild 0.006672098 0.3974163 0.1651037 0.09927661 0.2238669 0.9092129
## consult.number
## None 0.0004039611
## Mild 0.1688068986When looking at the p-values, we see that most do not reach a 0.05 level of statistical significance. The only one that does is consult number, telling us that there is a statistically significant difference in the “no depression” group (compared to “some depression” group, the baseline) based on the number of clinic appointments attended.
Patients aged 50-64 have a lower alpha level, as do patients who are female, but these are not statistically significant at alpha=0.05.
Risk Ratios
We can extract the coefficients from the model and exponentiate them to get risk ratios, shown below. These can be easier to interpret than the raw coefficient, which tells us the log odds.
exp(coef(test))## (Intercept) age_cat18-34 age_cat35-49 age_cat40-64 age_cat65+ genderM
## None 0.05416782 1.959653 3.536184 5.699448 2.269052 1.1159445
## Mild 0.12613671 1.945337 2.980362 3.878463 3.015412 0.9539035
## consult.number
## None 1.063245
## Mild 1.024625Explanation:
Create predicted probabilities based on the model
We exclude gender from this excercise - we are using one categorized variable (age categories) and one continous (consult number).
# first make a new model without gender
test2 <- multinom(none_mild_some1 ~ age_cat + consult.number, data = communities)## # weights: 21 (12 variable)
## initial value 427.360180
## iter 10 value 349.854007
## final value 348.847055
## converged# Now use the predicted probabilities to look at the averaged predicted probabilities
# for different values of the continuous predictor variable consult within each level of age_cat
dconsult <- data.frame(age_cat = rep(c("18-34", "35-49", "40-64", "65+"), each = 51), consult.number = rep(c(0:50), 4))
## store the predicted probabilities for each value of age and consult
pp.consult <- cbind(dconsult, predict(test2, newdata = dconsult, type = "probs", se = TRUE))
## calculate the mean probabilities within each level of age_cat
by(pp.consult[, 3:5], pp.consult$age_cat, colMeans)## pp.consult$age_cat: 18-34
## Some None Mild
## 0.5178325 0.2829233 0.1992441
## ------------------------------------------------------------
## pp.consult$age_cat: 35-49
## Some None Mild
## 0.4243370 0.3316632 0.2439998
## ------------------------------------------------------------
## pp.consult$age_cat: 40-64
## Some None Mild
## 0.3230017 0.4257870 0.2512113
## ------------------------------------------------------------
## pp.consult$age_cat: 65+
## Some None Mild
## 0.4935977 0.2493422 0.2570601# calculate the mean probabilities within each level of consult number
by(pp.consult[, 3:5], pp.consult$consult.number, colMeans)## pp.consult$consult.number: 0
## Some None Mild
## 0.6384736 0.1222375 0.2392889
## ------------------------------------------------------------
## pp.consult$consult.number: 1
## Some None Mild
## 0.6318666 0.1273965 0.2407369
## ------------------------------------------------------------
## pp.consult$consult.number: 2
## Some None Mild
## 0.6251476 0.1327317 0.2421207
## ------------------------------------------------------------
## pp.consult$consult.number: 3
## Some None Mild
## 0.6183167 0.1382459 0.2434374
## ------------------------------------------------------------
## pp.consult$consult.number: 4
## Some None Mild
## 0.6113743 0.1439419 0.2446838
## ------------------------------------------------------------
## pp.consult$consult.number: 5
## Some None Mild
## 0.6043211 0.1498221 0.2458568
## ------------------------------------------------------------
## pp.consult$consult.number: 6
## Some None Mild
## 0.5971578 0.1558888 0.2469534
## ------------------------------------------------------------
## pp.consult$consult.number: 7
## Some None Mild
## 0.5898856 0.1621439 0.2479705
## ------------------------------------------------------------
## pp.consult$consult.number: 8
## Some None Mild
## 0.5825057 0.1685892 0.2489051
## ------------------------------------------------------------
## pp.consult$consult.number: 9
## Some None Mild
## 0.5750197 0.1752260 0.2497543
## ------------------------------------------------------------
## pp.consult$consult.number: 10
## Some None Mild
## 0.5674293 0.1820556 0.2505150
## ------------------------------------------------------------
## pp.consult$consult.number: 11
## Some None Mild
## 0.5597367 0.1890787 0.2511846
## ------------------------------------------------------------
## pp.consult$consult.number: 12
## Some None Mild
## 0.5519441 0.1962957 0.2517602
## ------------------------------------------------------------
## pp.consult$consult.number: 13
## Some None Mild
## 0.5440542 0.2037066 0.2522392
## ------------------------------------------------------------
## pp.consult$consult.number: 14
## Some None Mild
## 0.5360699 0.2113111 0.2526190
## ------------------------------------------------------------
## pp.consult$consult.number: 15
## Some None Mild
## 0.5279942 0.2191086 0.2528972
## ------------------------------------------------------------
## pp.consult$consult.number: 16
## Some None Mild
## 0.5198307 0.2270977 0.2530716
## ------------------------------------------------------------
## pp.consult$consult.number: 17
## Some None Mild
## 0.5115832 0.2352769 0.2531399
## ------------------------------------------------------------
## pp.consult$consult.number: 18
## Some None Mild
## 0.5032555 0.2436441 0.2531003
## ------------------------------------------------------------
## pp.consult$consult.number: 19
## Some None Mild
## 0.4948521 0.2521969 0.2529509
## ------------------------------------------------------------
## pp.consult$consult.number: 20
## Some None Mild
## 0.4863776 0.2609323 0.2526901
## ------------------------------------------------------------
## pp.consult$consult.number: 21
## Some None Mild
## 0.4778368 0.2698467 0.2523165
## ------------------------------------------------------------
## pp.consult$consult.number: 22
## Some None Mild
## 0.4692349 0.2789362 0.2518288
## ------------------------------------------------------------
## pp.consult$consult.number: 23
## Some None Mild
## 0.4605774 0.2881965 0.2512262
## ------------------------------------------------------------
## pp.consult$consult.number: 24
## Some None Mild
## 0.4518698 0.2976225 0.2505077
## ------------------------------------------------------------
## pp.consult$consult.number: 25
## Some None Mild
## 0.4431181 0.3072089 0.2496730
## ------------------------------------------------------------
## pp.consult$consult.number: 26
## Some None Mild
## 0.4343285 0.3169498 0.2487216
## ------------------------------------------------------------
## pp.consult$consult.number: 27
## Some None Mild
## 0.4255074 0.3268389 0.2476537
## ------------------------------------------------------------
## pp.consult$consult.number: 28
## Some None Mild
## 0.4166613 0.3368692 0.2464695
## ------------------------------------------------------------
## pp.consult$consult.number: 29
## Some None Mild
## 0.4077970 0.3470337 0.2451693
## ------------------------------------------------------------
## pp.consult$consult.number: 30
## Some None Mild
## 0.3989214 0.3573245 0.2437542
## ------------------------------------------------------------
## pp.consult$consult.number: 31
## Some None Mild
## 0.3900416 0.3677335 0.2422249
## ------------------------------------------------------------
## pp.consult$consult.number: 32
## Some None Mild
## 0.3811647 0.3782523 0.2405830
## ------------------------------------------------------------
## pp.consult$consult.number: 33
## Some None Mild
## 0.3722981 0.3888720 0.2388299
## ------------------------------------------------------------
## pp.consult$consult.number: 34
## Some None Mild
## 0.3634491 0.3995833 0.2369676
## ------------------------------------------------------------
## pp.consult$consult.number: 35
## Some None Mild
## 0.3546251 0.4103769 0.2349981
## ------------------------------------------------------------
## pp.consult$consult.number: 36
## Some None Mild
## 0.3458335 0.4212427 0.2329238
## ------------------------------------------------------------
## pp.consult$consult.number: 37
## Some None Mild
## 0.3370818 0.4321708 0.2307474
## ------------------------------------------------------------
## pp.consult$consult.number: 38
## Some None Mild
## 0.3283773 0.4431510 0.2284717
## ------------------------------------------------------------
## pp.consult$consult.number: 39
## Some None Mild
## 0.3197274 0.4541727 0.2260999
## ------------------------------------------------------------
## pp.consult$consult.number: 40
## Some None Mild
## 0.3111393 0.4652255 0.2236352
## ------------------------------------------------------------
## pp.consult$consult.number: 41
## Some None Mild
## 0.3026201 0.4762986 0.2210813
## ------------------------------------------------------------
## pp.consult$consult.number: 42
## Some None Mild
## 0.2941769 0.4873813 0.2184418
## ------------------------------------------------------------
## pp.consult$consult.number: 43
## Some None Mild
## 0.2858164 0.4984628 0.2157208
## ------------------------------------------------------------
## pp.consult$consult.number: 44
## Some None Mild
## 0.2775453 0.5095325 0.2129222
## ------------------------------------------------------------
## pp.consult$consult.number: 45
## Some None Mild
## 0.2693699 0.5205796 0.2100505
## ------------------------------------------------------------
## pp.consult$consult.number: 46
## Some None Mild
## 0.2612965 0.5315937 0.2071099
## ------------------------------------------------------------
## pp.consult$consult.number: 47
## Some None Mild
## 0.2533309 0.5425642 0.2041049
## ------------------------------------------------------------
## pp.consult$consult.number: 48
## Some None Mild
## 0.2454788 0.5534810 0.2010402
## ------------------------------------------------------------
## pp.consult$consult.number: 49
## Some None Mild
## 0.2377455 0.5643340 0.1979205
## ------------------------------------------------------------
## pp.consult$consult.number: 50
## Some None Mild
## 0.2301361 0.5751134 0.1947505We see that there is not much difference in outcome by age category. However, within consult number, as number of consults increases, a person is more likely to be in the “no depression” category.
Plotting
## melt data set to long for ggplot2
plot1data <- melt(pp.consult, id.vars = c("age_cat", "consult.number"), value.name = "probability")
head(plot1data) # view first few rows## age_cat consult.number variable probability
## 1 18-34 0 Some 0.7179492
## 2 18-34 1 Some 0.7119687
## 3 18-34 2 Some 0.7058530
## 4 18-34 3 Some 0.6996005
## 5 18-34 4 Some 0.6932097
## 6 18-34 5 Some 0.6866794## plot predicted probabilities across write values for each level of ses
## facetted by program type
ggplot(plot1data, aes(x = consult.number, y = probability, colour = age_cat)) + geom_line() + facet_grid(variable ~ ., scales = "free")
This graph shows us that as consult number increases, the probability of being in “none” increases as well. Most significantly for patients between the ages of 50-64 and 35-49.