Good Friday
Simon Analysis
2024-03-29
[1] "C:/Users/PC/Downloads"
Top 5 Occupations (with Percentages) and Totals
===================================
occup n total Percentage
-----------------------------------
business 61 220 27.7%
student 38 220 17.3%
farmer 25 220 11.4%
housewife 17 220 7.7%
casual labourer 15 220 6.8%
-----------------------------------
Top 5 Occupations by Time and Treatment Control (with Percentages)
===========================================================
time treat_control occup n total Percentage
-----------------------------------------------------------
baseline control business 14 46 30.4%
baseline control student 8 46 17.4%
baseline control farmer 6 46 13.0%
baseline control casual labourer 5 46 10.9%
baseline control housewife 4 46 8.7%
baseline treatment business 13 57 22.8%
baseline treatment student 12 57 21.1%
baseline treatment farmer 10 57 17.5%
baseline treatment housewife 7 57 12.3%
baseline treatment casual labourer 4 57 7.0%
baseline treatment security officer 4 57 7.0%
baseline treatment teacher 4 57 7.0%
endline control business 19 64 29.7%
endline control student 7 64 10.9%
endline control housewife 5 64 7.8%
endline control casual labourer 4 64 6.2%
endline control farmer 4 64 6.2%
endline treatment business 15 53 28.3%
endline treatment student 11 53 20.8%
endline treatment farmer 5 53 9.4%
endline treatment teacher 3 53 5.7%
endline treatment casual labourer 2 53 3.8%
endline treatment security officer 2 53 3.8%
-----------------------------------------------------------
Participants' gender by Time and Treatment Control (with Percentages)
=================================================
time treat_control gender n total Percentage
-------------------------------------------------
baseline control female 20 46 43.5%
baseline control male 26 46 56.5%
baseline treatment female 24 57 42.1%
baseline treatment male 33 57 57.9%
endline control female 27 64 42.2%
endline control male 37 64 57.8%
endline treatment female 20 53 37.7%
endline treatment male 33 53 62.3%
-------------------------------------------------
Participants' marital status by Time and Treatment Control (with Percentages)
==================================================
time treat_control mstatus n total Percentage
--------------------------------------------------
baseline control married 30 46 65.2%
baseline control single 16 46 34.8%
baseline treatment married 39 57 68.4%
baseline treatment single 17 57 29.8%
baseline treatment widowed 1 57 1.8%
endline control married 53 64 82.8%
endline control single 11 64 17.2%
endline treatment married 37 53 69.8%
endline treatment single 15 53 28.3%
endline treatment widowed 1 53 1.9%
--------------------------------------------------
Participants' education level by Time and Treatment Control (with Percentages)
=====================================================
time treat_control edulevel n total Percentage
-----------------------------------------------------
baseline control college 8 46 17.4%
baseline control primary 18 46 39.1%
baseline control secondary 20 46 43.5%
baseline treatment college 13 57 22.8%
baseline treatment primary 25 57 43.9%
baseline treatment secondary 19 57 33.3%
endline control college 11 64 17.2%
endline control primary 16 64 25.0%
endline control secondary 34 64 53.1%
endline control university 3 64 4.7%
endline treatment college 12 53 22.6%
endline treatment primary 10 53 18.9%
endline treatment secondary 31 53 58.5%
-----------------------------------------------------
Participants' county by Time and Treatment Control (with Percentages)
===================================================
time treat_control county n total Percentage
---------------------------------------------------
baseline control migori 26 46 56.5%
baseline control siaya 20 46 43.5%
baseline treatment kakamega 34 57 59.6%
baseline treatment vihiga 23 57 40.4%
endline control migori 32 64 50.0%
endline control siaya 32 64 50.0%
endline treatment kakamega 27 53 50.9%
endline treatment vihiga 26 53 49.1%
---------------------------------------------------
Participants' age by Time and Treatment Control
======================================================================
time treat_control mean_duration ci_lower ci_upper min_age max_age
----------------------------------------------------------------------
baseline control 38.957 34.119 43.794 18 79
baseline treatment 39.982 35.611 44.354 17 79
endline control 41.906 38.658 45.154 21 67
endline treatment 39.340 35.132 43.547 18 72
----------------------------------------------------------------------
========================
MOI count
------------------------
1 rti 111
2 fall 32
3 stab 29
4 assault 18
5 burns 14
6 gunshot 5
7 cut 4
8 rti, motorbike 2
9 crushing 1
10 penetrating inj 1
11 rti,fall 1
12 stab, assault 1
13 stab/ cut 1
------------------------
=================================
time MOI count
---------------------------------
1 endline rti 58
2 baseline rti 53
3 endline fall 17
4 baseline fall 15
5 endline stab 15
6 baseline stab 14
7 baseline assault 10
8 endline assault 8
9 endline burns 8
10 baseline burns 6
11 endline gunshot 5
12 endline cut 3
13 baseline rti, motorbike 2
14 baseline cut 1
15 baseline rti,fall 1
16 baseline stab, assault 1
17 endline crushing 1
18 endline penetrating inj 1
19 endline stab/ cut 1
---------------------------------
==================================================================================
time treat_control mean_duration ci_lower ci_upper min_duration max_duration
----------------------------------------------------------------------------------
1 baseline control 68.370 50.265 86.474 30 360
2 baseline treatment 55.702 47.135 64.269 15 210
3 endline control 64.375 58.525 70.225 30 130
4 endline treatment 34.113 30.378 37.849 10 65
----------------------------------------------------------------------------------# Create new variable Occupational Therapy Admissions
survey_data1$OT_admissions <- ifelse(survey_data1$Admit == "icu", 1, 0)package 'tableone' successfully unpacked and MD5 sums checked
The downloaded binary packages are in
C:\Users\PC\AppData\Local\Temp\Rtmp2pV6Zk\downloaded_packages Stratified by treat_control
level control treatment p test
n 110 110
ICU_admissions (mean (SD)) 0.19 (0.39) 0.21 (0.41) 0.737
HRinitial (mean (SD)) 95.75 (14.61) 99.78 (15.41) 0.048
SBPinitial (mean (SD)) 109.61 (20.11) 104.94 (21.36) 0.096 
colnames(survey_data1)## [1] "code" "county" "time"
## [4] "treat_control" "age" "gender"
## [7] "mstatus" "edulevel" "occup"
## [10] "dateinjury" "timeinj" "MOI"
## [13] "anatlocation" "ISS" "admdate"
## [16] "arrivaltime" "arrivalmode" "refby"
## [19] "receiveby" "triagecategory" "timeassmnt"
## [22] "Durationmins" "state" "HRinitial"
## [25] "SBPinitial" "DBPinitial" "RRinitial"
## [28] "SPO2initial" "painscore" "GCS"
## [31] "AVPU" "airwaystatus" "airwayinterven"
## [34] "spinestabilized" "breathingstatus" "breathinginterven"
## [37] "bleedingcontrol" "Fluidadmin" "meds"
## [40] "Eddeparture" "Edduration_mins" "HRdispo"
## [43] "SBPdispo" "DBPdispo" "RRdispo"
## [46] "SPO2dispo" "Admit" "hospdischarge"
## [49] "LOSdays" "state.1" "outcome"
## [52] "ICU_LOS" "ward_LOS" "ICU_admissions"
## [55] "OT_admissions"# Load necessary package
library(base)
# Function to convert mixed date formats
convert_to_date <- function(date_column) {
# Replace "-" with "/" for consistency
date_column <- gsub("-", "/", date_column)
# Convert to Date format assuming day comes first
date_column <- as.Date(date_column, format = "%d/%m/%Y")
# Return as Date type
return(date_column)
}
# Apply function to date columns
survey_data1$dateinjury <- convert_to_date(survey_data1$dateinjury)
survey_data1$admdate <- convert_to_date(survey_data1$admdate)
glimpse(survey_data1)## Rows: 220
## Columns: 55
## $ code <int> 41, 57, 4, 92, 94, 16, 26, 61, 52, 62, 58, 14, 49, 7…
## $ county <chr> "siaya", "migori", "migori", "siaya", "siaya", "migo…
## $ time <chr> "endline", "baseline", "endline", "baseline", "basel…
## $ treat_control <chr> "control", "control", "control", "control", "control…
## $ age <int> 57, 20, 28, 22, 60, 39, 31, 28, 44, 18, 29, 53, 30, …
## $ gender <chr> "male", "male", "female", "female", "female", "male"…
## $ mstatus <chr> "married", "single", "single", "single", "married", …
## $ edulevel <chr> "secondary", "secondary", "college", "college", "pri…
## $ occup <chr> "security officer", "student", "business", "student"…
## $ dateinjury <date> NA, 0016-04-20, NA, NA, 0018-05-20, 0020-12-20, 002…
## $ timeinj <chr> "0900hrs", "1100hrs", "0900hrs", "1500hrs", "1900hrs…
## $ MOI <chr> "rti", "rti,fall", "stab", "burns", "fall", "rti", "…
## $ anatlocation <chr> "head", "head, face", "chest", "face, chest", "extre…
## $ ISS <int> 15, 16, 15, 16, 15, 15, 15, 16, 15, 16, 15, 15, 16, …
## $ admdate <date> NA, 0016-04-20, NA, NA, 0018-05-20, 0020-12-20, 002…
## $ arrivaltime <chr> "0950hrs", "1300hrs", "1000hrs", "1530hrs", "1930hrs…
## $ arrivalmode <chr> "motorbike", "motorbike", "taxi", "taxi", "motorbike…
## $ refby <chr> "self", "self", "nurse", "self", "self", "self", "se…
## $ receiveby <chr> "mo", "nurse", "mo", "nurse", "rco", "nurse", "mo", …
## $ triagecategory <chr> "immediate", "immediate", "immediate", "urgent", "ur…
## $ timeassmnt <chr> "1020hrs", "1330hrs", "1025hrs", "1600hrs", "2000hrs…
## $ Durationmins <int> 30, 30, 25, 30, 30, 30, 30, 20, 20, 10, 20, 40, 30, …
## $ state <chr> "alive", "alive", "alive", "alive", "alive", "alive"…
## $ HRinitial <int> 98, 94, 90, 105, 98, 102, 122, 110, 98, 118, 88, 99,…
## $ SBPinitial <int> 110, 124, 130, 90, 120, 90, 90, 115, 124, 90, 114, 9…
## $ DBPinitial <int> 80, 72, 80, 70, 70, 60, 75, 68, 75, 65, 78, 65, 76, …
## $ RRinitial <int> 14, 22, 20, 22, 13, 12, 24, 20, 13, 28, 18, 12, 22, …
## $ SPO2initial <int> 93, 94, 95, 88, 92, 90, 82, 95, 94, 90, 94, 84, 82, …
## $ painscore <int> NA, 10, 9, 10, 10, NA, NA, 10, 9, NA, 7, NA, NA, NA,…
## $ GCS <int> 13, 12, 15, 12, 14, 12, 13, 15, 13, 12, 15, 12, 13, …
## $ AVPU <chr> "u", "p", "a", "v", "a", "u", "u", "a", "p", "p", "a…
## $ airwaystatus <chr> "patent", "obstructed", "patent", "obstructed", "pat…
## $ airwayinterven <chr> "none", "none", "none", "suction", "none", "none", "…
## $ spinestabilized <chr> "no", "no", "no", "no", "no", "no", "no", "no", "no"…
## $ breathingstatus <chr> "spont", "spont", "spont", "supplemental", "spont", …
## $ breathinginterven <chr> "none", "none", "none", "ventilator", "none", "none"…
## $ bleedingcontrol <chr> "bandage", "none", "bandage", "N/A", "none", "pressu…
## $ Fluidadmin <chr> "NS 1L", "mannitol", "NS 500mls", "NS 1L", "NS 1L", …
## $ meds <chr> "mannitol", "tramadol", "analgesics", "tramadol", "p…
## $ Eddeparture <chr> "1100hrs", "1600hrs", "1110hrs", "1635hrs", "2040hrs…
## $ Edduration_mins <int> 40, 150, 45, 35, 40, 40, 45, 360, 90, 80, 120, 80, 6…
## $ HRdispo <int> 95, 96, 87, 100, 76, 108, 116, 69, 80, 105, 92, 94, …
## $ SBPdispo <int> 120, 103, 127, 105, 150, 90, 75, 120, 108, 100, 110,…
## $ DBPdispo <int> 70, 70, 69, 60, 60, 75, 50, 70, 75, 75, 80, 75, 50, …
## $ RRdispo <int> 15, 14, 19, 28, 12, 11, 12, 22, 18, 20, 15, 13, 10, …
## $ SPO2dispo <int> 93, 94, 94, 88, 92, 90, 83, 95, 94, 94, 95, 94, 80, …
## $ Admit <chr> "ward", "ward", "ward", "icu", "OT", "icu", "icu", "…
## $ hospdischarge <chr> "15-11-20", "7/5/2020", "23-09-20", "17-07-20", "13-…
## $ LOSdays <int> 35, 21, 15, 15, 25, 30, 1, 25, 20, 22, 24, 26, 2, 13…
## $ state.1 <chr> "fair", "stable", "stable", "rehab", "stable", "stab…
## $ outcome <chr> "alive", "alive", "alive", "alive", "alive", "alive"…
## $ ICU_LOS <int> NA, NA, NA, NA, NA, 8, NA, NA, NA, NA, NA, NA, 2, NA…
## $ ward_LOS <int> NA, NA, NA, NA, NA, 22, NA, NA, NA, NA, NA, NA, NA, …
## $ ICU_admissions <dbl> 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0…
## $ OT_admissions <dbl> 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0…#View(survey_data1)
==================================================================================
time treat_control mean_duration ci_lower ci_upper min_duration max_duration
----------------------------------------------------------------------------------
1 baseline control Inf -Inf
2 baseline treatment Inf -Inf
3 endline control 9.533 8.679 10.387 2 15
4 endline treatment 3.750 3.147 4.353 2 8
----------------------------------------------------------------------------------
=========================================================================
treat_control mean_duration ci_lower ci_upper min_duration max_duration
-------------------------------------------------------------------------
1 control 9.533 8.679 10.387 2 15
2 treatment 3.750 3.147 4.353 2 8
-------------------------------------------------------------------------
| Variable | p_value | Test |
|---|---|---|
| DBPinitial_mean | 1.0000000 | Wilcoxon signed-rank test (within-group) |
| GCS_mean | 1.0000000 | Wilcoxon signed-rank test (within-group) |
| HRinitial_mean | 1.0000000 | Wilcoxon signed-rank test (within-group) |
| RRinitial_mean | 1.0000000 | Wilcoxon signed-rank test (within-group) |
| SBPinitial_mean | 1.0000000 | Wilcoxon signed-rank test (within-group) |
| SPO2initial_mean | 1.0000000 | Wilcoxon signed-rank test (within-group) |
| painscore_mean | 1.0000000 | Wilcoxon signed-rank test (within-group) |
| DBPinitial_mean | 0.6673311 | Mann-Whitney U test (between-group) |
| GCS_mean | 0.6673311 | Mann-Whitney U test (between-group) |
| HRinitial_mean | 0.6673311 | Mann-Whitney U test (between-group) |
| RRinitial_mean | 0.6673311 | Mann-Whitney U test (between-group) |
| SBPinitial_mean | 0.6673311 | Mann-Whitney U test (between-group) |
| SPO2initial_mean | 0.6673311 | Mann-Whitney U test (between-group) |
| painscore_mean | 0.6673311 | Mann-Whitney U test (between-group) |
| Variable | p_value | Test |
|---|---|---|
| DBPinitial_mean | 1.0000000 | Wilcoxon signed-rank test (within-group) |
| GCS_mean | 1.0000000 | Wilcoxon signed-rank test (within-group) |
| HRinitial_mean | 1.0000000 | Wilcoxon signed-rank test (within-group) |
| RRinitial_mean | 1.0000000 | Wilcoxon signed-rank test (within-group) |
| SBPinitial_mean | 1.0000000 | Wilcoxon signed-rank test (within-group) |
| SPO2initial_mean | 1.0000000 | Wilcoxon signed-rank test (within-group) |
| painscore_mean | 1.0000000 | Wilcoxon signed-rank test (within-group) |
| DBPinitial_mean | 0.6673311 | Mann-Whitney U test (between-group) |
| GCS_mean | 0.6673311 | Mann-Whitney U test (between-group) |
| HRinitial_mean | 0.6673311 | Mann-Whitney U test (between-group) |
| RRinitial_mean | 0.6673311 | Mann-Whitney U test (between-group) |
| SBPinitial_mean | 0.6673311 | Mann-Whitney U test (between-group) |
| SPO2initial_mean | 0.6673311 | Mann-Whitney U test (between-group) |
| painscore_mean | 0.6673311 | Mann-Whitney U test (between-group) |
0 1
58 162
Summary of Hospital Admission by Time and Treatment Control
==============================================================================
time treat_control total_hosp_admission total_obs percent_hosp_admission
------------------------------------------------------------------------------
1 baseline control 35 46 76.087
2 baseline treatment 42 57 73.684
3 endline control 49 64 76.562
4 endline treatment 36 53 67.925
------------------------------------------------------------------------------
Summary of Hosp Length of Stay by Time and Treatment Control
==================================================================================
time treat_control mean_duration ci_lower ci_upper min_duration max_duration
----------------------------------------------------------------------------------
1 baseline control 22.239 18.465 26.014 1 60
2 baseline treatment 23.702 19.111 28.292 1 65
3 endline control 24.453 21.582 27.324 1 62
4 endline treatment 18.547 16.437 20.657 3 32
----------------------------------------------------------------------------------
Hospital Mortality by Time and Treatment Control (with Percentages)
==================================================
time treat_control outcome n total Percentage
--------------------------------------------------
baseline control alive 41 46 89.1%
baseline control dead 5 46 10.9%
baseline treatment alive 44 57 77.2%
baseline treatment dead 13 57 22.8%
endline control alive 58 64 90.6%
endline control dead 6 64 9.4%
endline treatment alive 51 53 96.2%
endline treatment dead 2 53 3.8%
--------------------------------------------------
Hospital Mortality by Treatment Control (with Percentages)
=========================================
treat_control outcome n total Percentage
-----------------------------------------
control alive 99 110 90.0%
control dead 11 110 10.0%
treatment alive 95 110 86.4%
treatment dead 15 110 13.6%
-----------------------------------------| treat_control | p_value |
|---|---|
| control | 0.532 |
| treatment | 0.532 |
| time | treat_control | airwaystatus | breathingstatus | count | total | Percentage |
|---|---|---|---|---|---|---|
| baseline | control | obstructed | spont | 1 | 46 | 2.2% |
| baseline | control | obstructed | supplemental | 10 | 46 | 21.7% |
| baseline | control | patent | spont | 31 | 46 | 67.4% |
| baseline | control | patent | supplemental | 4 | 46 | 8.7% |
| baseline | treatment | obstructed | supplemental | 21 | 57 | 36.8% |
| baseline | treatment | patent | spont | 30 | 57 | 52.6% |
| baseline | treatment | patent | supplemental | 6 | 57 | 10.5% |
| endline | control | obstructed | supplemental | 15 | 64 | 23.4% |
| endline | control | patent | spont | 44 | 64 | 68.8% |
| endline | control | patent | supplemental | 5 | 64 | 7.8% |
| endline | treatment | obstructed | supplemental | 20 | 53 | 37.7% |
| endline | treatment | patent | spont | 28 | 53 | 52.8% |
| endline | treatment | patent | supplemental | 5 | 53 | 9.4% |
| Variable | Normality_Assumption | Test_Statistic | P_Value |
|---|---|---|---|
| ICU_admissions | FALSE | 0.49 | 0.0000 |
| ICU_LOS | TRUE | 0.94 | 0.1548 |
| hosp_admission | FALSE | 0.55 | 0.0000 |
| LOSdays | FALSE | 0.94 | 0.0000 |
| final_outcome | FALSE | 0.38 | 0.0000 |
| ward_LOS | TRUE | 0.96 | 0.3957 |
| Test | Variable | Test_Statistic | P_Value |
|---|---|---|---|
| Mann-Whitney U | ICU_admissions | 5940.000000 | 0.7377728 |
| Independent t-test | ICU_LOS | 4.926856 | 0.0000804 |
| Mann-Whitney U | hosp_admission | 6380.000000 | 0.3603883 |
| Mann-Whitney U | LOSdays | 6983.000000 | 0.0479941 |
| Mann-Whitney U | final_outcome | 6270.000000 | 0.4056435 |
| Test | Estimate | p_value | Df | Sum Sq | Mean Sq | F value |
|---|---|---|---|---|---|---|
| Independent Samples t-test | 9.533333 | 0.0000804 | NA | NA | NA | NA |
| Independent Samples t-test | 3.750000 | 0.0000804 | NA | NA | NA | NA |
| treat_control | NA | 0.0003093 | 1 | 174.5058 | 174.505797 | 18.58013 |
| Residuals | NA | NA | 21 | 197.2333 | 9.392064 | NA |
Call:
randomForest(formula = ICU_LOS ~ ., data = subset_data, importance = TRUE, na.action = na.omit)
Type of random forest: regression
Number of trees: 500
No. of variables tried at each split: 2
Mean of squared residuals: 13.71438
% Var explained: 15.15
Call:
glm(formula = treat_control ~ ., family = "binomial", data = subset_data)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.657e+01 1.334e+05 0.000 1.000
timebaseline 2.472e-12 5.540e+04 0.000 1.000
genderfemale -3.569e-11 5.821e+04 0.000 1.000
occupstudent 2.671e-10 1.591e+05 0.000 1.000
occupbusiness 2.420e-10 1.287e+05 0.000 1.000
occuphousewife 2.533e-10 1.526e+05 0.000 1.000
occupcasual labourer 2.211e-10 1.498e+05 0.000 1.000
occupbar attendant 2.314e-10 2.837e+05 0.000 1.000
occuppubic servant 2.157e-10 3.874e+05 0.000 1.000
occupfarmer 2.926e-10 1.368e+05 0.000 1.000
occuptrader 3.258e-10 3.878e+05 0.000 1.000
occupteacher 2.996e-10 1.632e+05 0.000 1.000
occupshop attendant 5.626e-10 2.875e+05 0.000 1.000
occuptechnician 1.232e-11 3.144e+05 0.000 1.000
occupadministrator 3.605e-10 3.907e+05 0.000 1.000
occuptailor 2.461e-10 2.883e+05 0.000 1.000
occupsalonist 3.615e-10 3.889e+05 0.000 1.000
occupmason 2.322e-10 2.857e+05 0.000 1.000
occupaccountant 2.498e-10 2.855e+05 0.000 1.000
occupnone 2.317e-10 1.893e+05 0.000 1.000
occupbodarider 1.878e-10 2.189e+05 0.000 1.000
occupsaleslady 2.535e-10 3.945e+05 0.000 1.000
occupnurse 2.180e-10 3.886e+05 0.000 1.000
occupbarber 1.935e-10 3.803e+05 0.000 1.000
occupdriver 2.356e-10 2.441e+05 0.000 1.000
occuppreacher 1.852e-10 3.783e+05 0.000 1.000
occupfishing 1.935e-10 3.803e+05 0.000 1.000
occupvendor -4.413e-06 3.776e+05 0.000 1.000
occuprtd public servant -9.330e-10 3.788e+05 0.000 1.000
occupvet officer -1.984e-07 3.866e+05 0.000 1.000
occupmechanic -4.416e-06 3.788e+05 0.000 1.000
occuppastor 4.424e-06 3.931e+05 0.000 1.000
occupoffic assistant 4.619e-06 3.832e+05 0.000 1.000
occupcobbler 2.134e-07 3.776e+05 0.000 1.000
occupclerk 4.419e-06 3.967e+05 0.000 1.000
mstatussingle 1.014e-11 8.515e+04 0.000 1.000
mstatuswidowed -4.410e-06 2.625e+05 0.000 1.000
edulevelcollege -1.322e-10 9.308e+04 0.000 1.000
edulevelprimary 8.461e-12 6.498e+04 0.000 1.000
eduleveluniversity 2.815e-11 2.457e+05 0.000 1.000
countymigori -3.459e-11 7.521e+04 0.000 1.000
countykakamega 5.313e+01 7.333e+04 0.001 0.999
countyvihiga 5.313e+01 7.811e+04 0.001 0.999
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 3.0498e+02 on 219 degrees of freedom
Residual deviance: 1.2763e-09 on 177 degrees of freedom
AIC: 86
Number of Fisher Scoring iterations: 25# A tibble: 43 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) -2.66e+ 1 133398. -1.99e- 4 1.00
2 timebaseline 2.47e-12 55402. 4.46e-17 1
3 genderfemale -3.57e-11 58211. -6.13e-16 1
4 occupstudent 2.67e-10 159085. 1.68e-15 1.00
5 occupbusiness 2.42e-10 128652. 1.88e-15 1.00
6 occuphousewife 2.53e-10 152605. 1.66e-15 1.00
7 occupcasual labourer 2.21e-10 149816. 1.48e-15 1.00
8 occupbar attendant 2.31e-10 283717. 8.16e-16 1.00
9 occuppubic servant 2.16e-10 387356. 5.57e-16 1
10 occupfarmer 2.93e-10 136768. 2.14e-15 1.00
# ℹ 33 more rows
Call:
glm(formula = ICU_admissions ~ age + gender + MOI + triagecategory +
Edduration_mins, family = binomial, data = survey_data1)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.637e+01 2.790e+03 -0.006 0.9953
age -1.944e-02 1.259e-02 -1.543 0.1227
gendermale 2.378e-01 3.795e-01 0.627 0.5308
MOIburns 7.664e-01 1.025e+00 0.748 0.4547
MOIcrushing -1.513e+01 3.956e+03 -0.004 0.9969
MOIcut -1.504e+01 1.967e+03 -0.008 0.9939
MOIfall 5.429e-01 9.144e-01 0.594 0.5527
MOIgunshot 9.829e-01 1.371e+00 0.717 0.4734
MOIpenetrating inj 2.016e+01 3.956e+03 0.005 0.9959
MOIrti 8.753e-01 8.063e-01 1.086 0.2777
MOIrti, motorbike -1.596e+01 2.797e+03 -0.006 0.9954
MOIrti,fall -1.450e+01 3.956e+03 -0.004 0.9971
MOIstab 9.986e-01 8.956e-01 1.115 0.2649
MOIstab, assault -1.567e+01 3.956e+03 -0.004 0.9968
MOIstab/ cut -1.619e+01 3.956e+03 -0.004 0.9967
triagecategoryimmediate 1.565e+01 2.790e+03 0.006 0.9955
triagecategoryurgent 1.558e+01 2.790e+03 0.006 0.9955
Edduration_mins -1.460e-02 7.746e-03 -1.885 0.0594 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 220.18 on 219 degrees of freedom
Residual deviance: 201.00 on 202 degrees of freedom
AIC: 237
Number of Fisher Scoring iterations: 16
Call:
glm(formula = ICU_LOS ~ age + gender + MOI + triagecategory +
Edduration_mins, family = poisson, data = survey_data1)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.982985 0.529857 1.855 0.0636 .
age 0.007025 0.007037 0.998 0.3181
gendermale -0.048132 0.228867 -0.210 0.8334
MOIburns 0.152881 0.521656 0.293 0.7695
MOIfall 0.297835 0.356988 0.834 0.4041
MOIgunshot -0.741781 0.819136 -0.906 0.3652
MOIpenetrating inj 0.376254 0.405687 0.927 0.3537
MOIrti 0.253383 0.323121 0.784 0.4329
MOIstab -0.150723 0.425471 -0.354 0.7232
triagecategoryurgent 0.344267 0.193943 1.775 0.0759 .
Edduration_mins 0.007403 0.004400 1.683 0.0925 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 53.835 on 22 degrees of freedom
Residual deviance: 29.789 on 12 degrees of freedom
(197 observations deleted due to missingness)
AIC: 136.85
Number of Fisher Scoring iterations: 5
Call:
glm(formula = hosp_admission ~ age + gender + MOI + triagecategory +
Edduration_mins, family = binomial, data = survey_data1)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.78481 1.79113 -0.996 0.3190
age -0.01845 0.01141 -1.617 0.1059
gendermale 0.11801 0.35376 0.334 0.7387
MOIburns 17.10684 1028.48362 0.017 0.9867
MOIcrushing 16.18263 3956.18039 0.004 0.9967
MOIcut 0.89946 1.29612 0.694 0.4877
MOIfall 1.39684 0.67076 2.082 0.0373 *
MOIgunshot 0.30401 1.07141 0.284 0.7766
MOIpenetrating inj 17.29640 3956.18036 0.004 0.9965
MOIrti 0.98619 0.53956 1.828 0.0676 .
MOIrti, motorbike 17.90193 2791.11385 0.006 0.9949
MOIrti,fall 16.19949 3956.18044 0.004 0.9967
MOIstab 0.32757 0.64658 0.507 0.6124
MOIstab, assault 17.23183 3956.18037 0.004 0.9965
MOIstab/ cut -17.08353 3956.18037 -0.004 0.9966
triagecategoryimmediate 1.46665 1.56881 0.935 0.3498
triagecategoryurgent 2.16562 1.56424 1.384 0.1662
Edduration_mins 0.01290 0.00737 1.751 0.0799 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 253.80 on 219 degrees of freedom
Residual deviance: 219.42 on 202 degrees of freedom
AIC: 255.42
Number of Fisher Scoring iterations: 16Call:
coxph(formula = Surv(LOSdays, final_outcome) ~ age + gender +
MOI + triagecategory + Edduration_mins, data = survey_data1)
n= 220, number of events= 194
coef exp(coef) se(coef) z Pr(>|z|)
age -0.001402 0.998599 0.005713 -0.245 0.806086
gendermale -0.205337 0.814372 0.162341 -1.265 0.205924
MOIburns -0.190415 0.826616 0.370534 -0.514 0.607325
MOIcrushing 0.463716 1.589972 1.050479 0.441 0.658899
MOIcut 0.505075 1.657109 0.561529 0.899 0.368406
MOIfall -0.253544 0.776045 0.300521 -0.844 0.398847
MOIgunshot 1.138839 3.123140 0.641956 1.774 0.076061 .
MOIpenetrating inj -0.799076 0.449744 1.033650 -0.773 0.439485
MOIrti 0.124619 1.132716 0.262005 0.476 0.634334
MOIrti, motorbike 2.829903 16.943822 0.800128 3.537 0.000405 ***
MOIrti,fall 1.607653 4.991082 1.087901 1.478 0.139473
MOIstab 0.545869 1.726108 0.329622 1.656 0.097713 .
MOIstab, assault -0.175994 0.838623 1.044521 -0.168 0.866196
MOIstab/ cut 2.333049 10.309325 1.060419 2.200 0.027798 *
triagecategoryimmediate 0.066532 1.068795 0.759963 0.088 0.930237
triagecategoryurgent 0.340965 1.406303 0.747273 0.456 0.648190
Edduration_mins -0.006415 0.993606 0.002898 -2.214 0.026860 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
age 0.9986 1.00140 0.98748 1.0098
gendermale 0.8144 1.22794 0.59243 1.1195
MOIburns 0.8266 1.20975 0.39986 1.7088
MOIcrushing 1.5900 0.62894 0.20287 12.4612
MOIcut 1.6571 0.60346 0.55129 4.9811
MOIfall 0.7760 1.28858 0.43061 1.3986
MOIgunshot 3.1231 0.32019 0.88748 10.9907
MOIpenetrating inj 0.4497 2.22349 0.05931 3.4104
MOIrti 1.1327 0.88283 0.67780 1.8930
MOIrti, motorbike 16.9438 0.05902 3.53135 81.2983
MOIrti,fall 4.9911 0.20036 0.59180 42.0938
MOIstab 1.7261 0.57934 0.90468 3.2934
MOIstab, assault 0.8386 1.19243 0.10826 6.4963
MOIstab/ cut 10.3093 0.09700 1.29003 82.3873
triagecategoryimmediate 1.0688 0.93563 0.24100 4.7400
triagecategoryurgent 1.4063 0.71108 0.32509 6.0836
Edduration_mins 0.9936 1.00644 0.98798 0.9993
Concordance= 0.637 (se = 0.02 )
Likelihood ratio test= 29.24 on 17 df, p=0.03
Wald test = 33.76 on 17 df, p=0.009
Score (logrank) test = 43.54 on 17 df, p=4e-04
Call:
glm(formula = final_outcome ~ age + gender + MOI + triagecategory +
Durationmins, family = binomial, data = survey_data1)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 37.16721 4719.48517 0.008 0.9937
age 0.01551 0.01582 0.981 0.3268
gendermale 0.10787 0.47641 0.226 0.8209
MOIburns -15.86982 1505.27173 -0.011 0.9916
MOIcrushing 0.26637 6694.07621 0.000 1.0000
MOIcut -0.42861 3579.05387 0.000 0.9999
MOIfall -15.46169 1505.27171 -0.010 0.9918
MOIgunshot -18.68045 1505.27167 -0.012 0.9901
MOIpenetrating inj -0.02906 6694.07620 0.000 1.0000
MOIrti -17.03081 1505.27140 -0.011 0.9910
MOIrti, motorbike -0.24325 4849.66071 0.000 1.0000
MOIrti,fall 0.33895 6694.07623 0.000 1.0000
MOIstab -17.59175 1505.27146 -0.012 0.9907
MOIstab, assault -0.06641 6694.07623 0.000 1.0000
MOIstab/ cut -0.87010 6694.07623 0.000 0.9999
triagecategoryimmediate -18.14212 4472.99648 -0.004 0.9968
triagecategoryurgent -18.10056 4472.99643 -0.004 0.9968
Durationmins -0.04054 0.02098 -1.932 0.0533 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 159.85 on 219 degrees of freedom
Residual deviance: 138.31 on 202 degrees of freedom
AIC: 174.31
Number of Fisher Scoring iterations: 17# A tibble: 31 × 4
gender MOI triagecategory ICU_admissions
<chr> <chr> <chr> <dbl>
1 female assault immediate 0
2 female assault urgent 0.125
3 female burns urgent 0.333
4 female cut urgent 0
5 female fall immediate 0.333
6 female fall urgent 0.0769
7 female gunshot urgent 0
8 female penetrating inj urgent 1
9 female rti delayed 0
10 female rti immediate 0.286
# ℹ 21 more rows
Pearson's Chi-squared test with Yates' continuity correction
data: table(survey_data1$gender, survey_data1$ICU_admissions)
X-squared = 0.33851, df = 1, p-value = 0.5607 Df Sum Sq Mean Sq F value Pr(>F)
ICU_admissions 1 5146 5146 3.718 0.0551 .
Residuals 218 301761 1384
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Call:
lm(formula = ICU_LOS ~ age + gender + MOI + triagecategory +
Edduration_mins, data = survey_data1)
Residuals:
Min 1Q Median 3Q Max
-4.601 -2.167 0.000 1.537 7.386
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.71544 6.13548 -0.117 0.909
age 0.06842 0.08272 0.827 0.424
gendermale 0.06057 2.68300 0.023 0.982
MOIburns 1.00730 5.76987 0.175 0.864
MOIfall 2.86106 4.08864 0.700 0.497
MOIgunshot -1.76774 6.05557 -0.292 0.775
MOIpenetrating inj 3.86964 5.14916 0.752 0.467
MOIrti 1.88662 3.67280 0.514 0.617
MOIstab -0.40283 4.53041 -0.089 0.931
triagecategoryurgent 2.45636 2.11075 1.164 0.267
Edduration_mins 0.05176 0.04797 1.079 0.302
Residual standard error: 4.112 on 12 degrees of freedom
(197 observations deleted due to missingness)
Multiple R-squared: 0.4541, Adjusted R-squared: -0.0008499
F-statistic: 0.9981 on 10 and 12 DF, p-value: 0.4938
Call:
glm(formula = ICU_admissions ~ age + gender + MOI + triagecategory +
Edduration_mins, family = binomial, data = survey_data1)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.637e+01 2.790e+03 -0.006 0.9953
age -1.944e-02 1.259e-02 -1.543 0.1227
gendermale 2.378e-01 3.795e-01 0.627 0.5308
MOIburns 7.664e-01 1.025e+00 0.748 0.4547
MOIcrushing -1.513e+01 3.956e+03 -0.004 0.9969
MOIcut -1.504e+01 1.967e+03 -0.008 0.9939
MOIfall 5.429e-01 9.144e-01 0.594 0.5527
MOIgunshot 9.829e-01 1.371e+00 0.717 0.4734
MOIpenetrating inj 2.016e+01 3.956e+03 0.005 0.9959
MOIrti 8.753e-01 8.063e-01 1.086 0.2777
MOIrti, motorbike -1.596e+01 2.797e+03 -0.006 0.9954
MOIrti,fall -1.450e+01 3.956e+03 -0.004 0.9971
MOIstab 9.986e-01 8.956e-01 1.115 0.2649
MOIstab, assault -1.567e+01 3.956e+03 -0.004 0.9968
MOIstab/ cut -1.619e+01 3.956e+03 -0.004 0.9967
triagecategoryimmediate 1.565e+01 2.790e+03 0.006 0.9955
triagecategoryurgent 1.558e+01 2.790e+03 0.006 0.9955
Edduration_mins -1.460e-02 7.746e-03 -1.885 0.0594 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 220.18 on 219 degrees of freedom
Residual deviance: 201.00 on 202 degrees of freedom
AIC: 237
Number of Fisher Scoring iterations: 16 (Intercept) age gendermale
7.798393e-08 9.807503e-01 1.268500e+00
MOIburns MOIcrushing MOIcut
2.152075e+00 2.689826e-07 2.939642e-07
MOIfall MOIgunshot MOIpenetrating inj
1.720948e+00 2.672143e+00 5.708529e+08
MOIrti MOIrti, motorbike MOIrti,fall
2.399589e+00 1.171496e-07 5.021493e-07
MOIstab MOIstab, assault MOIstab/ cut
2.714438e+00 1.561369e-07 9.272204e-08
triagecategoryimmediate triagecategoryurgent Edduration_mins
6.238705e+06 5.834541e+06 9.855040e-01 # Extract odds ratios and p-values from the model summary
summary_model <- summary(logit_model)
odds_ratios <- exp(coef(logit_model))
p_values <- coef(summary_model)[, "Pr(>|z|)"]
# Create a data frame with odds ratios and p-values
results_df <- data.frame(
Variable = names(odds_ratios),
Odds_Ratio = odds_ratios,
P_Value = p_values
)
# Load knitr and gt for formatting
library(knitr)
library(gt)
# Option 1: Using kable from knitr to print a simple table
kable(results_df, digits = 3, caption = "Odds Ratios and P-values from Logistic Regression Model")| Variable | Odds_Ratio | P_Value | |
|---|---|---|---|
| (Intercept) | (Intercept) | 0.000000e+00 | 0.995 |
| age | age | 9.810000e-01 | 0.123 |
| gendermale | gendermale | 1.268000e+00 | 0.531 |
| MOIburns | MOIburns | 2.152000e+00 | 0.455 |
| MOIcrushing | MOIcrushing | 0.000000e+00 | 0.997 |
| MOIcut | MOIcut | 0.000000e+00 | 0.994 |
| MOIfall | MOIfall | 1.721000e+00 | 0.553 |
| MOIgunshot | MOIgunshot | 2.672000e+00 | 0.473 |
| MOIpenetrating inj | MOIpenetrating inj | 5.708529e+08 | 0.996 |
| MOIrti | MOIrti | 2.400000e+00 | 0.278 |
| MOIrti, motorbike | MOIrti, motorbike | 0.000000e+00 | 0.995 |
| MOIrti,fall | MOIrti,fall | 0.000000e+00 | 0.997 |
| MOIstab | MOIstab | 2.714000e+00 | 0.265 |
| MOIstab, assault | MOIstab, assault | 0.000000e+00 | 0.997 |
| MOIstab/ cut | MOIstab/ cut | 0.000000e+00 | 0.997 |
| triagecategoryimmediate | triagecategoryimmediate | 6.238705e+06 | 0.996 |
| triagecategoryurgent | triagecategoryurgent | 5.834541e+06 | 0.996 |
| Edduration_mins | Edduration_mins | 9.860000e-01 | 0.059 |
# Option 2: Using gt to create a more styled table (recommended for publication)
results_df %>%
gt() %>%
tab_header(
title = "Odds Ratios and P-values from Logistic Regression Model"
) %>%
fmt_number(
columns = vars(Odds_Ratio),
decimals = 2
) %>%
fmt_number(
columns = vars(P_Value),
decimals = 3
) %>%
cols_label(
Variable = "Variable",
Odds_Ratio = "Odds Ratio",
P_Value = "P-Value"
)## Warning: Since gt v0.3.0, `columns = vars(...)` has been deprecated.
## • Please use `columns = c(...)` instead.
## Since gt v0.3.0, `columns = vars(...)` has been deprecated.
## • Please use `columns = c(...)` instead.| Odds Ratios and P-values from Logistic Regression Model | ||
| Variable | Odds Ratio | P-Value |
|---|---|---|
| (Intercept) | 0.00 | 0.995 |
| age | 0.98 | 0.123 |
| gendermale | 1.27 | 0.531 |
| MOIburns | 2.15 | 0.455 |
| MOIcrushing | 0.00 | 0.997 |
| MOIcut | 0.00 | 0.994 |
| MOIfall | 1.72 | 0.553 |
| MOIgunshot | 2.67 | 0.473 |
| MOIpenetrating inj | 570,852,931.75 | 0.996 |
| MOIrti | 2.40 | 0.278 |
| MOIrti, motorbike | 0.00 | 0.995 |
| MOIrti,fall | 0.00 | 0.997 |
| MOIstab | 2.71 | 0.265 |
| MOIstab, assault | 0.00 | 0.997 |
| MOIstab/ cut | 0.00 | 0.997 |
| triagecategoryimmediate | 6,238,705.26 | 0.996 |
| triagecategoryurgent | 5,834,541.17 | 0.996 |
| Edduration_mins | 0.99 | 0.059 |
# Convert categorical variables to factors
survey_data <- survey_data1 %>%
mutate(
treat_control = factor(treat_control, levels = c("control", "intervention")),
time = factor(time, levels = c("baseline", "endline")),
county = factor(county),
gender = factor(gender),
outcome = factor(outcome, levels = c("alive", "dead")),
MOI = factor(MOI),
anatlocation = factor(anatlocation),
triagecategory = factor(triagecategory),
airwaystatus = factor(airwaystatus),
airwayinterven = factor(airwayinterven),
spinestabilized = factor(spinestabilized),
breathingstatus = factor(breathingstatus),
breathinginterven = factor(breathinginterven),
bleedingcontrol = factor(bleedingcontrol)
)# Handle missing data (example: remove rows with any NA)
survey_data <- survey_data %>% drop_na()
# 1. Descriptive Statistics
# --------------------------
# Load necessary libraries
library(dplyr)
library(stargazer)
library(dplyr)
library(stringr)
# Change all values in the 'county' variable to lowercase
survey_data <- survey_data %>%
mutate(county = str_to_lower(county))
# Create outcome variables
survey_data1$ICU_admissions <- ifelse(survey_data1$Admit == "icu", 1, 0)
survey_data1$hosp_admission <- ifelse(survey_data1$Admit == "ward", 1, 0)
survey_data1$OT_admission <- ifelse(survey_data1$Admit == "OT", 1, 0)
# Calculate descriptive statistics by treatment group, time (baseline/endline), and county
library(dplyr)
library(tidyr)
# Define the levels for each variable
county_levels <- c("Siaya", "Migori", "Kakamega", "Vihiga") # Replace with actual counties
time_levels <- c("baseline", "endline") # Replace with actual timepoints
treat_control_levels <- c("treatment", "control") # Replace with actual treatment/control groups
descriptive_stats <- survey_data %>%
# Ensure all combinations of 'treat_control', 'time', and 'county' are used
complete(treat_control = survey_data1$treat_control,
time = survey_data1$time,
county = survey_data1$county,
fill = list(
Age_Mean = NA,
Age_SD = NA,
Gender_Male_Percent = NA,
Gender_Female_Percent = NA,
Hosp_admission_Rate = NA,
ICU_admission_Rate = NA,
OT_admission_Rate = NA,
Mortality_Rate = NA
)) %>%
group_by(treat_control, time, county) %>%
summarize(
Age_Mean = mean(age, na.rm = TRUE),
Age_SD = sd(age, na.rm = TRUE),
Gender_Male_Percent = mean(gender == "male", na.rm = TRUE) * 100,
Gender_Female_Percent = mean(gender == "female", na.rm = TRUE) * 100,
Hosp_admission_Rate = mean(survey_data1$hosp_admission, na.rm = TRUE) * 100,
ICU_admission_Rate = mean(survey_data1$ICU_admissions, na.rm = TRUE) * 100,
OT_admission_Rate = mean(survey_data1$OT_admission, na.rm = TRUE) * 100,
Mortality_Rate = mean(outcome == "dead", na.rm = TRUE) * 100,
.groups = "drop" # This removes any unnecessary grouping after summarization
)%>%
# Remove rows with NaN values
filter(
!is.nan(Age_Mean) & !is.nan(Age_SD) &
!is.nan(Gender_Male_Percent) & !is.nan(Gender_Female_Percent) &
!is.nan(Hosp_admission_Rate) & !is.nan(ICU_admission_Rate) &
!is.nan(OT_admission_Rate) & !is.nan(Mortality_Rate)
)
# View the results
print(descriptive_stats)## # A tibble: 2 × 11
## treat_control time county Age_Mean Age_SD Gender_Male_Percent
## <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 control endline migori 35.5 7.78 100
## 2 control endline siaya 33 7.07 50
## # ℹ 5 more variables: Gender_Female_Percent <dbl>, Hosp_admission_Rate <dbl>,
## # ICU_admission_Rate <dbl>, OT_admission_Rate <dbl>, Mortality_Rate <dbl># Display Descriptive Statistics
descriptive_stats %>%
gt() %>%
tab_header(title = "Descriptive Statistics by Treatment Group")| Descriptive Statistics by Treatment Group | ||||||||||
| treat_control | time | county | Age_Mean | Age_SD | Gender_Male_Percent | Gender_Female_Percent | Hosp_admission_Rate | ICU_admission_Rate | OT_admission_Rate | Mortality_Rate |
|---|---|---|---|---|---|---|---|---|---|---|
| control | endline | migori | 35.5 | 7.778175 | 100 | 0 | 52.27273 | 20 | 26.36364 | 0 |
| control | endline | siaya | 33.0 | 7.071068 | 50 | 50 | 52.27273 | 20 | 26.36364 | 0 |
# Load required libraries
library(tidyverse)
library(stargazer)
library(ggplot2)
library(skimr)
library(broom)
# Convert categorical variables to factors
survey_data1 <- survey_data1 %>%
mutate(
treat_control = as.factor(treat_control),
gender = as.factor(gender),
mstatus = as.factor(mstatus),
edulevel = as.factor(edulevel),
occup = as.factor(occup),
MOI = as.factor(MOI),
anatlocation = as.factor(anatlocation),
triagecategory = as.factor(triagecategory),
arrivalmode = as.factor(arrivalmode),
refby = as.factor(refby),
airwaystatus = as.factor(airwaystatus),
breathingstatus = as.factor(breathingstatus),
Admit = as.factor(Admit),
outcome = as.factor(outcome)
)# Generate summary statistics for LOS and key variables
stargazer(survey_data1[, c("LOSdays", "Edduration_mins", "HRinitial")],
type = "text",
title = "Summary Statistics for Key Variables",
digits = 2,
summary.stat = c("mean", "sd", "min", "max", "n"))##
## Summary Statistics for Key Variables
## ==========================================
## Statistic Mean St. Dev. Min Max N
## ------------------------------------------
## LOSdays 22.37 12.92 1 65 220
## Edduration_mins 55.67 37.44 10 360 220
## HRinitial 97.77 15.11 63 150 220
## ------------------------------------------skim(survey_data1[, c("LOSdays", "Edduration_mins", "HRinitial", "SBPinitial", "DBPinitial")])| Name | …[] |
| Number of rows | 220 |
| Number of columns | 5 |
| _______________________ | |
| Column type frequency: | |
| numeric | 5 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| LOSdays | 0 | 1 | 22.37 | 12.92 | 1 | 15 | 21.5 | 29.25 | 65 | ▃▇▅▁▁ |
| Edduration_mins | 0 | 1 | 55.67 | 37.44 | 10 | 40 | 48.0 | 60.00 | 360 | ▇▁▁▁▁ |
| HRinitial | 0 | 1 | 97.77 | 15.11 | 63 | 90 | 98.0 | 108.50 | 150 | ▂▆▇▂▁ |
| SBPinitial | 0 | 1 | 107.27 | 20.83 | 45 | 92 | 105.0 | 124.00 | 175 | ▁▆▇▅▁ |
| DBPinitial | 0 | 1 | 72.91 | 11.97 | 30 | 68 | 74.0 | 80.00 | 100 | ▁▁▇▇▂ |
ggplot(survey_data1, aes(x = LOSdays)) +
geom_histogram(fill = "steelblue", bins = 30, alpha = 0.7) +
theme_minimal() +
labs(title = "Distribution of Length of Stay (LOS)",
x = "Length of Stay (Days)",
y = "Count")
ggplot(survey_data1, aes(x = treat_control, y = LOSdays, fill = treat_control)) +
geom_boxplot(alpha = 0.7) +
theme_minimal() +
labs(title = "Comparison of LOS Between Treatment and Control Groups",
x = "Group",
y = "Length of Stay (Days)")
library(ggplot2)
library(dplyr)
library(ggplot2)
library(dplyr)
# Check ISS range
summary(survey_data1$ISS)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 15.00 15.00 15.00 15.36 16.00 16.00# Create a more granular ISS grouping focusing on 15 and 16
survey_data1 <- survey_data1 %>%
mutate(ISS_category = case_when(
ISS < 15 ~ "<15",
ISS == 15 ~ "15",
ISS == 16 ~ "16",
ISS > 16 ~ ">16"
))
# Convert to factor to maintain order
survey_data1$ISS_category <- factor(survey_data1$ISS_category,
levels = c("<15", "15", "16", ">16"))
# Boxplot of LOS grouped by refined ISS categories
ggplot(survey_data1, aes(x = ISS_category, y = LOSdays, fill = treat_control)) +
geom_boxplot(alpha = 0.7, outlier.shape = NA) +
geom_jitter(position = position_jitter(0.2), alpha = 0.3) +
theme_minimal() +
labs(title = "LOS Across ISS Categories",
x = "ISS Category",
y = "Length of Stay (Days)",
fill = "Group") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
# Option 2: If ISS is more variable, improve scatter plot
ggplot(survey_data1, aes(x = ISS, y = LOSdays, color = treat_control)) +
geom_jitter(alpha = 0.5, width = 0.2) +
geom_smooth(method = "loess", se = TRUE) + # Use loess for non-linear trends
theme_minimal() +
labs(title = "LOS vs ISS (Non-Linear Fit)",
x = "Injury Severity Score (ISS)",
y = "Length of Stay (Days)",
color = "Group")## `geom_smooth()` using formula = 'y ~ x'
# Run ANOVA separately for Treatment and Control groups
anova_treat <- aov(LOSdays ~ gender + mstatus + ISS + MOI + triagecategory, data = subset(survey_data1, treat_control == "treatment"))
anova_control <- aov(LOSdays ~ gender + mstatus + ISS + MOI + triagecategory, data = subset(survey_data1, treat_control == "control"))
# Display ANOVA summaries
summary(anova_treat)## Df Sum Sq Mean Sq F value Pr(>F)
## gender 1 69 68.7 0.390 0.5338
## mstatus 2 16 8.2 0.046 0.9548
## ISS 1 2 1.6 0.009 0.9236
## MOI 7 2979 425.6 2.414 0.0255 *
## triagecategory 2 563 281.6 1.597 0.2078
## Residuals 96 16924 176.3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(anova_control)## Df Sum Sq Mean Sq F value Pr(>F)
## gender 1 3 2.9 0.020 0.8866
## mstatus 1 29 28.7 0.205 0.6520
## ISS 1 638 638.4 4.555 0.0354 *
## MOI 11 1872 170.2 1.215 0.2885
## triagecategory 2 159 79.4 0.567 0.5693
## Residuals 93 13034 140.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1install.packages("epitools")## Installing package into 'C:/Users/PC/AppData/Local/R/win-library/4.4'
## (as 'lib' is unspecified)## package 'epitools' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\PC\AppData\Local\Temp\Rtmp2pV6Zk\downloaded_packageslibrary(epitools)##
## Attaching package: 'epitools'## The following object is masked from 'package:survival':
##
## ratetable# Convert categorical variables to factors
survey_data1 <- survey_data1 %>%
mutate(gender = as.factor(gender),
mstatus = as.factor(mstatus),
MOI = as.factor(MOI),
triagecategory = as.factor(triagecategory),
treat_control = as.factor(treat_control))
# Mantel-Haenszel test to compare ISS categories across groups
mh_test <- mantelhaen.test(
x = survey_data1$triagecategory,
y = survey_data1$treat_control,
z = survey_data1$ISS
)
print(mh_test)##
## Cochran-Mantel-Haenszel test
##
## data: survey_data1$triagecategory and survey_data1$treat_control and survey_data1$ISS
## Cochran-Mantel-Haenszel M^2 = 2.5275, df = 2, p-value = 0.2826ggplot(survey_data1, aes(x = reorder(triagecategory, LOSdays, FUN = median), y = LOSdays, fill = treat_control)) +
geom_boxplot(alpha = 0.7, outlier.shape = NA) +
geom_jitter(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.8), alpha = 0.3, size = 1) +
theme_minimal() +
labs(title = "LOS Across Triage Categories by Group",
x = "Triage Category",
y = "Length of Stay (Days)",
fill = "Group") +
facet_wrap(~treat_control, ncol = 1) + # Stack the plots for better readability
theme(axis.text.x = element_text(angle = 30, hjust = 1), legend.position = "top") 
# Load necessary libraries for HCP analysis
library(dplyr)
library(ggplot2)
library(tidyr)
library(caret)
library(stargazer)# Load the dataset
hcp_data <- read.csv("C:/Users/PC/Downloads/hcp_data.csv")
head(hcp_data)## interv_ctrl code county age gender profession speciality experience
## 1 intervention 1 kakamega 33 female rn none 6
## 2 intervention 2 kakamega 25 female rco none 1
## 3 intervention 3 kakamega 26 male bsn none 4
## 4 intervention 4 kakamega 26 female rn none 3
## 5 intervention 5 kakamega 19 male rn midwife 1
## 6 intervention 6 kakamega 23 male rn midwife 1
## deployment cont_education lastattended SOP c.collar.applied
## 1 6 bls 4 yes correctly performed
## 2 1 none NA yes correctly performed
## 3 2 bls 1 yes correctly performed
## 4 1 bls 1 yes correctly performed
## 5 1 none NA yes correctly performed
## 6 1 none NA no correctly performed
## airway.insert suction ETT
## 1 correctly performed correctly performed not performed
## 2 correctly performed correctly performed correctly performed
## 3 correctly performed correctly performed correctly performed
## 4 correctly performed correctly performed not performed
## 5 not performed correctly performed not performed
## 6 correctly performed correctly performed incorrectly performed
## oxygensupp ventilator chesttube
## 1 correctly performed incorrectly performed incorrectly performed
## 2 correctly performed incorrectly performed correctly performed
## 3 correctly performed incorrectly performed not performed
## 4 correctly performed incorrectly performed not performed
## 5 correctly performed incorrectly performed incorrectly performed
## 6 correctly performed incorrectly performed incorrectly performed
## directpressure tourniquet splint
## 1 correctly performed incorrectly performed correctly performed
## 2 correctly performed incorrectly performed correctly performed
## 3 correctly performed incorrectly performed correctly performed
## 4 correctly performed incorrectly performed correctly performed
## 5 correctly performed incorrectly performed correctly performed
## 6 correctly performed incorrectly performed correctly performed
## pelvic.binder IV.Access IVF
## 1 incorrectly performed correctly performed correctly performed
## 2 incorrectly performed correctly performed correctly performed
## 3 not performed correctly performed correctly performed
## 4 incorrectly performed correctly performed correctly performed
## 5 incorrectly performed correctly performed correctly performed
## 6 incorrectly performed correctly performed correctly performed
## cardiac.monitoring resp.monitor neurological
## 1 correctly performed correctly performed correctly performed
## 2 correctly performed correctly performed correctly performed
## 3 correctly performed correctly performed correctly performed
## 4 correctly performed correctly performed correctly performed
## 5 correctly performed correctly performed correctly performed
## 6 correctly performed correctly performed correctly performed# Get the variable names (column names)
variable_names <- names(hcp_data)
# Print the variable names
print(variable_names)## [1] "interv_ctrl" "code" "county"
## [4] "age" "gender" "profession"
## [7] "speciality" "experience" "deployment"
## [10] "cont_education" "lastattended" "SOP"
## [13] "c.collar.applied" "airway.insert" "suction"
## [16] "ETT" "oxygensupp" "ventilator"
## [19] "chesttube" "directpressure" "tourniquet"
## [22] "splint" "pelvic.binder" "IV.Access"
## [25] "IVF" "cardiac.monitoring" "resp.monitor"
## [28] "neurological"# Load necessary libraries
library(dplyr)
library(ggplot2)
library(tidyr)
# Summary statistics for intervention and control groups
summary_intervention <- hcp_data %>%
filter(interv_ctrl == "intervention") %>%
summary()
summary_control <- hcp_data %>%
filter(interv_ctrl == "control") %>%
summary()
print("Summary for Intervention Group:")## [1] "Summary for Intervention Group:"print(summary_intervention)## interv_ctrl code county age
## Length:30 Length:30 Length:30 Min. :19.00
## Class :character Class :character Class :character 1st Qu.:25.25
## Mode :character Mode :character Mode :character Median :29.50
## Mean :30.43
## 3rd Qu.:35.00
## Max. :43.00
##
## gender profession speciality experience
## Length:30 Length:30 Length:30 Min. : 1.000
## Class :character Class :character Class :character 1st Qu.: 1.250
## Mode :character Mode :character Mode :character Median : 3.500
## Mean : 4.467
## 3rd Qu.: 6.000
## Max. :15.000
##
## deployment cont_education lastattended SOP
## Min. : 1.000 Length:30 Min. :1.000 Length:30
## 1st Qu.: 1.000 Class :character 1st Qu.:1.000 Class :character
## Median : 1.000 Mode :character Median :1.000 Mode :character
## Mean : 1.967 Mean :1.385
## 3rd Qu.: 2.000 3rd Qu.:1.000
## Max. :10.000 Max. :4.000
## NA's :17
## c.collar.applied airway.insert suction ETT
## Length:30 Length:30 Length:30 Length:30
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## oxygensupp ventilator chesttube directpressure
## Length:30 Length:30 Length:30 Length:30
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## tourniquet splint pelvic.binder IV.Access
## Length:30 Length:30 Length:30 Length:30
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## IVF cardiac.monitoring resp.monitor neurological
## Length:30 Length:30 Length:30 Length:30
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
## print("Summary for Control Group:")## [1] "Summary for Control Group:"print(summary_control)## interv_ctrl code county age
## Length:30 Length:30 Length:30 Min. :27.00
## Class :character Class :character Class :character 1st Qu.:30.00
## Mode :character Mode :character Mode :character Median :34.50
## Mean :34.27
## 3rd Qu.:37.50
## Max. :44.00
##
## gender profession speciality experience
## Length:30 Length:30 Length:30 Min. : 3.000
## Class :character Class :character Class :character 1st Qu.: 4.000
## Mode :character Mode :character Mode :character Median : 6.500
## Mean : 6.933
## 3rd Qu.:10.000
## Max. :14.000
##
## deployment cont_education lastattended SOP
## Min. :1.000 Length:30 Min. :1 Length:30
## 1st Qu.:1.000 Class :character 1st Qu.:1 Class :character
## Median :2.000 Mode :character Median :1 Mode :character
## Mean :2.133 Mean :1
## 3rd Qu.:3.000 3rd Qu.:1
## Max. :5.000 Max. :1
## NA's :28
## c.collar.applied airway.insert suction ETT
## Length:30 Length:30 Length:30 Length:30
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## oxygensupp ventilator chesttube directpressure
## Length:30 Length:30 Length:30 Length:30
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## tourniquet splint pelvic.binder IV.Access
## Length:30 Length:30 Length:30 Length:30
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## IVF cardiac.monitoring resp.monitor neurological
## Length:30 Length:30 Length:30 Length:30
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
## # Frequency tables for categorical variables
cat_vars <- c("gender", "profession", "speciality", "SOP",
"c.collar.applied", "airway.insert", "suction", "ETT",
"oxygensupp", "ventilator", "chesttube", "directpressure",
"tourniquet", "splint", "pelvic.binder", "IV.Access",
"IVF", "cardiac.monitoring", "resp.monitor", "neurological")
for (var in cat_vars) {
freq_table <- hcp_data %>%
dplyr:::group_by.data.frame(interv_ctrl, !!sym(var)) %>%
summarise(count = n(), .groups = 'drop')
print(paste("Frequency table for", var))
print(freq_table)
}## [1] "Frequency table for gender"
## # A tibble: 4 × 3
## interv_ctrl gender count
## <chr> <chr> <int>
## 1 control female 16
## 2 control male 14
## 3 intervention female 14
## 4 intervention male 16
## [1] "Frequency table for profession"
## # A tibble: 9 × 3
## interv_ctrl profession count
## <chr> <chr> <int>
## 1 control bsn 4
## 2 control mo 4
## 3 control rco 6
## 4 control rn 16
## 5 intervention bsn 3
## 6 intervention mo 1
## 7 intervention paramedic 1
## 8 intervention rco 4
## 9 intervention rn 21
## [1] "Frequency table for speciality"
## # A tibble: 6 × 3
## interv_ctrl speciality count
## <chr> <chr> <int>
## 1 control none 30
## 2 intervention A&E 1
## 3 intervention critical care 1
## 4 intervention eye 1
## 5 intervention midwife 2
## 6 intervention none 25
## [1] "Frequency table for SOP"
## # A tibble: 3 × 3
## interv_ctrl SOP count
## <chr> <chr> <int>
## 1 control no 30
## 2 intervention no 19
## 3 intervention yes 11
## [1] "Frequency table for c.collar.applied"
## # A tibble: 5 × 3
## interv_ctrl c.collar.applied count
## <chr> <chr> <int>
## 1 control "correctly performed" 6
## 2 control "incorrectly performed" 16
## 3 control "not performed " 8
## 4 intervention "correctly performed" 20
## 5 intervention "incorrectly performed" 10
## [1] "Frequency table for airway.insert"
## # A tibble: 6 × 3
## interv_ctrl airway.insert count
## <chr> <chr> <int>
## 1 control correctly performed 4
## 2 control incorrectly performed 17
## 3 control not performed 9
## 4 intervention correctly performed 15
## 5 intervention incorrectly performed 11
## 6 intervention not performed 4
## [1] "Frequency table for suction"
## # A tibble: 4 × 3
## interv_ctrl suction count
## <chr> <chr> <int>
## 1 control correctly performed 27
## 2 control incorrectly performed 3
## 3 intervention correctly performed 27
## 4 intervention incorrectly performed 3
## [1] "Frequency table for ETT"
## # A tibble: 5 × 3
## interv_ctrl ETT count
## <chr> <chr> <int>
## 1 control incorrectly performed 4
## 2 control not performed 26
## 3 intervention correctly performed 8
## 4 intervention incorrectly performed 12
## 5 intervention not performed 10
## [1] "Frequency table for oxygensupp"
## # A tibble: 2 × 3
## interv_ctrl oxygensupp count
## <chr> <chr> <int>
## 1 control correctly performed 30
## 2 intervention correctly performed 30
## [1] "Frequency table for ventilator"
## # A tibble: 6 × 3
## interv_ctrl ventilator count
## <chr> <chr> <int>
## 1 control correctly performed 4
## 2 control incorrectly performed 10
## 3 control not performed 16
## 4 intervention correctly performed 10
## 5 intervention incorrectly performed 16
## 6 intervention not performed 4
## [1] "Frequency table for chesttube"
## # A tibble: 6 × 3
## interv_ctrl chesttube count
## <chr> <chr> <int>
## 1 control correctly performed 3
## 2 control incorrectly performed 16
## 3 control not performed 11
## 4 intervention correctly performed 5
## 5 intervention incorrectly performed 17
## 6 intervention not performed 8
## [1] "Frequency table for directpressure"
## # A tibble: 2 × 3
## interv_ctrl directpressure count
## <chr> <chr> <int>
## 1 control correctly performed 30
## 2 intervention correctly performed 30
## [1] "Frequency table for tourniquet"
## # A tibble: 6 × 3
## interv_ctrl tourniquet count
## <chr> <chr> <int>
## 1 control correctly performed 5
## 2 control incorrectly performed 21
## 3 control not performed 4
## 4 intervention correctly performed 4
## 5 intervention incorrectly performed 21
## 6 intervention not performed 5
## [1] "Frequency table for splint"
## # A tibble: 4 × 3
## interv_ctrl splint count
## <chr> <chr> <int>
## 1 control correctly performed 14
## 2 control incorrectly performed 16
## 3 intervention correctly performed 20
## 4 intervention incorrectly performed 10
## [1] "Frequency table for pelvic.binder"
## # A tibble: 5 × 3
## interv_ctrl pelvic.binder count
## <chr> <chr> <int>
## 1 control incorrectly performed 16
## 2 control not performed 14
## 3 intervention correctly performed 4
## 4 intervention incorrectly performed 17
## 5 intervention not performed 9
## [1] "Frequency table for IV.Access"
## # A tibble: 2 × 3
## interv_ctrl IV.Access count
## <chr> <chr> <int>
## 1 control correctly performed 30
## 2 intervention correctly performed 30
## [1] "Frequency table for IVF"
## # A tibble: 4 × 3
## interv_ctrl IVF count
## <chr> <chr> <int>
## 1 control correctly performed 22
## 2 control incorrectly performed 8
## 3 intervention correctly performed 27
## 4 intervention incorrectly performed 3
## [1] "Frequency table for cardiac.monitoring"
## # A tibble: 2 × 3
## interv_ctrl cardiac.monitoring count
## <chr> <chr> <int>
## 1 control correctly performed 30
## 2 intervention correctly performed 30
## [1] "Frequency table for resp.monitor"
## # A tibble: 2 × 3
## interv_ctrl resp.monitor count
## <chr> <chr> <int>
## 1 control correctly performed 30
## 2 intervention correctly performed 30
## [1] "Frequency table for neurological"
## # A tibble: 3 × 3
## interv_ctrl neurological count
## <chr> <chr> <int>
## 1 control correctly performed 12
## 2 control incorrectly performed 18
## 3 intervention correctly performed 30# Visualize differences between intervention and control groups
for (var in cat_vars) {
plot <- ggplot(hcp_data, aes(x = !!sym(var), fill = interv_ctrl)) +
geom_bar(position = "dodge") +
labs(title = paste("Distribution of", var, "by Group"),
x = var, y = "Count") +
theme_minimal()
print(plot)
}
# Create a dummy variable for SOP (yes = 1, no = 0)
hcp_data$SOP_dummy <- ifelse(hcp_data$SOP == "yes", 1, 0)
# Run the logistic regression
logit_model <- glm(SOP_dummy ~ interv_ctrl + age + county + gender + profession + speciality, data = hcp_data, family = binomial)
# Summarize the model
summary(logit_model)##
## Call:
## glm(formula = SOP_dummy ~ interv_ctrl + age + county + gender +
## profession + speciality, family = binomial, data = hcp_data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.698e+14 2.835e+14 0.599 0.549
## interv_ctrlintervention -1.698e+14 2.835e+14 -0.599 0.549
## age -9.920e-02 1.363e-01 -0.728 0.467
## countymigori -1.698e+14 2.835e+14 -0.599 0.549
## countysiaya -1.698e+14 2.835e+14 -0.599 0.549
## countyvihiga -4.781e+01 2.574e+05 0.000 1.000
## gendermale -1.508e+00 1.780e+00 -0.847 0.397
## professionmo 2.014e+01 4.647e+05 0.000 1.000
## professionparamedic 1.958e+01 1.059e+06 0.000 1.000
## professionrco 2.278e+01 1.775e+04 0.001 0.999
## professionrn 6.206e-01 2.058e+00 0.302 0.763
## specialitycritical care 4.667e+01 3.077e+05 0.000 1.000
## specialityeye -2.986e-01 1.706e+05 0.000 1.000
## specialitymidwife -2.517e+01 1.210e+05 0.000 1.000
## specialitynone -2.484e+01 1.210e+05 0.000 1.000
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 57.169 on 59 degrees of freedom
## Residual deviance: 14.356 on 45 degrees of freedom
## AIC: 44.356
##
## Number of Fisher Scoring iterations: 25### TODAY 24, 03, 2025
# Function to generate summary tables with p-values and handle missing values
generate_summary_table <- function(df, outcome_var) {
df_filtered <- df %>% filter(!is.na(.data[[outcome_var]])) # Remove NA values
if (nrow(df_filtered) == 0) {
return(NULL) # Return NULL if no valid data for the outcome
}
summary_table <- df_filtered %>%
group_by(MOI, treat_control) %>%
summarise(
count = n(),
percentage = round(100 * n() / nrow(df_filtered), 1),
outcome_mean = signif(mean(.data[[outcome_var]], na.rm = TRUE), 3), # Round to 3 significant figures
.groups = "drop"
) %>%
mutate(
n_pct = paste0(count, " (", percentage, "%)"),
outcome = outcome_var # Add outcome name as a column
)
# Check if there are enough valid observations for a t-test
if (length(unique(df_filtered$treat_control)) > 1) {
t_test_result <- t.test(df_filtered[[outcome_var]] ~ df_filtered$treat_control, na.action = na.omit)
p_value <- signif(t_test_result$p.value, 3) # Round p-value to 3 significant figures
} else {
p_value <- NA # Not enough data for t-test
}
summary_table <- summary_table %>%
mutate(p_value = p_value)
return(summary_table)
}
# Generate tables for each outcome variable
los_summary <- generate_summary_table(survey_data1, "LOSdays")
icu_summary <- generate_summary_table(survey_data1, "ICU_LOS")
ward_summary <- generate_summary_table(survey_data1, "ward_LOS")
# Combine all tables into one
combined_summary <- bind_rows(los_summary, icu_summary, ward_summary)
# Filter only significant results (p < 0.05) and non-missing p-values
significant_results <- combined_summary %>%
filter(!is.na(p_value) & p_value < 0.05)
# Print final filtered table
print("Filtered Significant Results:")## [1] "Filtered Significant Results:"print(significant_results)## # A tibble: 10 × 8
## MOI treat_control count percentage outcome_mean n_pct outcome p_value
## <fct> <fct> <int> <dbl> <dbl> <chr> <chr> <dbl>
## 1 assault control 1 4.3 12 1 (4… ICU_LOS 8.04e-5
## 2 assault treatment 1 4.3 2 1 (4… ICU_LOS 8.04e-5
## 3 burns control 1 4.3 7 1 (4… ICU_LOS 8.04e-5
## 4 fall control 4 17.4 10.5 4 (1… ICU_LOS 8.04e-5
## 5 gunshot treatment 1 4.3 2 1 (4… ICU_LOS 8.04e-5
## 6 penetratin… control 1 4.3 12 1 (4… ICU_LOS 8.04e-5
## 7 rti control 6 26.1 9.67 6 (2… ICU_LOS 8.04e-5
## 8 rti treatment 4 17.4 5 4 (1… ICU_LOS 8.04e-5
## 9 stab control 2 8.7 6 2 (8… ICU_LOS 8.04e-5
## 10 stab treatment 2 8.7 3 2 (8… ICU_LOS 8.04e-5# Load required libraries
library(ggplot2)
library(dplyr)
library(tidyr)
# Filter data to include only MOI values found in both groups
valid_moi <- survey_data1 %>%
filter(!is.na(ICU_LOS) & !is.na(LOSdays) & !is.na(ward_LOS)) %>%
group_by(MOI, treat_control) %>%
summarise(n = n(), .groups = "drop") %>%
count(MOI) %>%
filter(n == 2) %>%
pull(MOI)
# Aggregate the mean values for each outcome
heatmap_data <- survey_data1 %>%
filter(MOI %in% valid_moi) %>%
group_by(MOI, treat_control) %>%
summarise(
LOSdays_mean = round(mean(LOSdays, na.rm = TRUE), 3),
ICU_LOS_mean = round(mean(ICU_LOS, na.rm = TRUE), 3),
ward_LOS_mean = round(mean(ward_LOS, na.rm = TRUE), 3),
.groups = "drop"
) %>%
pivot_longer(cols = c(LOSdays_mean, ICU_LOS_mean, ward_LOS_mean),
names_to = "Outcome", values_to = "Mean_Value") %>%
drop_na(Mean_Value) # Ensure no missing values
# Define color scale for heatmap
max_val <- max(heatmap_data$Mean_Value, na.rm = TRUE)
min_val <- min(heatmap_data$Mean_Value, na.rm = TRUE)
# Create heatmap
ggplot(heatmap_data, aes(x = treat_control, y = MOI, fill = Mean_Value)) +
geom_tile(color = "white", size = 0.6) +
scale_fill_gradientn(
colors = c("#2166AC", "#67A9CF", "#D1E5F0", "#FDDDBC", "#D6604D", "#B2182B"),
limits = c(min_val, max_val),
name = "Mean Value"
) +
facet_wrap(~ Outcome, ncol = 1) + # Combine all outcomes in one heatmap
labs(
title = "🔥 Heatmap of Length of Stay by MOI & Treatment Group",
subtitle = "Darker Red = Higher Mean Value | Darker Blue = Lower Mean Value",
x = "Group",
y = "Mechanism of Injury (MOI)"
) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(face = "bold", size = 16),
plot.subtitle = element_text(size = 12, face = "italic"),
axis.text.x = element_text(face = "bold", size = 12),
axis.text.y = element_text(face = "bold", size = 12),
legend.title = element_text(face = "bold"),
strip.text = element_text(face = "bold", size = 14),
legend.position = "right"
)
# Load required libraries
library(dplyr)
library(broom)
library(ggpubr)##
## Attaching package: 'ggpubr'## The following objects are masked from 'package:flextable':
##
## border, font, rotate# Define categorical variables
categorical_vars <- c("gender", "mstatus", "edulevel")
outcomes <- c("LOSdays", "ICU_LOS", "ward_LOS")
# Ensure categorical variables are factors
survey_data1 <- survey_data1 %>%
mutate(across(all_of(categorical_vars), as.factor))
# Function for ANOVA or Kruskal-Wallis Test (categorical vs numeric)
run_anova_or_kruskal <- function(var, outcome_var) {
if (length(unique(survey_data1[[var]])) > 1) {
# Check normality with Shapiro-Wilk test
p_normality <- shapiro.test(survey_data1[[outcome_var]])$p.value
if (p_normality > 0.05) {
# Use ANOVA if normality assumption is met
test_result <- aov(as.formula(paste(outcome_var, "~", var)), data = survey_data1, na.action = na.omit)
p_value <- summary(test_result)[[1]][["Pr(>F)"]][1]
test_name <- "ANOVA"
} else {
# Use Kruskal-Wallis test if data is non-normal
test_result <- kruskal.test(as.formula(paste(outcome_var, "~", var)), data = survey_data1)
p_value <- test_result$p.value
test_name <- "Kruskal-Wallis"
}
return(data.frame(Variable = var, Outcome = outcome_var, Test = test_name, P_Value = p_value))
} else {
return(data.frame(Variable = var, Outcome = outcome_var, Test = "Skipped (only 1 level)", P_Value = NA))
}
}
# Run ANOVA/Kruskal-Wallis for each categorical variable and outcome
anova_kruskal_results <- bind_rows(lapply(categorical_vars, function(var) {
bind_rows(lapply(outcomes, function(outcome_var) run_anova_or_kruskal(var, outcome_var)))
}))
# Convert outcome to binary (1 = alive, 0 = dead)
survey_data1 <- survey_data1 %>%
mutate(outcome_binary = ifelse(outcome == "alive", 1, 0))
# Logistic regression for categorical variables predicting mortality
run_logistic_regression <- function(var) {
model <- glm(outcome_binary ~ survey_data1[[var]], data = survey_data1, family = binomial, na.action = na.omit)
model_summary <- broom::tidy(model)
data.frame(Variable = var, Outcome = "Mortality", Test = "Logistic Regression", P_Value = model_summary$p.value[2])
}
# Run logistic regression for categorical variables
logit_results <- bind_rows(lapply(categorical_vars, run_logistic_regression))
# Combine all results into a single table
association_table <- bind_rows(anova_kruskal_results, logit_results) %>%
arrange(P_Value)
# Print results
print(association_table)## Variable Outcome Test P_Value
## 1 mstatus ward_LOS ANOVA 0.2672840
## 2 gender ICU_LOS ANOVA 0.3061479
## 3 gender ward_LOS ANOVA 0.3848395
## 4 edulevel LOSdays Kruskal-Wallis 0.6256050
## 5 edulevel ward_LOS ANOVA 0.6457966
## 6 edulevel Mortality Logistic Regression 0.6743351
## 7 edulevel ICU_LOS ANOVA 0.7470896
## 8 gender Mortality Logistic Regression 0.7491159
## 9 mstatus ICU_LOS ANOVA 0.8198396
## 10 gender LOSdays Kruskal-Wallis 0.8219949
## 11 mstatus LOSdays Kruskal-Wallis 0.9712521
## 12 mstatus Mortality Logistic Regression 0.9862287# Load required libraries
library(dplyr)
# Ensure MOI, time, treat_control, and Admit are factors
survey_data1 <- survey_data1 %>%
mutate(across(c("MOI", "time", "treat_control", "Admit"), as.factor)) %>%
mutate(ICU_Admit = ifelse(Admit == "icu", 1, 0)) # Create binary ICU admission variable
# Filter for endline data only
endline_data <- survey_data1 %>%
filter(time == "endline") %>%
drop_na(MOI, treat_control, Admit) # Remove NAs in key variables
# Compute ICU admission proportions by MOI category and treatment group
moi_descriptive_stats <- endline_data %>%
group_by(MOI, treat_control) %>%
summarise(
Total = n(),
ICU_Admissions = sum(ICU_Admit),
ICU_Proportion = round(mean(ICU_Admit) * 100, 1), # Convert to percentage
.groups = "drop"
)
# Count ICU vs. Other admissions at endline by treatment group
admit_counts <- endline_data %>%
group_by(treat_control, Admit) %>%
summarise(Count = n(), .groups = "drop")
# Print results
print(moi_descriptive_stats)## # A tibble: 17 × 5
## MOI treat_control Total ICU_Admissions ICU_Proportion
## <fct> <fct> <int> <dbl> <dbl>
## 1 assault control 5 0 0
## 2 assault treatment 3 1 33.3
## 3 burns control 3 1 33.3
## 4 burns treatment 5 0 0
## 5 crushing control 1 0 0
## 6 cut control 1 0 0
## 7 cut treatment 2 0 0
## 8 fall control 8 5 62.5
## 9 fall treatment 9 0 0
## 10 gunshot control 3 0 0
## 11 gunshot treatment 2 1 50
## 12 penetrating inj control 1 1 100
## 13 rti control 34 6 17.6
## 14 rti treatment 24 4 16.7
## 15 stab control 8 2 25
## 16 stab treatment 7 2 28.6
## 17 stab/ cut treatment 1 0 0print(admit_counts)## # A tibble: 6 × 3
## treat_control Admit Count
## <fct> <fct> <int>
## 1 control icu 15
## 2 control OT 15
## 3 control ward 34
## 4 treatment icu 8
## 5 treatment OT 17
## 6 treatment ward 28# Load required libraries
library(dplyr)
library(ggplot2)
# Ensure MOI, time, treat_control, and Admit are factors
survey_data1 <- survey_data1 %>%
mutate(across(c("MOI", "time", "treat_control", "Admit"), as.factor)) %>%
mutate(ICU_Admit = ifelse(Admit == "icu", 1, 0)) # Create binary ICU admission variable
# Filter for endline data only
endline_data <- survey_data1 %>%
filter(time == "endline") %>%
drop_na(treat_control, Admit) # Remove NAs in key variables
# Count ICU vs. Other admissions at endline by treatment group
admit_counts <- endline_data %>%
group_by(treat_control, Admit) %>%
summarise(Count = n(), .groups = "drop")
# Create a bar plot
ggplot(admit_counts, aes(x = treat_control, y = Count, fill = Admit)) +
geom_bar(stat = "identity", position = "dodge") +
labs(
title = "ICU vs. Other Admissions by Treatment Group (Endline)",
x = "Treatment Group",
y = "Number of Admissions",
fill = "Admission Type"
) +
theme_minimal()
# Load required libraries
library(dplyr)
# Ensure relevant variables are in correct formats
survey_data1 <- survey_data1 %>%
mutate(across(c("MOI", "time"), as.factor))
# Filter for endline data only
endline_data <- survey_data1 %>%
filter(time == "endline")
# Handle missing values: Replace missing LOSdays with group mean (imputation)
endline_data <- endline_data %>%
group_by(MOI) %>%
mutate(LOSdays = ifelse(is.na(LOSdays), mean(LOSdays, na.rm = TRUE), LOSdays)) %>%
ungroup()
# Compute Mean, Min, and Max for LOSdays by MOI, removing groups with NaN values
los_summary <- endline_data %>%
group_by(MOI) %>%
summarise(
Mean_LOS = mean(LOSdays, na.rm = TRUE),
Min_LOS = min(LOSdays, na.rm = TRUE),
Max_LOS = max(LOSdays, na.rm = TRUE),
.groups = "drop"
) %>%
filter(!is.nan(Mean_LOS)) # Remove rows where mean is NaN
# Print results
print(los_summary)## # A tibble: 10 × 4
## MOI Mean_LOS Min_LOS Max_LOS
## <fct> <dbl> <int> <int>
## 1 assault 27.5 14 62
## 2 burns 28.2 16 50
## 3 crushing 26 26 26
## 4 cut 18.3 15 23
## 5 fall 22.1 1 34
## 6 gunshot 11 1 23
## 7 penetrating inj 33 33 33
## 8 rti 22.0 1 58
## 9 stab 17.9 1 30
## 10 stab/ cut 14 14 14# Load required libraries
library(dplyr)
# Ensure relevant variables are factors
survey_data1 <- survey_data1 %>%
mutate(across(c("MOI", "time", "state.1", "outcome"), as.factor))
# Filter for endline data only
endline_data <- survey_data1 %>%
filter(time == "endline")
# Compute frequency & proportion of state.1 by MOI
state_summary <- endline_data %>%
group_by(MOI, state.1) %>%
summarise(Count = n(), .groups = "drop") %>%
group_by(MOI) %>%
mutate(Proportion = Count / sum(Count) * 100) %>%
ungroup()
# Compute frequency & proportion of outcome by MOI
outcome_summary <- endline_data %>%
group_by(MOI, outcome) %>%
summarise(Count = n(), .groups = "drop") %>%
group_by(MOI) %>%
mutate(Proportion = Count / sum(Count) * 100) %>%
ungroup()
# Print results
print(state_summary)## # A tibble: 19 × 4
## MOI state.1 Count Proportion
## <fct> <fct> <int> <dbl>
## 1 assault stable 8 100
## 2 burns rehab 1 12.5
## 3 burns stable 7 87.5
## 4 crushing stable 1 100
## 5 cut stable 3 100
## 6 fall dead 1 5.88
## 7 fall fair 1 5.88
## 8 fall stable 15 88.2
## 9 gunshot dead 2 40
## 10 gunshot stable 3 60
## 11 penetrating inj stable 1 100
## 12 rti dead 2 3.45
## 13 rti fair 4 6.90
## 14 rti rehab 2 3.45
## 15 rti stable 50 86.2
## 16 stab crutches 1 6.67
## 17 stab dead 3 20
## 18 stab stable 11 73.3
## 19 stab/ cut stable 1 100print(outcome_summary)## # A tibble: 14 × 4
## MOI outcome Count Proportion
## <fct> <fct> <int> <dbl>
## 1 assault alive 8 100
## 2 burns alive 8 100
## 3 crushing alive 1 100
## 4 cut alive 3 100
## 5 fall alive 16 94.1
## 6 fall dead 1 5.88
## 7 gunshot alive 3 60
## 8 gunshot dead 2 40
## 9 penetrating inj alive 1 100
## 10 rti alive 56 96.6
## 11 rti dead 2 3.45
## 12 stab alive 12 80
## 13 stab dead 3 20
## 14 stab/ cut alive 1 100