Good Friday!
Simon Analysis
2024-03-29
[1] "C:/Users/Afg-iMMAP-75/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\Afg-iMMAP-75\AppData\Local\Temp\Rtmpm0AkZV\downloaded_packages Stratified by treat_control
level control treatment p test
n 110 110
ICU_admissions (mean (SD)) 0.19 (0.39) 0.20 (0.40) 0.866
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"
==================================================================================
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
59 161
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 41 57 71.930
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.48 | 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 | 5995.000000 | 0.8665146 |
| Independent t-test | ICU_LOS | 4.926856 | 0.0000804 |
| Mann-Whitney U | hosp_admission | 6435.000000 | 0.2884663 |
| 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.93611
% Var explained: 13.78
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.635e+01 2.785e+03 -0.006 0.995
age -2.048e-02 1.276e-02 -1.605 0.108
gendermale 3.204e-01 3.857e-01 0.831 0.406
MOIburns 7.562e-01 1.026e+00 0.737 0.461
MOIcrushing -1.519e+01 3.956e+03 -0.004 0.997
MOIcut -1.502e+01 1.970e+03 -0.008 0.994
MOIfall 5.375e-01 9.164e-01 0.586 0.558
MOIgunshot 9.890e-01 1.374e+00 0.720 0.472
MOIpenetrating inj 2.022e+01 3.956e+03 0.005 0.996
MOIrti 8.110e-01 8.089e-01 1.003 0.316
MOIrti, motorbike -1.602e+01 2.797e+03 -0.006 0.995
MOIrti,fall -1.458e+01 3.956e+03 -0.004 0.997
MOIstab 9.650e-01 8.966e-01 1.076 0.282
MOIstab, assault -1.575e+01 3.956e+03 -0.004 0.997
MOIstab/ cut -1.625e+01 3.956e+03 -0.004 0.997
triagecategoryimmediate 1.563e+01 2.785e+03 0.006 0.996
triagecategoryurgent 1.555e+01 2.785e+03 0.006 0.996
Edduration_mins -1.451e-02 7.745e-03 -1.873 0.061 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 217.38 on 219 degrees of freedom
Residual deviance: 197.87 on 202 degrees of freedom
AIC: 233.87
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.783e+00 1.787e+00 -0.998 0.3183
age -1.865e-02 1.135e-02 -1.643 0.1004
gendermale 1.803e-01 3.511e-01 0.514 0.6076
MOIburns 1.710e+01 1.029e+03 0.017 0.9867
MOIcrushing 1.614e+01 3.956e+03 0.004 0.9967
MOIcut 9.084e-01 1.297e+00 0.700 0.4838
MOIfall 1.396e+00 6.711e-01 2.081 0.0374 *
MOIgunshot 3.047e-01 1.071e+00 0.285 0.7760
MOIpenetrating inj 1.733e+01 3.956e+03 0.004 0.9965
MOIrti 9.277e-01 5.384e-01 1.723 0.0848 .
MOIrti, motorbike 1.785e+01 2.791e+03 0.006 0.9949
MOIrti,fall 1.613e+01 3.956e+03 0.004 0.9967
MOIstab 2.961e-01 6.464e-01 0.458 0.6469
MOIstab, assault 1.717e+01 3.956e+03 0.004 0.9965
MOIstab/ cut -1.713e+01 3.956e+03 -0.004 0.9965
triagecategoryimmediate 1.455e+00 1.564e+00 0.930 0.3524
triagecategoryurgent 2.133e+00 1.559e+00 1.368 0.1713
Edduration_mins 1.305e-02 7.374e-03 1.770 0.0767 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 255.83 on 219 degrees of freedom
Residual deviance: 221.46 on 202 degrees of freedom
AIC: 257.46
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.623, df = 1, p-value = 0.4299 Df Sum Sq Mean Sq F value Pr(>F)
ICU_admissions 1 4977 4977 3.593 0.0593 .
Residuals 218 301930 1385
---
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.635e+01 2.785e+03 -0.006 0.995
age -2.048e-02 1.276e-02 -1.605 0.108
gendermale 3.204e-01 3.857e-01 0.831 0.406
MOIburns 7.562e-01 1.026e+00 0.737 0.461
MOIcrushing -1.519e+01 3.956e+03 -0.004 0.997
MOIcut -1.502e+01 1.970e+03 -0.008 0.994
MOIfall 5.375e-01 9.164e-01 0.586 0.558
MOIgunshot 9.890e-01 1.374e+00 0.720 0.472
MOIpenetrating inj 2.022e+01 3.956e+03 0.005 0.996
MOIrti 8.110e-01 8.089e-01 1.003 0.316
MOIrti, motorbike -1.602e+01 2.797e+03 -0.006 0.995
MOIrti,fall -1.458e+01 3.956e+03 -0.004 0.997
MOIstab 9.650e-01 8.966e-01 1.076 0.282
MOIstab, assault -1.575e+01 3.956e+03 -0.004 0.997
MOIstab/ cut -1.625e+01 3.956e+03 -0.004 0.997
triagecategoryimmediate 1.563e+01 2.785e+03 0.006 0.996
triagecategoryurgent 1.555e+01 2.785e+03 0.006 0.996
Edduration_mins -1.451e-02 7.745e-03 -1.873 0.061 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 217.38 on 219 degrees of freedom
Residual deviance: 197.87 on 202 degrees of freedom
AIC: 233.87
Number of Fisher Scoring iterations: 16 (Intercept) age gendermale
7.950773e-08 9.797247e-01 1.377658e+00
MOIburns MOIcrushing MOIcut
2.130117e+00 2.535836e-07 3.002222e-07
MOIfall MOIgunshot MOIpenetrating inj
1.711702e+00 2.688601e+00 6.028755e+08
MOIrti MOIrti, motorbike MOIrti,fall
2.250114e+00 1.105878e-07 4.633237e-07
MOIstab MOIstab, assault MOIstab/ cut
2.624905e+00 1.451573e-07 8.751255e-08
triagecategoryimmediate triagecategoryurgent Edduration_mins
6.147839e+06 5.665643e+06 9.855971e-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.800000e-01 | 0.108 |
| gendermale | gendermale | 1.378000e+00 | 0.406 |
| MOIburns | MOIburns | 2.130000e+00 | 0.461 |
| MOIcrushing | MOIcrushing | 0.000000e+00 | 0.997 |
| MOIcut | MOIcut | 0.000000e+00 | 0.994 |
| MOIfall | MOIfall | 1.712000e+00 | 0.558 |
| MOIgunshot | MOIgunshot | 2.689000e+00 | 0.472 |
| MOIpenetrating inj | MOIpenetrating inj | 6.028755e+08 | 0.996 |
| MOIrti | MOIrti | 2.250000e+00 | 0.316 |
| MOIrti, motorbike | MOIrti, motorbike | 0.000000e+00 | 0.995 |
| MOIrti,fall | MOIrti,fall | 0.000000e+00 | 0.997 |
| MOIstab | MOIstab | 2.625000e+00 | 0.282 |
| MOIstab, assault | MOIstab, assault | 0.000000e+00 | 0.997 |
| MOIstab/ cut | MOIstab/ cut | 0.000000e+00 | 0.997 |
| triagecategoryimmediate | triagecategoryimmediate | 6.147839e+06 | 0.996 |
| triagecategoryurgent | triagecategoryurgent | 5.665643e+06 | 0.996 |
| Edduration_mins | Edduration_mins | 9.860000e-01 | 0.061 |
# 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.108 |
| gendermale | 1.38 | 0.406 |
| MOIburns | 2.13 | 0.461 |
| MOIcrushing | 0.00 | 0.997 |
| MOIcut | 0.00 | 0.994 |
| MOIfall | 1.71 | 0.558 |
| MOIgunshot | 2.69 | 0.472 |
| MOIpenetrating inj | 602,875,534.81 | 0.996 |
| MOIrti | 2.25 | 0.316 |
| MOIrti, motorbike | 0.00 | 0.995 |
| MOIrti,fall | 0.00 | 0.997 |
| MOIstab | 2.62 | 0.282 |
| MOIstab, assault | 0.00 | 0.997 |
| MOIstab/ cut | 0.00 | 0.997 |
| triagecategoryimmediate | 6,147,838.81 | 0.996 |
| triagecategoryurgent | 5,665,642.95 | 0.996 |
| Edduration_mins | 0.99 | 0.061 |
# Load data
survey_data <- read_csv("C:/Users/Afg-iMMAP-75/Downloads/survey_data1.csv")## Rows: 220 Columns: 53
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (33): county, time, treat_control, gender, mstatus, edulevel, occup, dat...
## dbl (20): code, age, ISS, Durationmins, HRinitial, SBPinitial, DBPinitial, R...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.# Convert categorical variables to factors
survey_data <- survey_data %>%
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 39.7 9.07 66.7
## 2 control endline siaya 32.3 5.13 66.7
## # ℹ 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 | 39.66667 | 9.073772 | 66.66667 | 33.33333 | 52.27273 | 19.54545 | 26.36364 | 0 |
| control | endline | siaya | 32.33333 | 5.131601 | 66.66667 | 33.33333 | 52.27273 | 19.54545 | 26.36364 | 0 |
# 2. Intervention Impact Analysis (Controlling for Confounders)
# -------------------------------------------------------------
# Logistic Regression for mortality controlling for confounders
mortality_model <- glm(survey_data1$final_outcome ~ survey_data1$treat_control + survey_data1$gender + survey_data1$MOI +
survey_data1$HRinitial + survey_data1$SBPinitial + survey_data1$DBPinitial + survey_data1$RRinitial + survey_data1$SPO2initial,
data = survey_data, family = binomial)
tidy(mortality_model )## # A tibble: 20 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 9.83 1452. 0.00677 0.995
## 2 survey_data1$treat_controltreatment -0.436 0.476 -0.915 0.360
## 3 survey_data1$gendermale 0.301 0.511 0.589 0.556
## 4 survey_data1$MOIburns -15.7 1452. -0.0108 0.991
## 5 survey_data1$MOIcrushing -0.658 6682. -0.0000984 1.00
## 6 survey_data1$MOIcut 0.183 3447. 0.0000531 1.00
## 7 survey_data1$MOIfall -15.2 1452. -0.0105 0.992
## 8 survey_data1$MOIgunshot -19.0 1452. -0.0131 0.990
## 9 survey_data1$MOIpenetrating inj -0.491 6682. -0.0000735 1.00
## 10 survey_data1$MOIrti -16.5 1452. -0.0114 0.991
## 11 survey_data1$MOIrti, motorbike -0.300 4835. -0.0000619 1.00
## 12 survey_data1$MOIrti,fall -0.0621 6682. -0.00000930 1.00
## 13 survey_data1$MOIstab -16.9 1452. -0.0116 0.991
## 14 survey_data1$MOIstab, assault -0.569 6682. -0.0000852 1.00
## 15 survey_data1$MOIstab/ cut -0.747 6682. -0.000112 1.00
## 16 survey_data1$HRinitial -0.00405 0.0164 -0.247 0.805
## 17 survey_data1$SBPinitial -0.0151 0.0176 -0.858 0.391
## 18 survey_data1$DBPinitial -0.00637 0.0285 -0.223 0.823
## 19 survey_data1$RRinitial -0.120 0.0509 -2.36 0.0180
## 20 survey_data1$SPO2initial 0.147 0.0473 3.11 0.00184# Extract tidy results and calculate odds ratios
odds_ratios <- tidy(mortality_model) %>%
mutate(odds_ratio = exp(estimate)) # Exponentiate to get odds ratios
# View the odds ratios
print(odds_ratios)## # A tibble: 20 × 6
## term estimate std.error statistic p.value odds_ratio
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 9.83e+0 1452. 6.77e-3 0.995 1.86e+4
## 2 survey_data1$treat_controltr… -4.36e-1 0.476 -9.15e-1 0.360 6.47e-1
## 3 survey_data1$gendermale 3.01e-1 0.511 5.89e-1 0.556 1.35e+0
## 4 survey_data1$MOIburns -1.57e+1 1452. -1.08e-2 0.991 1.45e-7
## 5 survey_data1$MOIcrushing -6.58e-1 6682. -9.84e-5 1.00 5.18e-1
## 6 survey_data1$MOIcut 1.83e-1 3447. 5.31e-5 1.00 1.20e+0
## 7 survey_data1$MOIfall -1.52e+1 1452. -1.05e-2 0.992 2.56e-7
## 8 survey_data1$MOIgunshot -1.90e+1 1452. -1.31e-2 0.990 5.52e-9
## 9 survey_data1$MOIpenetrating … -4.91e-1 6682. -7.35e-5 1.00 6.12e-1
## 10 survey_data1$MOIrti -1.65e+1 1452. -1.14e-2 0.991 6.62e-8
## 11 survey_data1$MOIrti, motorbi… -3.00e-1 4835. -6.19e-5 1.00 7.41e-1
## 12 survey_data1$MOIrti,fall -6.21e-2 6682. -9.30e-6 1.00 9.40e-1
## 13 survey_data1$MOIstab -1.69e+1 1452. -1.16e-2 0.991 4.63e-8
## 14 survey_data1$MOIstab, assault -5.69e-1 6682. -8.52e-5 1.00 5.66e-1
## 15 survey_data1$MOIstab/ cut -7.47e-1 6682. -1.12e-4 1.00 4.74e-1
## 16 survey_data1$HRinitial -4.05e-3 0.0164 -2.47e-1 0.805 9.96e-1
## 17 survey_data1$SBPinitial -1.51e-2 0.0176 -8.58e-1 0.391 9.85e-1
## 18 survey_data1$DBPinitial -6.37e-3 0.0285 -2.23e-1 0.823 9.94e-1
## 19 survey_data1$RRinitial -1.20e-1 0.0509 -2.36e+0 0.0180 8.87e-1
## 20 survey_data1$SPO2initial 1.47e-1 0.0473 3.11e+0 0.00184 1.16e+0# Display Impact of Intervention on Mortality
# Load the required library
library(gtsummary)##
## Attaching package: 'gtsummary'## The following object is masked from 'package:flextable':
##
## continuous_summary## The following object is masked from 'package:MASS':
##
## selectmodel1 <- survey_data1 |> select(final_outcome, gender, MOI)
model1 |> tbl_summary()| Characteristic | N = 2201 |
|---|---|
| final_outcome | 194 (88%) |
| gender | |
| female | 91 (41%) |
| male | 129 (59%) |
| MOI | |
| assault | 18 (8.2%) |
| burns | 14 (6.4%) |
| crushing | 1 (0.5%) |
| cut | 4 (1.8%) |
| fall | 32 (15%) |
| gunshot | 5 (2.3%) |
| penetrating inj | 1 (0.5%) |
| rti | 111 (50%) |
| rti, motorbike | 2 (0.9%) |
| rti,fall | 1 (0.5%) |
| stab | 29 (13%) |
| stab, assault | 1 (0.5%) |
| stab/ cut | 1 (0.5%) |
| 1 n (%) | |
library(tidyr)
library(dplyr)
library(broom)
library(gt)
# Clean the data by removing rows with NA values in the relevant columns
survey_data1_clean <- survey_data1 %>%
drop_na(treat_control, gender, MOI, HRinitial, SBPinitial, DBPinitial, RRinitial, SPO2initial)
# Logistic Regression for mortality controlling for confounders
mortality_model <- glm(final_outcome ~ treat_control + gender + MOI +
HRinitial + SBPinitial + DBPinitial + RRinitial + SPO2initial,
data = survey_data1_clean, family = binomial)
# Extract tidy results and calculate odds ratios
mortality_results <- tidy(mortality_model)
# Exponentiate the estimate to get odds ratios
mortality_results <- mortality_results %>%
mutate(odds_ratio = exp(estimate))
# Calculate Confidence Intervals for odds ratios
conf_intervals <- confint(mortality_model)## Waiting for profiling to be done...## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredconf_intervals <- exp(conf_intervals) # Convert CI from log-odds to odds ratios
# Add Confidence Interval columns to the results
mortality_results <- mortality_results %>%
mutate(
CI_Lower = conf_intervals[, 1],
CI_Upper = conf_intervals[, 2]
)
# Display the results in a nice table format
mortality_results %>%
gt() %>%
tab_header(title = "Impact of Intervention on Mortality")| Impact of Intervention on Mortality | |||||||
| term | estimate | std.error | statistic | p.value | odds_ratio | CI_Lower | CI_Upper |
|---|---|---|---|---|---|---|---|
| (Intercept) | 9.830459042 | 1.452438e+03 | 6.768249e-03 | 0.994599760 | 1.859149e+04 | 1.716746e-23 | 4.489749e+211 |
| treat_controltreatment | -0.435823144 | 4.761837e-01 | -9.152417e-01 | 0.360064758 | 6.467321e-01 | 2.481482e-01 | 1.634543e+00 |
| gendermale | 0.300939618 | 5.113207e-01 | 5.885535e-01 | 0.556160819 | 1.351128e+00 | 4.881282e-01 | 3.711690e+00 |
| MOIburns | -15.747210329 | 1.452433e+03 | -1.084195e-02 | 0.991349541 | 1.449017e-07 | NA | 2.268984e+28 |
| MOIcrushing | -0.657756657 | 6.682393e+03 | -9.843131e-05 | 0.999921463 | 5.180121e-01 | 0.000000e+00 | 3.756866e+146 |
| MOIcut | 0.183025046 | 3.447001e+03 | 5.309690e-05 | 0.999957635 | 1.200844e+00 | 0.000000e+00 | 8.317789e-42 |
| MOIfall | -15.179964837 | 1.452433e+03 | -1.045141e-02 | 0.991661136 | 2.555201e-07 | 3.278474e-197 | 4.292790e+23 |
| MOIgunshot | -19.015283269 | 1.452433e+03 | -1.309202e-02 | 0.989554375 | 5.517818e-09 | 6.019299e-193 | 1.118945e+23 |
| MOIpenetrating inj | -0.491245473 | 6.682393e+03 | -7.351341e-05 | 0.999941345 | 6.118639e-01 | 0.000000e+00 | 6.928723e+125 |
| MOIrti | -16.530338854 | 1.452432e+03 | -1.138114e-02 | 0.990919359 | 6.621632e-08 | NA | 2.642809e+27 |
| MOIrti, motorbike | -0.299504773 | 4.834720e+03 | -6.194873e-05 | 0.999950572 | 7.411852e-01 | 0.000000e+00 | 2.744056e+96 |
| MOIrti,fall | -0.062143976 | 6.682393e+03 | -9.299659e-06 | 0.999992580 | 9.397476e-01 | 0.000000e+00 | 1.662122e+149 |
| MOIstab | -16.888115956 | 1.452432e+03 | -1.162747e-02 | 0.990722830 | 4.630037e-08 | NA | 7.249948e+27 |
| MOIstab, assault | -0.569223934 | 6.682393e+03 | -8.518265e-05 | 0.999932034 | 5.659645e-01 | 0.000000e+00 | 1.433008e+139 |
| MOIstab/ cut | -0.746647140 | 6.682393e+03 | -1.117335e-04 | 0.999910850 | 4.739530e-01 | 0.000000e+00 | 1.091618e+135 |
| HRinitial | -0.004046115 | 1.640810e-02 | -2.465925e-01 | 0.805223616 | 9.959621e-01 | 9.641578e-01 | 1.028500e+00 |
| SBPinitial | -0.015139689 | 1.763886e-02 | -8.583145e-01 | 0.390718836 | 9.849743e-01 | 9.507935e-01 | 1.019641e+00 |
| DBPinitial | -0.006366841 | 2.848883e-02 | -2.234856e-01 | 0.823157616 | 9.936534e-01 | 9.393564e-01 | 1.051247e+00 |
| RRinitial | -0.120343019 | 5.088664e-02 | -2.364924e+00 | 0.018033793 | 8.866163e-01 | 7.984703e-01 | 9.767569e-01 |
| SPO2initial | 0.147451372 | 4.734757e-02 | 3.114233e+00 | 0.001844235 | 1.158877e+00 | 1.060905e+00 | 1.281055e+00 |
## TODAY'S 07/02/25 Analysis
library(dplyr)
library(ggplot2)
# Subsetting relevant columns for analysis
relevant_data <- survey_data1 %>%
select(treat_control, MOI, EDduration = Durationmins,
LOSdays, outcome, OT_admissions, state.1, HRinitial, SBPinitial, DBPinitial
)
# Checking data structure
str(relevant_data)## 'data.frame': 220 obs. of 10 variables:
## $ treat_control: chr "control" "control" "control" "control" ...
## $ MOI : chr "rti" "rti,fall" "stab" "burns" ...
## $ EDduration : int 30 30 25 30 30 30 30 20 20 10 ...
## $ LOSdays : int 35 21 15 15 25 30 1 25 20 22 ...
## $ outcome : chr "alive" "alive" "alive" "alive" ...
## $ OT_admissions: num 0 0 0 1 0 1 1 0 0 0 ...
## $ state.1 : chr "fair" "stable" "stable" "rehab" ...
## $ HRinitial : int 98 94 90 105 98 102 122 110 98 118 ...
## $ SBPinitial : int 110 124 130 90 120 90 90 115 124 90 ...
## $ DBPinitial : int 80 72 80 70 70 60 75 68 75 65 ...# Descriptive statistics for intervention and control groups
summary_stats <- relevant_data %>%
group_by(treat_control) %>%
summarise(
avg_EDduration = mean(EDduration, na.rm = TRUE),
avg_LOS = mean(LOSdays, na.rm = TRUE),
avg_HRinitial = mean(HRinitial, na.rm = TRUE),
avg_SBPinitial = mean(SBPinitial, na.rm = TRUE),
avg_DBPinitial = mean(DBPinitial, na.rm = TRUE)
)
print(summary_stats)## # A tibble: 2 × 6
## treat_control avg_EDduration avg_LOS avg_HRinitial avg_SBPinitial
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 control 25 23.5 95.8 110.
## 2 treatment 18.8 21.2 99.8 105.
## # ℹ 1 more variable: avg_DBPinitial <dbl># Scatter plot for ED duration vs LOS for individual patients
ggplot(relevant_data, aes(x = EDduration, y = LOSdays, color = treat_control)) +
geom_point(alpha = 0.5, size = 2) + # Adding transparency for better visibility
geom_smooth(method = "lm", se = TRUE) + # Linear regression with confidence interval
theme_minimal() +
labs(title = "ED Duration vs Length of Stay (LOS) by Group",
x = "ED Duration (Minutes)",
y = "Length of Stay (Days)",
color = "Group") +
scale_color_manual(values = c("control" = "red", "intervention" = "blue"))## `geom_smooth()` using formula = 'y ~ x'
# Load dplyr
library(dplyr)
# Summary statistics by group
relevant_data %>%
group_by(treat_control) %>%
summarise(
Mean_EDduration = mean(EDduration, na.rm = TRUE),
SD_EDduration = sd(EDduration, na.rm = TRUE),
Mean_LOS = mean(LOSdays, na.rm = TRUE),
SD_LOS = sd(LOSdays, na.rm = TRUE),
n = n()
)## # A tibble: 2 × 6
## treat_control Mean_EDduration SD_EDduration Mean_LOS SD_LOS n
## <chr> <dbl> <dbl> <dbl> <dbl> <int>
## 1 control 25 9.11 23.5 12.0 110
## 2 treatment 18.8 12.3 21.2 13.7 110# Multiple regression to assess impact of multiple factors on LOS
# Compare mean LOS between control and intervention groups
anova_los <- aov(LOSdays ~ treat_control, data = relevant_data)
summary(anova_los)## Df Sum Sq Mean Sq F value Pr(>F)
## treat_control 1 293 293.2 1.762 0.186
## Residuals 218 36288 166.5# Post-hoc test if needed
posthoc_los <- TukeyHSD(anova_los)
print(posthoc_los)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = LOSdays ~ treat_control, data = relevant_data)
##
## $treat_control
## diff lwr upr p adj
## treatment-control -2.309091 -5.737865 1.119683 0.1857989# Create a binary outcome variable (e.g., LOS > 5 days)
relevant_data$outcome_los <- ifelse(relevant_data$LOSdays > 5, 1, 0)
# Logistic regression model
logistic_model <- glm(outcome_los ~ EDduration + treat_control + MOI,
family = binomial(link = "logit"), data = relevant_data)
summary(logistic_model)##
## Call:
## glm(formula = outcome_los ~ EDduration + treat_control + MOI,
## family = binomial(link = "logit"), data = relevant_data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 20.05168 1488.11864 0.013 0.9892
## EDduration -0.04082 0.01949 -2.094 0.0362 *
## treat_controltreatment -0.57581 0.47840 -1.204 0.2287
## MOIburns -15.85292 1488.11882 -0.011 0.9915
## MOIcrushing -0.26111 6690.23997 0.000 1.0000
## MOIcut -0.28276 3564.59707 0.000 0.9999
## MOIfall -16.07005 1488.11863 -0.011 0.9914
## MOIgunshot -18.72863 1488.11877 -0.013 0.9900
## MOIpenetrating inj -0.26111 6690.23997 0.000 1.0000
## MOIrti -16.78259 1488.11848 -0.011 0.9910
## MOIrti, motorbike -0.66928 4846.32890 0.000 0.9999
## MOIrti,fall -0.26111 6690.23997 0.000 1.0000
## MOIstab -17.34479 1488.11854 -0.012 0.9907
## MOIstab, assault -0.66928 6690.23999 0.000 0.9999
## MOIstab/ cut -0.74653 6690.23999 0.000 0.9999
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 147.38 on 219 degrees of freedom
## Residual deviance: 131.03 on 205 degrees of freedom
## AIC: 161.03
##
## Number of Fisher Scoring iterations: 17# 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/Afg-iMMAP-75/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\Afg-iMMAP-75\AppData\Local\Temp\Rtmpm0AkZV\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/Afg-iMMAP-75/Downloads/nick/hcp_data.csv")
head(hcp_data)## intervention_control county age gender profession speciality experience
## 1 intervention kakamega 33 female rn none 6
## 2 intervention kakamega 25 female rco none 1
## 3 intervention kakamega 26 male bsn none 4
## 4 intervention kakamega 26 female rn none 3
## 5 intervention kakamega 19 male rn midwife 1
## 6 intervention kakamega 23 male rn midwife 1
## deployment cont_education lastattended SOP ccollar_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 X
## 1 correctly performed correctly performed correctly performed NA
## 2 correctly performed correctly performed correctly performed NA
## 3 correctly performed correctly performed correctly performed NA
## 4 correctly performed correctly performed correctly performed NA
## 5 correctly performed correctly performed correctly performed NA
## 6 correctly performed correctly performed correctly performed NA# Get the variable names (column names)
variable_names <- names(hcp_data)
# Print the variable names
print(variable_names)## [1] "intervention_control" "county" "age"
## [4] "gender" "profession" "speciality"
## [7] "experience" "deployment" "cont_education"
## [10] "lastattended" "SOP" "ccollar_applied"
## [13] "airway_insert" "suction" "ETT"
## [16] "oxygensupp" "ventilator" "chesttube"
## [19] "directpressure" "tourniquet" "splint"
## [22] "pelvic_binder" "IV_Access" "IVF"
## [25] "cardiac_monitoring" "resp_monitor" "neurological"
## [28] "X"# Load necessary libraries
library(dplyr)
library(ggplot2)
library(tidyr)
# View the structure of the dataset
str(hcp_data)## 'data.frame': 60 obs. of 28 variables:
## $ intervention_control: chr "intervention" "intervention" "intervention" "intervention" ...
## $ county : chr "kakamega" "kakamega" "kakamega" "kakamega" ...
## $ age : int 33 25 26 26 19 23 40 26 26 40 ...
## $ gender : chr "female" "female" "male" "female" ...
## $ profession : chr "rn" "rco" "bsn" "rn" ...
## $ speciality : chr "none" "none" "none" "none" ...
## $ experience : int 6 1 4 3 1 1 8 2 1 1 ...
## $ deployment : int 6 1 2 1 1 1 2 1 1 1 ...
## $ cont_education : chr "bls" "none" "bls" "bls" ...
## $ lastattended : int 4 NA 1 1 NA NA 1 1 1 1 ...
## $ SOP : chr "yes" "yes" "yes" "yes" ...
## $ ccollar_applied : chr "correctly performed" "correctly performed" "correctly performed" "correctly performed" ...
## $ airway_insert : chr "correctly performed" "correctly performed" "correctly performed" "correctly performed" ...
## $ suction : chr "correctly performed" "correctly performed" "correctly performed" "correctly performed" ...
## $ ETT : chr "not performed" "correctly performed" "correctly performed" "not performed" ...
## $ oxygensupp : chr "correctly performed" "correctly performed" "correctly performed" "correctly performed" ...
## $ ventilator : chr "incorrectly performed" "incorrectly performed" "incorrectly performed" "incorrectly performed" ...
## $ chesttube : chr "incorrectly performed" "correctly performed" "not performed" "not performed" ...
## $ directpressure : chr "correctly performed" "correctly performed" "correctly performed" "correctly performed" ...
## $ tourniquet : chr "incorrectly performed" "incorrectly performed" "incorrectly performed" "incorrectly performed" ...
## $ splint : chr "correctly performed" "correctly performed" "correctly performed" "correctly performed" ...
## $ pelvic_binder : chr "incorrectly performed" "incorrectly performed" "not performed" "incorrectly performed" ...
## $ IV_Access : chr "correctly performed" "correctly performed" "correctly performed" "correctly performed" ...
## $ IVF : chr "correctly performed" "correctly performed" "correctly performed" "correctly performed" ...
## $ cardiac_monitoring : chr "correctly performed" "correctly performed" "correctly performed" "correctly performed" ...
## $ resp_monitor : chr "correctly performed" "correctly performed" "correctly performed" "correctly performed" ...
## $ neurological : chr "correctly performed" "correctly performed" "correctly performed" "correctly performed" ...
## $ X : logi NA NA NA NA NA NA ...# Summary statistics for intervention and control groups
summary_intervention <- hcp_data %>%
filter(intervention_control == "intervention") %>%
summary()
summary_control <- hcp_data %>%
filter(intervention_control == "control") %>%
summary()
print("Summary for Intervention Group:")## [1] "Summary for Intervention Group:"print(summary_intervention)## intervention_control county age gender
## Length:30 Length:30 Min. :19.00 Length:30
## Class :character Class :character 1st Qu.:25.25 Class :character
## Mode :character Mode :character Median :29.50 Mode :character
## Mean :30.43
## 3rd Qu.:35.00
## Max. :43.00
##
## profession speciality experience deployment
## Length:30 Length:30 Min. : 1.000 Min. : 1.000
## Class :character Class :character 1st Qu.: 1.250 1st Qu.: 1.000
## Mode :character Mode :character Median : 3.500 Median : 1.000
## Mean : 4.467 Mean : 1.967
## 3rd Qu.: 6.000 3rd Qu.: 2.000
## Max. :15.000 Max. :10.000
##
## cont_education lastattended SOP ccollar_applied
## Length:30 Min. :1.000 Length:30 Length:30
## Class :character 1st Qu.:1.000 Class :character Class :character
## Mode :character Median :1.000 Mode :character Mode :character
## Mean :1.385
## 3rd Qu.:1.000
## Max. :4.000
## NA's :17
## airway_insert suction ETT oxygensupp
## Length:30 Length:30 Length:30 Length:30
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## ventilator chesttube directpressure tourniquet
## Length:30 Length:30 Length:30 Length:30
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## splint pelvic_binder IV_Access IVF
## Length:30 Length:30 Length:30 Length:30
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## cardiac_monitoring resp_monitor neurological X
## Length:30 Length:30 Length:30 Mode:logical
## Class :character Class :character Class :character NA's:30
## Mode :character Mode :character Mode :character
##
##
##
## print("Summary for Control Group:")## [1] "Summary for Control Group:"print(summary_control)## intervention_control county age gender
## Length:30 Length:30 Min. :27.00 Length:30
## Class :character Class :character 1st Qu.:30.00 Class :character
## Mode :character Mode :character Median :34.50 Mode :character
## Mean :34.27
## 3rd Qu.:37.50
## Max. :44.00
##
## profession speciality experience deployment
## Length:30 Length:30 Min. : 3.000 Min. :1.000
## Class :character Class :character 1st Qu.: 4.000 1st Qu.:1.000
## Mode :character Mode :character Median : 6.500 Median :2.000
## Mean : 6.933 Mean :2.133
## 3rd Qu.:10.000 3rd Qu.:3.000
## Max. :14.000 Max. :5.000
##
## cont_education lastattended SOP ccollar_applied
## Length:30 Min. :1 Length:30 Length:30
## Class :character 1st Qu.:1 Class :character Class :character
## Mode :character Median :1 Mode :character Mode :character
## Mean :1
## 3rd Qu.:1
## Max. :1
## NA's :28
## airway_insert suction ETT oxygensupp
## Length:30 Length:30 Length:30 Length:30
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## ventilator chesttube directpressure tourniquet
## Length:30 Length:30 Length:30 Length:30
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## splint pelvic_binder IV_Access IVF
## Length:30 Length:30 Length:30 Length:30
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## cardiac_monitoring resp_monitor neurological X
## Length:30 Length:30 Length:30 Mode:logical
## Class :character Class :character Class :character NA's:30
## Mode :character Mode :character Mode :character
##
##
##
## # Frequency tables for categorical variables
cat_vars <- c("gender", "profession", "speciality", "SOP",
"ccollar_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(intervention_control, !!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
## intervention_control 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
## intervention_control 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
## intervention_control 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
## intervention_control SOP count
## <chr> <chr> <int>
## 1 control no 30
## 2 intervention no 19
## 3 intervention yes 11
## [1] "Frequency table for ccollar_applied"
## # A tibble: 5 × 3
## intervention_control ccollar_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
## intervention_control 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
## intervention_control 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
## intervention_control 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
## intervention_control oxygensupp count
## <chr> <chr> <int>
## 1 control correctly performed 30
## 2 intervention correctly performed 30
## [1] "Frequency table for ventilator"
## # A tibble: 6 × 3
## intervention_control 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
## intervention_control 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
## intervention_control directpressure count
## <chr> <chr> <int>
## 1 control correctly performed 30
## 2 intervention correctly performed 30
## [1] "Frequency table for tourniquet"
## # A tibble: 6 × 3
## intervention_control 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
## intervention_control 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
## intervention_control 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
## intervention_control 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
## intervention_control 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
## intervention_control 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
## intervention_control 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
## intervention_control 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 = intervention_control)) +
geom_bar(position = "dodge") +
labs(title = paste("Distribution of", var, "by Group"),
x = var, y = "Count") +
theme_minimal()
print(plot)
}
# Chi-square test for each categorical variable
for (var in cat_vars) {
chi_test <- chisq.test(table(hcp_data[[var]], hcp_data$intervention_control))
print(paste("Chi-square test for", var))
print(chi_test)
}## [1] "Chi-square test for gender"
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: table(hcp_data[[var]], hcp_data$intervention_control)
## X-squared = 0.066667, df = 1, p-value = 0.7963
##
## [1] "Chi-square test for profession"
##
## Pearson's Chi-squared test
##
## data: table(hcp_data[[var]], hcp_data$intervention_control)
## X-squared = 4.0185, df = 4, p-value = 0.4035
##
## [1] "Chi-square test for speciality"
##
## Pearson's Chi-squared test
##
## data: table(hcp_data[[var]], hcp_data$intervention_control)
## X-squared = 5.4545, df = 4, p-value = 0.2438
##
## [1] "Chi-square test for SOP"
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: table(hcp_data[[var]], hcp_data$intervention_control)
## X-squared = 11.132, df = 1, p-value = 0.0008486
##
## [1] "Chi-square test for ccollar_applied"
##
## Pearson's Chi-squared test
##
## data: table(hcp_data[[var]], hcp_data$intervention_control)
## X-squared = 16.923, df = 2, p-value = 0.0002114
##
## [1] "Chi-square test for airway_insert"
##
## Pearson's Chi-squared test
##
## data: table(hcp_data[[var]], hcp_data$intervention_control)
## X-squared = 9.5772, df = 2, p-value = 0.008324
##
## [1] "Chi-square test for suction"
##
## Pearson's Chi-squared test
##
## data: table(hcp_data[[var]], hcp_data$intervention_control)
## X-squared = 0, df = 1, p-value = 1
##
## [1] "Chi-square test for ETT"
##
## Pearson's Chi-squared test
##
## data: table(hcp_data[[var]], hcp_data$intervention_control)
## X-squared = 19.111, df = 2, p-value = 7.081e-05
##
## [1] "Chi-square test for oxygensupp"
##
## Chi-squared test for given probabilities
##
## data: table(hcp_data[[var]], hcp_data$intervention_control)
## X-squared = 0, df = 1, p-value = 1
##
## [1] "Chi-square test for ventilator"
##
## Pearson's Chi-squared test
##
## data: table(hcp_data[[var]], hcp_data$intervention_control)
## X-squared = 11.156, df = 2, p-value = 0.00378
##
## [1] "Chi-square test for chesttube"
##
## Pearson's Chi-squared test
##
## data: table(hcp_data[[var]], hcp_data$intervention_control)
## X-squared = 1.004, df = 2, p-value = 0.6053
##
## [1] "Chi-square test for directpressure"
##
## Chi-squared test for given probabilities
##
## data: table(hcp_data[[var]], hcp_data$intervention_control)
## X-squared = 0, df = 1, p-value = 1
##
## [1] "Chi-square test for tourniquet"
##
## Pearson's Chi-squared test
##
## data: table(hcp_data[[var]], hcp_data$intervention_control)
## X-squared = 0.22222, df = 2, p-value = 0.8948
##
## [1] "Chi-square test for splint"
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: table(hcp_data[[var]], hcp_data$intervention_control)
## X-squared = 1.6968, df = 1, p-value = 0.1927
##
## [1] "Chi-square test for pelvic_binder"
##
## Pearson's Chi-squared test
##
## data: table(hcp_data[[var]], hcp_data$intervention_control)
## X-squared = 5.1173, df = 2, p-value = 0.07741
##
## [1] "Chi-square test for IV_Access"
##
## Chi-squared test for given probabilities
##
## data: table(hcp_data[[var]], hcp_data$intervention_control)
## X-squared = 0, df = 1, p-value = 1
##
## [1] "Chi-square test for IVF"
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: table(hcp_data[[var]], hcp_data$intervention_control)
## X-squared = 1.7811, df = 1, p-value = 0.182
##
## [1] "Chi-square test for cardiac_monitoring"
##
## Chi-squared test for given probabilities
##
## data: table(hcp_data[[var]], hcp_data$intervention_control)
## X-squared = 0, df = 1, p-value = 1
##
## [1] "Chi-square test for resp_monitor"
##
## Chi-squared test for given probabilities
##
## data: table(hcp_data[[var]], hcp_data$intervention_control)
## X-squared = 0, df = 1, p-value = 1
##
## [1] "Chi-square test for neurological"
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: table(hcp_data[[var]], hcp_data$intervention_control)
## X-squared = 22.937, df = 1, p-value = 1.674e-06# T-test for continuous variables
cont_vars <- c("age", "experience")
for (var in cont_vars) {
t_test <- t.test(hcp_data[[var]] ~ hcp_data$intervention_control)
print(paste("T-test for", var))
print(t_test)
}## [1] "T-test for age"
##
## Welch Two Sample t-test
##
## data: hcp_data[[var]] by hcp_data$intervention_control
## t = 2.5544, df = 51.242, p-value = 0.01365
## alternative hypothesis: true difference in means between group control and group intervention is not equal to 0
## 95 percent confidence interval:
## 0.8209704 6.8456962
## sample estimates:
## mean in group control mean in group intervention
## 34.26667 30.43333
##
## [1] "T-test for experience"
##
## Welch Two Sample t-test
##
## data: hcp_data[[var]] by hcp_data$intervention_control
## t = 2.8292, df = 57.212, p-value = 0.006425
## alternative hypothesis: true difference in means between group control and group intervention is not equal to 0
## 95 percent confidence interval:
## 0.7209221 4.2124112
## sample estimates:
## mean in group control mean in group intervention
## 6.933333 4.466667# 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 ~ intervention_control + age + county + gender + profession + speciality, data = hcp_data, family = binomial)
# Summarize the model
summary(logit_model)##
## Call:
## glm(formula = SOP_dummy ~ intervention_control + 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
## intervention_controlintervention -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# Load necessary libraries
library(dplyr)
library(caret)
library(DescTools) # For effect size calculation (Cramér's V)##
## Attaching package: 'DescTools'## The following objects are masked from 'package:caret':
##
## MAE, RMSE## The following objects are masked from 'package:psych':
##
## AUC, ICC, SDlibrary(stargazer) # For clean table output
# Assuming hcp_data is already loaded and dummified
# List of variables to analyze
variables_of_interest <- c("ccollar_applied", "airway_insert", "suction", "ETT", "oxygensupp",
"ventilator", "chesttube", "directpressure", "tourniquet", "splint",
"pelvic_binder", "IV_Access", "IVF", "cardiac_monitoring",
"resp_monitor", "neurological")
# Function to apply Chi-Square or Fisher's test for each variable
test_statistics <- sapply(variables_of_interest, function(var) {
# Create a contingency table for the variable against intervention_control
contingency_table <- table(hcp_data[[var]], hcp_data$intervention_control)
# Check if Chi-square test is applicable (i.e., no expected frequency < 5)
if (all(dim(contingency_table) > 1)) {
chi_sq_test <- chisq.test(contingency_table)
# If any expected frequency is less than 5, switch to Fisher's test
if (any(chi_sq_test$expected < 5)) {
fisher_test <- fisher.test(contingency_table)
return(c(test = "Fisher's Exact Test", p_value = fisher_test$p.value, statistic = fisher_test$estimate))
} else {
return(c(test = "Chi-square Test", p_value = chi_sq_test$p.value, statistic = chi_sq_test$statistic))
}
} else {
# If no valid data, return NA
return(c(test = "No valid data", p_value = NA, statistic = NA))
}
}, simplify = FALSE)
# Convert results into a data frame and bind them together
results_df <- bind_rows(lapply(test_statistics, function(x) as.data.frame(t(x), stringsAsFactors = FALSE)))
# Add variable names as a column
results_df$variable <- variables_of_interest
# Print the results in a cleaner format using stargazer
stargazer(results_df,
type = "text",
summary = FALSE,
digits = 3,
header = FALSE,
rownames = FALSE,
title = "Chi-square / Fisher's Test Results for Variables between Intervention and Control Groups")##
## Chi-square / Fisher's Test Results for Variables between Intervention and Control Groups
## ==============================================================================================================
## test p_value statistic.X-squared statistic.odds ratio statistic variable
## --------------------------------------------------------------------------------------------------------------
## Fisher's Exact Test 0.000160221245787189 ccollar_applied
## Chi-square Test 0.00832405193975244 9.57721226142279 airway_insert
## Fisher's Exact Test 1 1 suction
## Fisher's Exact Test 1.83880218106987e-05 ETT
## No valid data oxygensupp
## Chi-square Test 0.00378003512138551 11.156043956044 ventilator
## Fisher's Exact Test 0.586792883071943 chesttube
## No valid data directpressure
## Fisher's Exact Test 0.999999999999918 tourniquet
## Chi-square Test 0.192702728161021 1.69683257918552 splint
## Fisher's Exact Test 0.0866849171471914 pelvic_binder
## No valid data IV_Access
## Chi-square Test 0.182017247136709 1.78107606679035 IVF
## No valid data cardiac_monitoring
## No valid data resp_monitor
## Chi-square Test 1.67441317666643e-06 22.9365079365079 neurological
## --------------------------------------------------------------------------------------------------------------# Assuming test_statistics from the previous code is generated
# Convert results into a data frame
results_df <- bind_rows(lapply(test_statistics, function(x) as.data.frame(t(x), stringsAsFactors = FALSE)))
# Add column names for clarity
colnames(results_df) <- c("Test Type", "P-value", "Statistic")
# Add variable names as a column
results_df$Variable <- variables_of_interest
# Reorder columns for readability
results_df <- results_df[, c("Variable", "Test Type", "P-value", "Statistic")]
# Use stargazer to create a neat, readable table in text format
stargazer(results_df,
type = "text",
summary = FALSE,
digits = 3,
header = FALSE,
rownames = FALSE,
title = "Chi-Square / Fisher's Test Results: Intervention vs Control")##
## Chi-Square / Fisher's Test Results: Intervention vs Control
## ============================================================================
## Variable Test Type P-value Statistic
## ----------------------------------------------------------------------------
## ccollar_applied Fisher's Exact Test 0.000160221245787189
## airway_insert Chi-square Test 0.00832405193975244 9.57721226142279
## suction Fisher's Exact Test 1
## ETT Fisher's Exact Test 1.83880218106987e-05
## oxygensupp No valid data
## ventilator Chi-square Test 0.00378003512138551 11.156043956044
## chesttube Fisher's Exact Test 0.586792883071943
## directpressure No valid data
## tourniquet Fisher's Exact Test 0.999999999999918
## splint Chi-square Test 0.192702728161021 1.69683257918552
## pelvic_binder Fisher's Exact Test 0.0866849171471914
## IV_Access No valid data
## IVF Chi-square Test 0.182017247136709 1.78107606679035
## cardiac_monitoring No valid data
## resp_monitor No valid data
## neurological Chi-square Test 1.67441317666643e-06 22.9365079365079
## ----------------------------------------------------------------------------# Load necessary libraries
library(dplyr)
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:DescTools':
##
## Recode## The following object is masked from 'package:likert':
##
## recode## The following object is masked from 'package:purrr':
##
## some## The following object is masked from 'package:psych':
##
## logit## The following object is masked from 'package:dplyr':
##
## recodelibrary(MASS)
library(stargazer)
# List of variables to analyze
variables_of_interest <- c("ccollar_applied", "airway_insert", "suction", "ETT", "oxygensupp",
"ventilator", "chesttube", "directpressure", "tourniquet", "splint",
"pelvic_binder", "IV_Access", "IVF", "cardiac_monitoring",
"resp_monitor", "neurological")
# Create a data frame to store the results
results_df <- data.frame(
variable = character(),
test = character(),
p_value = numeric(),
statistic = numeric(),
stringsAsFactors = FALSE
)
# Function to perform Chi-Square Test for each variable and capture results
for (var in variables_of_interest) {
# Ensure the variable is treated as a factor
hcp_data[[var]] <- factor(hcp_data[[var]], levels = c("not performed", "correctly performed", "incorrectly performed"))
# Create a contingency table for the Chi-square test
contingency_table <- table(hcp_data[[var]], hcp_data$intervention_control)
# Perform the Chi-Square test of independence
chi_sq_test <- chisq.test(contingency_table)
# If the expected counts are low, apply Fisher's Exact Test instead
if (any(chi_sq_test$expected < 5)) {
fisher_test <- fisher.test(contingency_table)
results_df <- rbind(results_df, data.frame(
variable = var,
test = "Fisher's Exact Test",
p_value = fisher_test$p.value,
statistic = NA # Fisher's test does not return a statistic like Chi-square
))
} else {
results_df <- rbind(results_df, data.frame(
variable = var,
test = "Chi-Square Test",
p_value = chi_sq_test$p.value,
statistic = chi_sq_test$statistic
))
}
}
# Use stargazer to generate the table in long format
stargazer(results_df, type = "text", summary = FALSE, rownames = FALSE,
column.labels = c("Variable", "Test", "P-Value", "Statistic"),
digits = 3)##
## ========================================================
## variable test p_value statistic
## --------------------------------------------------------
## ccollar_applied Fisher's Exact Test 0.011
## airway_insert Chi-Square Test 0.008 9.577
## suction Fisher's Exact Test 1
## ETT Fisher's Exact Test 0.00002
## oxygensupp Fisher's Exact Test 1
## ventilator Chi-Square Test 0.004 11.156
## chesttube Fisher's Exact Test 0.587
## directpressure Fisher's Exact Test 1
## tourniquet Fisher's Exact Test 1
## splint Fisher's Exact Test 0.192
## pelvic_binder Fisher's Exact Test 0.087
## IV_Access Fisher's Exact Test 1
## IVF Fisher's Exact Test 0.181
## cardiac_monitoring Fisher's Exact Test 1
## resp_monitor Fisher's Exact Test 1
## neurological Fisher's Exact Test 0.00000
## --------------------------------------------------------#Skills in performance of procedures
# Load necessary library
library(dplyr)
# List of variables to analyze
variables_of_interest <- c("ccollar_applied", "airway_insert", "suction", "ETT", "oxygensupp",
"ventilator", "chesttube", "directpressure", "tourniquet", "splint",
"pelvic_binder", "IV_Access", "IVF", "cardiac_monitoring",
"resp_monitor", "neurological")
# Create an empty data frame to store the results
table_df <- data.frame(
variable = character(),
intervention_control = character(),
correctly_performed = integer(),
incorrectly_performed = integer(),
not_performed = integer(),
Total = integer(),
stringsAsFactors = FALSE
)
# Loop through each variable to create contingency tables
for (var in variables_of_interest) {
# Ensure the variable is treated as a factor
hcp_data[[var]] <- factor(hcp_data[[var]], levels = c("not performed", "correctly performed", "incorrectly performed"))
# Summarize counts for each intervention group
counts <- hcp_data %>%
group_by(intervention_control, !!sym(var)) %>%
summarise(count = n(), .groups = "drop") %>%
tidyr::pivot_wider(names_from = !!sym(var), values_from = count, values_fill = 0)
# Add missing columns if necessary
for (col in c("correctly performed", "incorrectly performed", "not performed")) {
if (!(col %in% colnames(counts))) {
counts[[col]] <- 0
}
}
# Compute total
counts$Total <- rowSums(counts[, c("correctly performed", "incorrectly performed", "not performed")], na.rm = TRUE)
# Add variable name
counts$variable <- var
# Append to the results data frame
table_df <- rbind(table_df, counts[, c("variable", "intervention_control", "correctly performed", "incorrectly performed", "not performed", "Total")])
}
# Print the final table
table_df## # A tibble: 32 × 6
## variable intervention_control `correctly performed` incorrectly performe…¹
## <chr> <chr> <int> <dbl>
## 1 ccollar_ap… control 6 16
## 2 ccollar_ap… intervention 20 10
## 3 airway_ins… control 4 17
## 4 airway_ins… intervention 15 11
## 5 suction control 27 3
## 6 suction intervention 27 3
## 7 ETT control 0 4
## 8 ETT intervention 8 12
## 9 oxygensupp control 30 0
## 10 oxygensupp intervention 30 0
## # ℹ 22 more rows
## # ℹ abbreviated name: ¹`incorrectly performed`
## # ℹ 2 more variables: `not performed` <dbl>, Total <dbl># Load necessary libraries
library(dplyr)
library(tidyr)
# List of variables to analyze
variables_of_interest <- c("ccollar_applied", "airway_insert", "suction", "ETT", "oxygensupp",
"ventilator", "chesttube", "directpressure", "tourniquet", "splint",
"pelvic_binder", "IV_Access", "IVF", "cardiac_monitoring",
"resp_monitor", "neurological")
# Create an empty data frame to store the results
table_df <- data.frame(
variable = character(),
profession = character(),
intervention_control = character(),
correctly_performed = integer(),
incorrectly_performed = integer(),
not_performed = integer(),
Total = integer(),
stringsAsFactors = FALSE
)
# Loop through each variable to create contingency tables
for (var in variables_of_interest) {
# Ensure the variable is treated as a factor
hcp_data[[var]] <- factor(hcp_data[[var]], levels = c("not performed", "correctly performed", "incorrectly performed"))
# Summarize counts for each profession and intervention group
counts <- hcp_data %>%
group_by(profession, intervention_control, !!sym(var)) %>%
summarise(count = n(), .groups = "drop") %>%
pivot_wider(names_from = !!sym(var), values_from = count, values_fill = list(count = 0))
# Add missing columns if necessary
for (col in c("correctly performed", "incorrectly performed", "not performed")) {
if (!(col %in% colnames(counts))) {
counts[[col]] <- 0
}
}
# Compute total
counts$Total <- rowSums(counts[, c("correctly performed", "incorrectly performed", "not performed")], na.rm = TRUE)
# Add variable name
counts$variable <- var
# Append to the results data frame
table_df <- rbind(table_df, counts[, c("variable", "profession", "intervention_control", "correctly performed", "incorrectly performed", "not performed", "Total")])
}
# Print the final table
table_df## # A tibble: 144 × 7
## variable profession intervention_control `correctly performed`
## <chr> <chr> <chr> <int>
## 1 ccollar_applied bsn control 0
## 2 ccollar_applied bsn intervention 2
## 3 ccollar_applied mo control 0
## 4 ccollar_applied mo intervention 0
## 5 ccollar_applied paramedic intervention 1
## 6 ccollar_applied rco control 2
## 7 ccollar_applied rco intervention 2
## 8 ccollar_applied rn control 4
## 9 ccollar_applied rn intervention 15
## 10 airway_insert bsn control 0
## # ℹ 134 more rows
## # ℹ 3 more variables: `incorrectly performed` <dbl>, `not performed` <dbl>,
## # Total <dbl># Load necessary library
library(dplyr)
# List of variables to analyze
variables_of_interest <- c("ccollar_applied", "airway_insert", "suction", "ETT", "oxygensupp",
"ventilator", "chesttube", "directpressure", "tourniquet", "splint",
"pelvic_binder", "IV_Access", "IVF", "cardiac_monitoring",
"resp_monitor", "neurological")
# Create an empty data frame to store the results
table_df <- data.frame(
variable = character(),
profession = character(),
intervention_control = character(),
performed = integer(),
not_performed_or_incorrect = integer(),
Total = integer(),
stringsAsFactors = FALSE
)
# Loop through each variable to create contingency tables
for (var in variables_of_interest) {
# Ensure the variable is treated as a factor
hcp_data[[var]] <- factor(hcp_data[[var]], levels = c("not performed", "correctly performed", "incorrectly performed"))
# Summarize counts for each profession and intervention group
counts <- hcp_data %>%
group_by(profession, intervention_control) %>%
summarise(
performed = sum(!!sym(var) == "correctly performed", na.rm = TRUE),
not_performed_or_incorrect = sum(!!sym(var) %in% c("not performed", "incorrectly performed"), na.rm = TRUE),
Total = n(),
.groups = "drop"
)
# Add variable name
counts$variable <- var
# Append to the results data frame
table_df <- rbind(table_df, counts[, c("variable", "profession", "intervention_control", "performed", "not_performed_or_incorrect", "Total")])
}
# Print the final table
table_df## # A tibble: 144 × 6
## variable profession intervention_control performed not_performed_or_inc…¹
## <chr> <chr> <chr> <int> <int>
## 1 ccollar_app… bsn control 0 2
## 2 ccollar_app… bsn intervention 2 1
## 3 ccollar_app… mo control 0 3
## 4 ccollar_app… mo intervention 0 1
## 5 ccollar_app… paramedic intervention 1 0
## 6 ccollar_app… rco control 2 3
## 7 ccollar_app… rco intervention 2 2
## 8 ccollar_app… rn control 4 8
## 9 ccollar_app… rn intervention 15 6
## 10 airway_inse… bsn control 0 4
## # ℹ 134 more rows
## # ℹ abbreviated name: ¹not_performed_or_incorrect
## # ℹ 1 more variable: Total <int># Save the final table as a CSV file
write.csv(table_df, "hcp_intervention_summary.csv", row.names = FALSE)# Load necessary library
library(dplyr)
# List of variables to analyze
variables_of_interest <- c("ccollar_applied", "airway_insert", "suction", "ETT", "oxygensupp",
"ventilator", "chesttube", "directpressure", "tourniquet", "splint",
"pelvic_binder", "IV_Access", "IVF", "cardiac_monitoring",
"resp_monitor", "neurological")
# Create an empty data frame to store the results
table_df <- data.frame(
variable = character(),
cont_education = character(),
intervention_control = character(),
performed = integer(),
not_performed_or_incorrect = integer(),
Total = integer(),
stringsAsFactors = FALSE
)
# Loop through each variable to create contingency tables
for (var in variables_of_interest) {
# Ensure the variable is treated as a factor
hcp_data[[var]] <- factor(hcp_data[[var]], levels = c("not performed", "correctly performed", "incorrectly performed"))
# Summarize counts for each profession and intervention group
counts <- hcp_data %>%
group_by(cont_education, intervention_control) %>%
summarise(
performed = sum(!!sym(var) == "correctly performed", na.rm = TRUE),
not_performed_or_incorrect = sum(!!sym(var) %in% c("not performed", "incorrectly performed"), na.rm = TRUE),
Total = n(),
.groups = "drop"
)
# Add variable name
counts$variable <- var
# Append to the results data frame
table_df <- rbind(table_df, counts[, c("variable", "cont_education", "intervention_control", "performed", "not_performed_or_incorrect", "Total")])
}
# Print the final table
table_df## # A tibble: 96 × 6
## variable cont_education intervention_control performed not_performed_or_inc…¹
## <chr> <chr> <chr> <int> <int>
## 1 ccollar… bls control 1 1
## 2 ccollar… bls intervention 8 2
## 3 ccollar… bls, acls intervention 1 0
## 4 ccollar… bls, acls, at… intervention 2 0
## 5 ccollar… none control 5 15
## 6 ccollar… none intervention 9 8
## 7 airway_… bls control 1 1
## 8 airway_… bls intervention 7 3
## 9 airway_… bls, acls intervention 1 0
## 10 airway_… bls, acls, at… intervention 2 0
## # ℹ 86 more rows
## # ℹ abbreviated name: ¹not_performed_or_incorrect
## # ℹ 1 more variable: Total <int># Save the final table as a CSV file
write.csv(table_df, "hcp_intervention_cont_education.csv", row.names = FALSE)# Load necessary libraries
library(dplyr)
# Create a dummy variable for SOP (yes = 1, no = 0)
hcp_data$SOP_dummy <- ifelse(hcp_data$SOP == "yes", 1, 0)
# Summarize counts for correctly performed, incorrectly performed, and not performed
summary_table <- hcp_data %>%
group_by(SOP_dummy, intervention_control) %>%
summarise(
correctly_performed = sum(hcp_data == "correctly performed", na.rm = TRUE),
incorrectly_performed = sum(hcp_data == "incorrectly performed", na.rm = TRUE),
not_performed = sum(hcp_data == "not performed", na.rm = TRUE),
Total = n(),
.groups = "drop"
) %>%
mutate(SOP = ifelse(SOP_dummy == 1, "Yes", "No")) %>%
select(SOP, intervention_control, correctly_performed, incorrectly_performed, not_performed, Total)
# Save the final table as a CSV file
write.csv(summary_table, "hcp_intervention_SOP_summary.csv", row.names = FALSE)
# Print the final table
summary_table## # A tibble: 3 × 6
## SOP intervention_control correctly_performed incorrectly_performed
## <chr> <chr> <int> <int>
## 1 No control 567 265
## 2 No intervention 567 265
## 3 Yes intervention 567 265
## # ℹ 2 more variables: not_performed <int>, Total <int># Save the final table as a CSV file
write.csv(summary_table, "hcp_intervention_SOP_summary.csv", row.names = FALSE)