Analysis on Radiological Research

Load Excel Data

library(readxl)
library(plyr)

## Warning: package 'plyr' was built under R version 3.4.4

library(pROC)

## Warning: package 'pROC' was built under R version 3.4.4

## Type 'citation("pROC")' for a citation.

## 
## Attaching package: 'pROC'

## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var

Data <- read_excel("C:/Users/Roel Ceballos/Downloads/Data.xlsx")
View(Data)

Subsetting Data (Vertebral Fracture Vs No Vertebral Fracture) Mann-Whitney U test

newdata1=subset(Data, Vertebra=="No",)
newdata2=subset(Data, Vertebra=="Yes",)
G1mean=mean(newdata1$Age)
G2mean=mean(newdata2$Age)

wilcox.test(newdata1$Age,newdata2$Age)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  newdata1$Age and newdata2$Age
## W = 68796, p-value = 2.535e-12
## alternative hypothesis: true location shift is not equal to 0

hipdata1=subset(Data,Hip=="No",)
hipdata2=subset(Data,Hip=="Yes",)
H1mean=mean(hipdata1$Age)
H2mean=mean(hipdata2$Age)

wilcox.test(hipdata1$Age,hipdata2$Age)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  hipdata1$Age and hipdata2$Age
## W = 18758, p-value = 0.603
## alternative hypothesis: true location shift is not equal to 0

hmean=mean(hipdata1$`L5 Trabecular Attenuation (HU)`)
hmean2=mean(hipdata2$`L5 Trabecular Attenuation (HU)`)
print(hmean)

## [1] 123.6264

print(hmean2)

## [1] 52.60976

Chisquare test

chisq.test(Data$Sex,Data$Vertebra)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  Data$Sex and Data$Vertebra
## X-squared = 19.261, df = 1, p-value = 1.14e-05

chisq.test(Data$`CT Study Type`,Data$Vertebra)

## Warning in chisq.test(Data$`CT Study Type`, Data$Vertebra): Chi-squared
## approximation may be incorrect

## 
##  Pearson's Chi-squared test
## 
## data:  Data$`CT Study Type` and Data$Vertebra
## X-squared = 10.308, df = 5, p-value = 0.06697

chisq.test(Data$Vertebra,Data$`Use of Contrast Material`)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  Data$Vertebra and Data$`Use of Contrast Material`
## X-squared = 0.096831, df = 1, p-value = 0.7557

chisq.test(Data$Sex, Data$Hip)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  Data$Sex and Data$Hip
## X-squared = 0.031502, df = 1, p-value = 0.8591

chisq.test(Data$`CT Study Type`,Data$Hip)

## Warning in chisq.test(Data$`CT Study Type`, Data$Hip): Chi-squared
## approximation may be incorrect

## 
##  Pearson's Chi-squared test
## 
## data:  Data$`CT Study Type` and Data$Hip
## X-squared = 19.161, df = 5, p-value = 0.001794

chisq.test(Data$`Use of Contrast Material`,Data$Hip)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  Data$`Use of Contrast Material` and Data$Hip
## X-squared = 3.5403, df = 1, p-value = 0.05989

Table 2 L5 VertFRact NoVertFRact TRUE A=TP B=FP FALSE C=FN D= TN Total n1=262 n2=740 Sensitivity = A/(A+C) Specificity = D/(B+D) Accuracy=(a+d/1002) Prevalence = (A+B)/(A+B+C+D) PPV = TP/(TP+FP) NPV = TN/(TN+FN)

for(i in c(50, 70, 90,96.5,100,110,130,150,170)){
    dat=subset(Data, Data$`L5 Trabecular Attenuation (HU)`<=i)
    b=sum(dat$Vertebra=="No")
    a=sum(dat$Vertebra=="Yes")
    dat2=subset(Data, Data$`L5 Trabecular Attenuation (HU)`>i)
    d=sum(dat2$Vertebra=="No")
    c=sum(dat2$Vertebra=="Yes")
    
    print(paste("Group",i))
    print(c(a,b,c,d))
    Sen=a/(a+c)
    Spec=d/(b+d)
    Prev=(a+b)/(1002)
    acc=(a+d)/1002
    PPV=(Sen*Prev)/((Sen*Prev)+(1-Spec)*(1-Prev))
    NPV=(Spec*(1-Prev))/((1-Sen*Prev)+(Spec)*(1-Prev))
    print(paste("Sensitivity is ", Sen))
    print(paste("Specificity is ", Spec))
    print(paste("Prevalence is", Prev))
    print(paste("accuracy is ",acc))
    print(paste("PPV",a/(a+b)))
    print(paste("NPV",d/(c+d)))

}

## [1] "Group 50"
## [1]  63  12 199 728
## [1] "Sensitivity is  0.240458015267176"
## [1] "Specificity is  0.983783783783784"
## [1] "Prevalence is 0.0748502994011976"
## [1] "accuracy is  0.789421157684631"
## [1] "PPV 0.84"
## [1] "NPV 0.785329018338727"
## [1] "Group 70"
## [1] 140  31 122 709
## [1] "Sensitivity is  0.534351145038168"
## [1] "Specificity is  0.958108108108108"
## [1] "Prevalence is 0.170658682634731"
## [1] "accuracy is  0.847305389221557"
## [1] "PPV 0.818713450292398"
## [1] "NPV 0.853188929001203"
## [1] "Group 90"
## [1] 219  91  43 649
## [1] "Sensitivity is  0.83587786259542"
## [1] "Specificity is  0.877027027027027"
## [1] "Prevalence is 0.30938123752495"
## [1] "accuracy is  0.86626746506986"
## [1] "PPV 0.706451612903226"
## [1] "NPV 0.937861271676301"
## [1] "Group 96.5"
## [1] 230 115  32 625
## [1] "Sensitivity is  0.877862595419847"
## [1] "Specificity is  0.844594594594595"
## [1] "Prevalence is 0.344311377245509"
## [1] "accuracy is  0.853293413173653"
## [1] "PPV 0.666666666666667"
## [1] "NPV 0.951293759512938"
## [1] "Group 100"
## [1] 234 131  28 609
## [1] "Sensitivity is  0.893129770992366"
## [1] "Specificity is  0.822972972972973"
## [1] "Prevalence is 0.364271457085828"
## [1] "accuracy is  0.841317365269461"
## [1] "PPV 0.641095890410959"
## [1] "NPV 0.956043956043956"
## [1] "Group 110"
## [1] 242 198  20 542
## [1] "Sensitivity is  0.923664122137405"
## [1] "Specificity is  0.732432432432432"
## [1] "Prevalence is 0.439121756487026"
## [1] "accuracy is  0.782435129740519"
## [1] "PPV 0.55"
## [1] "NPV 0.9644128113879"
## [1] "Group 130"
## [1] 254 335   8 405
## [1] "Sensitivity is  0.969465648854962"
## [1] "Specificity is  0.547297297297297"
## [1] "Prevalence is 0.587824351297405"
## [1] "accuracy is  0.657684630738523"
## [1] "PPV 0.431239388794567"
## [1] "NPV 0.980629539951574"
## [1] "Group 150"
## [1] 258 462   4 278
## [1] "Sensitivity is  0.984732824427481"
## [1] "Specificity is  0.375675675675676"
## [1] "Prevalence is 0.718562874251497"
## [1] "accuracy is  0.534930139720559"
## [1] "PPV 0.358333333333333"
## [1] "NPV 0.985815602836879"
## [1] "Group 170"
## [1] 261 571   1 169
## [1] "Sensitivity is  0.99618320610687"
## [1] "Specificity is  0.228378378378378"
## [1] "Prevalence is 0.830339321357285"
## [1] "accuracy is  0.429141716566866"
## [1] "PPV 0.313701923076923"
## [1] "NPV 0.994117647058824"

print("Hip Fracture -----------------------------------------------------------------------")

## [1] "Hip Fracture -----------------------------------------------------------------------"

for(i in c(50, 70, 90,67.5,100,110,130,150,170)){
    dat=subset(Data, Data$`L5 Trabecular Attenuation (HU)`<=i)
    b=sum(dat$Hip=="No")
    a=sum(dat$Hip=="Yes")
    dat2=subset(Data, Data$`L5 Trabecular Attenuation (HU)`>i)
    d=sum(dat2$Hip=="No")
    c=sum(dat2$Hip=="Yes")
    
    print(paste("Group",i))
    print(c(a,b,c,d))
    Sen=a/(a+c)
    Spec=d/(b+d)
    Prev=(a+b)/(1002)
    acc=(a+d)/1002
    PPV=(Sen*Prev)/((Sen*Prev)+(1-Spec)*(1-Prev))
    NPV=(Spec*(1-Prev))/((1-Sen*Prev)+(Spec)*(1-Prev))
    print(paste("Sensitivity is ", Sen))
    print(paste("Specificity is ", Spec))
    print(paste("Prevalence is", Prev))
    print(paste("accuracy is ",acc))
    print(paste("PPV",a/(a+b)))
    print(paste("NPV",d/(c+d)))

}

## [1] "Group 50"
## [1]  22  53  19 908
## [1] "Sensitivity is  0.536585365853659"
## [1] "Specificity is  0.944849115504683"
## [1] "Prevalence is 0.0748502994011976"
## [1] "accuracy is  0.92814371257485"
## [1] "PPV 0.293333333333333"
## [1] "NPV 0.97950377562028"
## [1] "Group 70"
## [1]  38 133   3 828
## [1] "Sensitivity is  0.926829268292683"
## [1] "Specificity is  0.861602497398543"
## [1] "Prevalence is 0.170658682634731"
## [1] "accuracy is  0.864271457085828"
## [1] "PPV 0.222222222222222"
## [1] "NPV 0.996389891696751"
## [1] "Group 90"
## [1]  38 272   3 689
## [1] "Sensitivity is  0.926829268292683"
## [1] "Specificity is  0.716961498439126"
## [1] "Prevalence is 0.30938123752495"
## [1] "accuracy is  0.725548902195609"
## [1] "PPV 0.12258064516129"
## [1] "NPV 0.995664739884393"
## [1] "Group 67.5"
## [1]  38 124   3 837
## [1] "Sensitivity is  0.926829268292683"
## [1] "Specificity is  0.870967741935484"
## [1] "Prevalence is 0.161676646706587"
## [1] "accuracy is  0.873253493013972"
## [1] "PPV 0.234567901234568"
## [1] "NPV 0.996428571428571"
## [1] "Group 100"
## [1]  38 327   3 634
## [1] "Sensitivity is  0.926829268292683"
## [1] "Specificity is  0.659729448491155"
## [1] "Prevalence is 0.364271457085828"
## [1] "accuracy is  0.67065868263473"
## [1] "PPV 0.104109589041096"
## [1] "NPV 0.995290423861852"
## [1] "Group 110"
## [1]  40 400   1 561
## [1] "Sensitivity is  0.975609756097561"
## [1] "Specificity is  0.583766909469303"
## [1] "Prevalence is 0.439121756487026"
## [1] "accuracy is  0.599800399201597"
## [1] "PPV 0.0909090909090909"
## [1] "NPV 0.998220640569395"
## [1] "Group 130"
## [1]  40 549   1 412
## [1] "Sensitivity is  0.975609756097561"
## [1] "Specificity is  0.428720083246618"
## [1] "Prevalence is 0.587824351297405"
## [1] "accuracy is  0.451097804391218"
## [1] "PPV 0.067911714770798"
## [1] "NPV 0.997578692493947"
## [1] "Group 150"
## [1]  41 679   0 282
## [1] "Sensitivity is  1"
## [1] "Specificity is  0.293444328824142"
## [1] "Prevalence is 0.718562874251497"
## [1] "accuracy is  0.322355289421158"
## [1] "PPV 0.0569444444444444"
## [1] "NPV 1"
## [1] "Group 170"
## [1]  41 791   0 170
## [1] "Sensitivity is  1"
## [1] "Specificity is  0.176899063475546"
## [1] "Prevalence is 0.830339321357285"
## [1] "accuracy is  0.210578842315369"
## [1] "PPV 0.0492788461538462"
## [1] "NPV 1"

Boxplot

boxplot(Data$`L5 Trabecular Attenuation (HU)`~Data$Vertebra,xlab="Vertebral Fracture",ylab="L5 Attenuation (HU)")

boxplot(Data$`L5 Trabecular Attenuation (HU)`~Data$Hip,xlab="Hip Fracture", ylab="L5 Attenuation (HU)")

ROC Curve

plot.roc(Data$Vertebra,Data$`L5 Trabecular Attenuation (HU)`, main="Vertebral Fracture Attenuation (HU) threshold", percent=TRUE, ci=TRUE, of="thresholds",thresholds="best", print.thres="best")

plot.roc(Data$Hip,Data$`L5 Trabecular Attenuation (HU)`, main="Hip Fracture Attenuation (HU) threshold", percent=TRUE, ci=TRUE, of="thresholds", thresholds="best", print.thres="best")

Relative Frequency

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 3.4.4

library(reshape2)

## Warning: package 'reshape2' was built under R version 3.4.4

dat=list(Vertebral_Fracture=newdata2$`L5 Trabecular Attenuation (HU)`,No_Fracture=newdata1$`L5 Trabecular Attenuation (HU)`)
data<-melt(dat)
ggplot(data,aes(x=value, fill=L1)) + geom_density(alpha=0.25)

dat2=list(Hip_Fracture=hipdata2$`L5 Trabecular Attenuation (HU)`,No_Fracture=hipdata1$`L5 Trabecular Attenuation (HU)`)
data1<-melt(dat2)
ggplot(data1,aes(x=value,fill=L1))+geom_density(alpha=0.25)

Dot Plot

library(ggplot2)
kl=data.frame(Age=Data$Age, L5=Data$`L5 Trabecular Attenuation (HU)`, Vertebral_Fracture=Data$Vertebra, ncol=3)
ggplot(kl,aes(x=Age,y=L5, color=Vertebral_Fracture))+geom_point()+labs(x="Age",y="L5 Attenuation")+geom_smooth(method='lm',se=FALSE, fullrange=TRUE)

lp=data.frame(Age=Data$Age, L5=Data$`L5 Trabecular Attenuation (HU)`,Hip_Fracture=Data$Hip,ncol=3)
ggplot(lp,aes(x=Age,y=L5,color=Hip_Fracture))+geom_point()+labs(x="Age",y="L5 Attenuation (HU)")+geom_smooth(method='lm',se=FALSE, fullrange=TRUE)

Hip Fracture

Analysis on Radiological Research

Roel Ceballos

5/30/2019