marty temporal bone

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  df2$dlptot and df3$dlptot
## W = 13894, p-value = 2.848e-11
## alternative hypothesis: true location shift is not equal to 0

#check distribution of DLP data
x <- df$dlptot
summary(x)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0   239.5   459.1   593.3   778.0  2863.5

ggqqplot(x)

shapiro.test(x)

## 
##  Shapiro-Wilk normality test
## 
## data:  x
## W = 0.85003, p-value < 2.2e-16

#DLP is hella skewed so need to use non parametric test

#summarize data
tib <- df %>%
  group_by(tbf) %>%
            summarise(n_patients=n(),mean_DLP = mean(dlptot),
            n_patients=n(),median_DLP = median(dlptot),
            std_error = sd(dlptot) / sqrt(n_patients))

tib

## # A tibble: 2 × 5
##     tbf n_patients mean_DLP median_DLP std_error
##   <dbl>      <int>    <dbl>      <dbl>     <dbl>
## 1     0        364     517.       414.      23.4
## 2     1         48    1172.      1016.     104.

table(df$mechanism_cat)

## 
##                            ATV Bike, Scooter, Skateboard, etc 
##                             19                             35 
##                   Blunt Injury         Car against Pedestrian 
##                             38                             15 
##                           Fall        Motor Vehicle Collision 
##                            198                             65 
##                          Other 
##                             42

#visual summariy boxplots
ggboxplot(df, x = "tbf", y = "dlptot", lwd=.6, show.legend = F,
          color = "tbf", palette = c("firebrick2","royalblue"),
          ylab = "DLP (mGy*cm)", xlab = "TBF",width=.7,bxp.errorbar.width = 4)+geom_point(size=2,alpha=.5,position = position_jitterdodge(jitter.width = .6), aes(color = factor(tbf)))+ scale_x_discrete(name ="", labels=c("0" = "Other \nfracture\ntype", "1" = "Temporal \nbone\nFracture"))+ggtitle("Total Radiation Exposure \n(in DLP) by Fracture Type") +   scale_y_continuous(breaks = round(seq(min(df$dlptot), max(df$dlptot), by = 600),1))+theme_gray(base_size = 14)+theme(legend.position = "none")+gghisto

ggplot(df, aes(x = tbf, y = dlptot), add = c("mean_se")) + 
  geom_point(size=5,shape=16, position = position_jitterdodge(jitter.width = .12),aes(color=mechanism_cat),alpha=.4, palette = c("firebrick2","royalblue")) +gghisto+
theme(plot.caption = element_text(hjust = 0))+ ggtitle("Total Radiation Exposure \n(in DLP) by Fracture Type") +
stat_summary(aes(tbf, dlptot), data = df, fun = "mean", geom = "crossbar", size=0.2, width=0.8, color="brown")+ scale_x_discrete(name ="", labels=c("0" = "Other \nfracture\ntype", "1" = "Temporal \nbone\nFracture"))+   scale_y_continuous(breaks = round(seq(min(df$dlptot), max(df$dlptot), by = 400),1))+
  ylab("DLP (mGy*cm)")+ xlab("TBF")

mech_labs <- c("ATV","Bike, \nScooter, \nSkateboard, \netc","Blunt \nInjury","Car v \nPedestrian","Fall","Motor \nVehicle \nCollision", "Other")

ggplot(df) + geom_bar(color="black",aes(mechanism_cat, dlptot, fill = as.factor(tbf_cat)), 
           position = "dodge", stat = "summary", fun = "mean")+gghisto+ggtitle("Total Radiation by Injury Mechanism \nand Fracture Type")+ scale_fill_manual(name = "Fracture Type",values=c("indianred3", "royalblue1"))+theme(axis.text.x = element_text(face="bold", color="royalblue4", size=10))+ scale_x_discrete(labels= mech_labs)+ylab("DLP (mGy*cm)")+xlab("Mechanism of Injury")

ggplot(df) + geom_bar(color="black",aes(GCS_cat, dlptot, fill = as.factor(tbf_cat)), 
           position = "dodge", stat = "summary", fun = "mean")+gghisto+ggtitle("Total Radiation by Injury Mechanism \nand Fracture Type")+ scale_fill_manual(name = "Fracture Type",values=c("indianred3", "royalblue1"))+theme(axis.text.x = element_text(face="bold", color="royalblue4", size=14))+ylab("DLP (mGy*cm)")+xlab("GCS Score")

ggplot(df, aes(x = Age, y = dlptot))+ geom_point()+
  geom_smooth(method = "lm") +   scale_y_continuous(breaks = round(seq(min(df$dlptot), max(df$dlptot), by = 300),1))+   scale_x_continuous(breaks = round(seq(min(df$Age), max(df$Age), by = 1),1))+gghisto+ggtitle("Total Radiation by Age")

ggplot(df, aes(x = Age, y = dlptot)) + geom_point(aes(color = tbf_cat))+
  geom_smooth(aes(color = tbf_cat, fill = tbf_cat), method = "lm") +
  scale_color_manual(name = "Fracture Type",values = c("indianred3", "royalblue1"))+
  scale_fill_manual(name = "Fracture Type",values = c("indianred3", "royalblue1"))+   scale_y_continuous(breaks = round(seq(min(df$dlptot), max(df$dlptot), by = 300),1))+   scale_x_continuous(breaks = round(seq(min(df$Age), max(df$Age), by = 1),1))+gghisto+ggtitle("Total Radiation by Age \nand Fracture Type")

ggplot(df, aes(x = Age, y = dlptot)) + geom_point(aes(color = tbf_cat)) + 
  geom_smooth(aes(color = tbf_cat, fill = tbf_cat), method = "lm",fullrange=TRUE) +
  facet_wrap(~tbf_cat) +
  scale_color_manual(name = "Fracture Type",values = c("indianred3", "royalblue1"))+
  scale_fill_manual(name = "Fracture Type",values = c("indianred3", "royalblue1"))+
  scale_y_continuous(breaks = round(seq(min(df$dlptot), max(df$dlptot), by = 300),1))+   scale_x_continuous(breaks = round(seq(min(df$Age), max(df$Age), by = 1),1))+gghisto+ggtitle("Total Radiation by Age and Fracture Type")+
  theme_bw()

ggboxplot(df, x = "tbf", y = "Age", lwd=.6, show.legend = F,
          color = "tbf", palette = c("indianred3", "royalblue1"),
          ylab = "Age", xlab = "TBF",width=.35,bxp.errorbar.width = 4)+geom_point(size=2,alpha=.5,position = position_jitterdodge(jitter.width = 0.6),aes(color = factor(tbf)))+ scale_x_discrete(name ="", labels=c("0" = "Other \nfracture\ntype", "1" = "Temporal \nbone\nFracture"))+ggtitle("Age by Fracture Type") +   scale_y_continuous(breaks = round(seq(min(df$Age), max(df$Age), by = 2),1))+theme_gray(base_size = 14)+theme(legend.position = "none")+gghisto

#visual summaries
no_y <- list(theme(axis.ticks.y = element_blank(),
        axis.text.y = element_blank()))

ggplot(df, aes(x=dlptot)) + 
  geom_histogram(aes(y=..density..),      
                 binwidth=110,
                color="gray18",fill="white")+
  geom_density(alpha=.2, fill="indianred2",size=.5, color="gray14") +
  ggtitle("Total DLP distribution for all patients")+
  gghisto+no_y+   scale_x_continuous(breaks = round(seq(min(df$dlptot), max(df$dlptot), by = 400),1))

ggplot(df, aes(x = dlptot, fill = tbf_cat)) +geom_density(alpha = .3)+ ggtitle("Total DLP distribution by Fracture Type")+gghisto+ no_y+   scale_x_continuous(breaks = round(seq(min(df$dlptot), max(df$dlptot), by = 400),1))+ theme(legend.position = c(.8, .7))

#wilcox is non parametric ttest for your predictor vs outcome 
w <- wilcox.test(df2$dlptot, df3$dlptot, alternative = "two.sided")
w

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  df2$dlptot and df3$dlptot
## W = 13894, p-value = 2.848e-11
## alternative hypothesis: true location shift is not equal to 0

#data summary
tbl_summary(df[,c(1:4,6,8,11:13)],label = dlptot ~ "Total DLP",by=tbf_cat)%>%
  bold_labels()  %>%add_p()

Characteristic	Other Fracture Type, N = 364¹	Temp Bone Fracture, N = 48¹	p-value²
Age	1.7 (0.4, 7.6)	12.0 (8.4, 14.9)	<0.001
sex	227 (62%)	31 (65%)	0.9
Transfer	312 (86%)	32 (67%)	0.002
LOS (d)	2.0 (0.9, 3.0)	4.5 (2.0, 11.0)	<0.001
Unknown	1	0
mechanism_cat
ATV	15 (4.1%)	4 (8.3%)
Bike, Scooter, Skateboard, etc	24 (6.6%)	11 (23%)
Blunt Injury	33 (9.1%)	5 (10%)
Car against Pedestrian	8 (2.2%)	7 (15%)
Fall	192 (53%)	6 (12%)
Motor Vehicle Collision	51 (14%)	14 (29%)
Other	41 (11%)	1 (2.1%)
GCS_cat			<0.001
15	299 (82%)	23 (48%)
Below 15	65 (18%)	25 (52%)
Hemorrhage Type			<0.001
0	180 (49%)	21 (44%)
1	69 (19%)	27 (56%)
2	69 (19%)	0 (0%)
3	12 (3.3%)	0 (0%)
4	34 (9.3%)	0 (0%)
Total DLP	414 (222, 697)	1,016 (700, 1,531)	<0.001
¹ Statistics presented: median (IQR); n (%) ² Statistical tests performed: Wilcoxon rank-sum test; chi-square test of independence; Fisher's exact test

#summarize covariates

df_cov <- c(names(df[c(1:5,9,11,12,13)]))

cov.summary <- df %>%
  group_by(tbf_cat) %>%
  select(one_of(df_cov)) %>%
  summarise_all(funs(median(., na.rm = T)))

## Warning: `funs()` was deprecated in dplyr 0.8.0.
## Please use a list of either functions or lambdas: 
## 
##   # Simple named list: 
##   list(mean = mean, median = median)
## 
##   # Auto named with `tibble::lst()`: 
##   tibble::lst(mean, median)
## 
##   # Using lambdas
##   list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))

cov_before <- cov.summary

options(scipen=999)
dfx <- df[,c(1:5,9,10,12,13)]

#wilcoxon tests to see if these medians are statistically distinguishable
df_cov <- c(names(dfx))

cov_w_tests <- lapply(dfx[,df_cov], function(x) wilcox.test(x ~ dfx$tbf))

cov.pvals <- as.matrix(round(sapply(cov_w_tests, function(x) {
    p <- x$p.value
    n <- outer(rownames(p), colnames(p), paste, sep='v')
    p <- as.vector(p)
    names(p) <- n
    p
}),6))

cov.pvals

##                        [,1]
## Age.NA             0.000000
## sex.NA             0.765874
## Transfer.NA        0.000849
## LOS (d).NA         0.000000
## Mechanism.NA       0.000420
## GCS_bin.NA         0.000000
## tbf.NA             0.000000
## Hemorrhage Type.NA 0.144119
## dlptot.NA          0.000000

#propensity score estimation
#create regression model
ps1 <- glm(data=dfx,tbf ~.,family = "binomial")

summary(ps1)

## 
## Call:
## glm(formula = tbf ~ ., family = "binomial", data = dfx)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.1455  -0.3820  -0.2446  -0.1384   2.7528  
## 
## Coefficients:
##                     Estimate Std. Error z value  Pr(>|z|)    
## (Intercept)       -2.8384640  0.7072839  -4.013 0.0000599 ***
## Age                0.1238928  0.0401010   3.090   0.00200 ** 
## sex               -0.0471953  0.3943150  -0.120   0.90473    
## Transfer          -0.3787337  0.4522372  -0.837   0.40233    
## `LOS (d)`          0.0328483  0.0223900   1.467   0.14235    
## Mechanism          0.0922001  0.0974052   0.947   0.34386    
## GCS_bin           -0.8844372  0.4363961  -2.027   0.04269 *  
## `Hemorrhage Type` -0.7686272  0.2414927  -3.183   0.00146 ** 
## dlptot             0.0011137  0.0003595   3.098   0.00195 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 296.31  on 410  degrees of freedom
## Residual deviance: 201.35  on 402  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: 219.35
## 
## Number of Fisher Scoring iterations: 6

#calculate propensity score
#prs_df in the guide

scores <- data.frame(score = predict(ps1, type ="response"), tbf = ps1$model$tbf)

#examine region of common support

labs <- paste0("Fracture Type: ",c("Temporal","Other"))

scores %>%
  mutate(tbf = ifelse(tbf == 1, labs[1], labs[2])) %>%
  ggplot(aes(x=score)) +
  geom_histogram(color="white",bins = 14) +
  facet_wrap(~tbf) +
  xlab("Probability of being admitted for a Temporal Bone Fracture") +
  theme_bw()

#remove all missing values
apply(is.na(dfx), 2, which)

## $Age
## integer(0)
## 
## $sex
## integer(0)
## 
## $Transfer
## integer(0)
## 
## $`LOS (d)`
## [1] 147
## 
## $Mechanism
## integer(0)
## 
## $GCS_bin
## integer(0)
## 
## $tbf
## integer(0)
## 
## $`Hemorrhage Type`
## integer(0)
## 
## $dlptot
## integer(0)

dfx <- na.omit(dfx)

#implement MatchIt package
mod_match <- matchit(tbf ~ ., method="nearest",data=dfx)

#create a df containing only the matched observations

df_m <- match.data(mod_match) #distance is the propensity score
dim(df_m) #dimensions of the data to see how it changed

## [1] 96 12

#96 is the 48 tbf samples matched to 48 similar non tbf samples

#visual inspection
#write a function that will plot each propensity score
                 
fn_bal <- function(df, v) {
  df$v <- df[, v]
  df$tbf <- as.factor(df$tbf)
  support <- c(min(df$v), max(df$v))
  ggplot(df, aes(x = distance, y = unlist(v), color = tbf)) +
    geom_point(alpha = 0.2, size = 1.3) +
    geom_smooth(method = "loess", se = F) +
    xlab("Propensity score") +
    ylab(v) +
    theme_bw() +
    ylim(support)
}

#run the function and plot the scores

grid.arrange(
   fn_bal(df_m, "Age")+ theme(legend.position = "none"),
   fn_bal(df_m, "sex") + theme(legend.position = "none"),
   fn_bal(df_m, "Transfer"),
   fn_bal(df_m, "LOS (d)") + theme(legend.position = "none"),
   fn_bal(df_m, "Mechanism")+ theme(legend.position = "none"),
   fn_bal(df_m, "GCS_bin"),
   fn_bal(df_m, "Hemorrhage Type"),ncol=3
)

## # A tibble: 2 × 9
##   tbf_cat        Age   sex Transfer `LOS (d)` Mechanism GCS_bin `Hemorrhage Typ…
##   <chr>        <dbl> <dbl>    <dbl>     <dbl>     <dbl>   <dbl>            <dbl>
## 1 Other Fract…  1.66     1        1       2           1       1                1
## 2 Temp Bone F… 12.0      1        1       4.5         3       0                1
## # … with 1 more variable: dlptot <dbl>

## # A tibble: 2 × 9
##     tbf   Age   sex Transfer `LOS (d)` Mechanism GCS_bin `Hemorrhage Typ… dlptot
##   <dbl> <dbl> <dbl>    <dbl>     <dbl>     <dbl>   <dbl>            <dbl>  <dbl>
## 1     0  12.9     1        1       3           2       1                0   873.
## 2     1  12.0     1        1       4.5         3       0                1  1016.

## Warning in wilcox.test.default(x = c(1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, :
## cannot compute exact p-value with ties

## Warning in wilcox.test.default(x = c(0, 0, 0, 1, 1, 0, 1, 1, 1, 0, 0, 0, :
## cannot compute exact p-value with ties

## Warning in wilcox.test.default(x = c(0.9, 4, 0.9, 35, 0.9, 0.9, 0.9, 0.9, :
## cannot compute exact p-value with ties

## Warning in wilcox.test.default(x = c(3, 2, 5, 1, 6, 1, 6, 1, 5, 2, 2, 2, :
## cannot compute exact p-value with ties

## Warning in wilcox.test.default(x = c(1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, :
## cannot compute exact p-value with ties

## Warning in wilcox.test.default(x = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, :
## cannot compute exact p-value with ties

## Warning in wilcox.test.default(x = c(0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, :
## cannot compute exact p-value with ties

## Warning in wilcox.test.default(x = c(1354, 383, 732, 790, 1475.12, 576, : cannot
## compute exact p-value with ties

##                        [,1]
## Age.NA             0.806966
## sex.NA             0.517714
## Transfer.NA        0.831632
## LOS (d).NA         0.250655
## Mechanism.NA       0.544776
## GCS_bin.NA         0.687768
## tbf.NA             0.000000
## Hemorrhage Type.NA 0.252425
## dlptot.NA          0.229458

#estimate treatment effects with a ttest
#there are ties in data, so use exact=false
with(df_m, wilcox.test(dlptot ~ tbf,exact=FALSE))

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  dlptot by tbf
## W = 987.5, p-value = 0.2295
## alternative hypothesis: true location shift is not equal to 0

#or use Ordinary Least Squares without covariates
m2 <- lm(dlptot ~ tbf, data = df_m)
summary(m2)

## 
## Call:
## lm(formula = dlptot ~ tbf, data = df_m)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1172.3  -441.3  -154.1   338.7  1828.2 
## 
## Coefficients:
##             Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)   1004.4      100.9   9.950 0.000000000000000234 ***
## tbf            167.9      142.7   1.176                0.242    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 699.3 on 94 degrees of freedom
## Multiple R-squared:  0.01451,    Adjusted R-squared:  0.004024 
## F-statistic: 1.384 on 1 and 94 DF,  p-value: 0.2424

#linear model with covs
m3 <- lm(dlptot ~ .-subclass -weights -distance, data = df_m)
summary(m3)

## 
## Call:
## lm(formula = dlptot ~ . - subclass - weights - distance, data = df_m)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1102.88  -474.56   -79.59   451.38  1600.06 
## 
## Coefficients:
##                   Estimate Std. Error t value   Pr(>|t|)    
## (Intercept)       -147.774    251.332  -0.588     0.5581    
## Age                 57.513     12.116   4.747 0.00000805 ***
## sex                -63.634    132.785  -0.479     0.6330    
## Transfer           314.671    138.411   2.273     0.0255 *  
## `LOS (d)`            6.316      6.099   1.035     0.3033    
## Mechanism           64.806     36.064   1.797     0.0758 .  
## GCS_bin             14.386    149.081   0.096     0.9233    
## tbf                132.225    123.882   1.067     0.2888    
## `Hemorrhage Type`  191.780    106.291   1.804     0.0746 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 602.3 on 87 degrees of freedom
## Multiple R-squared:  0.3235, Adjusted R-squared:  0.2613 
## F-statistic:   5.2 on 8 and 87 DF,  p-value: 0.00002541

tib2 <- df_m %>%
  group_by(tbf) %>%
            summarise(n_patients=n(),mean_DLP = mean(dlptot),
            n_patients=n(),median_DLP = median(dlptot),
            std_error = sd(dlptot) / sqrt(n_patients))

## `summarise()` ungrouping output (override with `.groups` argument)

tib2

## # A tibble: 2 × 5
##     tbf n_patients mean_DLP median_DLP std_error
##   <dbl>      <int>    <dbl>      <dbl>     <dbl>
## 1     0         48    1004.       873.      97.4
## 2     1         48    1172.      1016.     104.

#visual summariy boxplots
ggboxplot(df_m, x = "tbf", y = "dlptot", lwd=.6, show.legend = F,
          color = "tbf", palette = c("pink2","paleturquoise3"),
          ylab = "DLP (mGy*cm)", xlab = "TBF",width=.7,bxp.errorbar.width = 4)+geom_point(size=2,alpha=.5,position = position_jitterdodge(jitter.width = .6), aes(color = factor(tbf)))+ scale_x_discrete(name ="", labels=c("0" = "Other \nfracture\ntype", "1" = "Temporal \nbone\nFracture"))+ggtitle("DLP by Fracture Type \nafter Score Matching") +   scale_y_continuous(breaks = round(seq(min(df$dlptot), max(df$dlptot), by = 400),1))+theme_gray(base_size = 14)+theme(legend.position = "none")+gghisto

ggboxplot(df_m, x = "tbf", y = "Age", lwd=.6, show.legend = F,
          color = "tbf", palette = c("pink2","paleturquoise3"),
          ylab = "Age", xlab = "TBF",width=.7,bxp.errorbar.width = 4)+geom_point(size=2,alpha=.5,position = position_jitterdodge(jitter.width = 0.6),aes(color = factor(tbf)))+ scale_x_discrete(name ="", labels=c("0" = "Other \nfracture\ntype", "1" = "Temporal \nbone\nFracture"))+ggtitle("Age by Fracture Type") +   scale_y_continuous(breaks = round(seq(min(df$Age), max(df$Age), by = 2),1))+theme_gray(base_size = 14)+theme(legend.position = "none")+gghisto

ggboxplot(df_m, x = "tbf", y = "dlptot", lwd=.6, show.legend = F,
          color = "tbf", palette = c("indianred3", "royalblue1"),
          ylab = "DLP", xlab = "TBF",width=.35,bxp.errorbar.width = 4)+geom_point(size=2,alpha=.5,position = position_jitterdodge(jitter.width = 0.6),aes(color = factor(tbf)))+ scale_x_discrete(name ="", labels=c("0" = "Other \nfracture\ntype", "1" = "Temporal \nbone\nFracture"))+ggtitle("DLP by Fracture Type") +theme_gray(base_size = 14)+theme(legend.position = "none")+gghisto

df_m <- match.data(mod_match)

df_m <- add_column(df_m, GCS_cat = ifelse((df_m$GCS_bin==15), "15", "Below 15"), .after = 6)


df_m  <- add_column(df_m, mechanism_cat = df_m$Mechanism, .after = 5)

df_m$mechanism_cat[df_m$Mechanism == 1] <- "Fall"
df_m$mechanism_cat[df_m$Mechanism == 2] <- "Motor Vehicle Collision"
df_m$mechanism_cat[df_m$Mechanism == 3] <- "Car against Pedestrian"
df_m$mechanism_cat[df_m$Mechanism == 4] <- "Bike, Scooter, Skateboard, etc"
df_m$mechanism_cat[df_m$Mechanism == 5] <- "Blunt Injury"
df_m$mechanism_cat[df_m$Mechanism == 6] <- "ATV"
df_m$mechanism_cat[df_m$Mechanism == 7] <- "Other"

df_m$mechanism_cat <- as.factor(df_m$mechanism_cat)

kruskal.test(dlptot ~ mechanism_cat, data = df_m)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  dlptot by mechanism_cat
## Kruskal-Wallis chi-squared = 19.242, df = 6, p-value = 0.003775

pairwise.wilcox.test(df_m$dlptot, df_m$mechanism_cat,
                 p.adjust.method = "BH", exact=FALSE)

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  df_m$dlptot and df_m$mechanism_cat 
## 
##                                ATV   Bike, Scooter, Skateboard, etc
## Bike, Scooter, Skateboard, etc 0.404 -                             
## Blunt Injury                   0.277 0.460                         
## Car against Pedestrian         0.248 0.277                         
## Fall                           0.035 0.052                         
## Motor Vehicle Collision        0.114 0.400                         
## Other                          0.035 0.039                         
##                                Blunt Injury Car against Pedestrian Fall 
## Bike, Scooter, Skateboard, etc -            -                      -    
## Blunt Injury                   -            -                      -    
## Car against Pedestrian         0.967        -                      -    
## Fall                           0.675        0.572                  -    
## Motor Vehicle Collision        0.572        0.646                  0.157
## Other                          0.221        0.039                  0.157
##                                Motor Vehicle Collision
## Bike, Scooter, Skateboard, etc -                      
## Blunt Injury                   -                      
## Car against Pedestrian         -                      
## Fall                           -                      
## Motor Vehicle Collision        -                      
## Other                          0.035                  
## 
## P value adjustment method: BH

cor.test(df_m$dlptot, df_m$Age, 
                    method = "pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  df_m$dlptot and df_m$Age
## t = 4.4678, df = 94, p-value = 0.00002201
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2380047 0.5710877
## sample estimates:
##       cor 
## 0.4185203

ggplot(df_m, aes(x = Age, y = dlptot)) + geom_point(aes(color = tbf))+
  geom_smooth(method = "lm") +   scale_y_continuous(breaks = round(seq(min(df$dlptot), max(df$dlptot), by = 400),1))+   scale_x_continuous(breaks = round(seq(min(df$Age), max(df$Age), by = 1),1))+
  gghisto+ggtitle("Total Radiation by Age \nand Fracture Type")

## `geom_smooth()` using formula 'y ~ x'

marty temporal bone

alice

10/22/2021