Introduction

This vignette introduces the use of FRESA.CAD to create Models of Kellgren and Lawrence Scores (KL) and predictors of total Knee replacements (TKR)

Preparation

Loading the libraries



library("epiR")
library("FRESA.CAD")
library(network)
library(GGally)
library("R.matlab")
library("gplots")
library(rpart)
library("randomForest")
library("e1071")

a=as.numeric(Sys.time());
set.seed(a);

error.bar <- function(x, y, upper, lower=upper, length=0.05,...){
if(length(x) != length(y) | length(y) !=length(lower) | length(lower) != length(upper))
stop("vectors must be same length")
arrows(x,y+upper, x, y-lower, angle=90, code=3, length=length, ...)
}

barPlotCiError<-  function(citable,metricname,thesets,themethod,main,...)
{
colnames(citable) <- c(metricname,"lower","upper")
rownames(citable) <- rep(thesets,length(themethod))
pander::pander(citable,caption=main,round = 3)
citable <- citable[order(rep(1:length(thesets),length(themethod))),]
barmatrix <- matrix(citable[,1],length(themethod),length(thesets))
colnames(barmatrix) <- thesets
rownames(barmatrix) <- themethod
pander::pander(barmatrix,caption=main,round = 3)
barp <- barplot(barmatrix,cex.names=0.7,las=2,ylim=c(0.0,1.0),main=main,ylab=metricname,beside=TRUE,legend = themethod,...)
error.bar(barp,citable[,1],citable[,3]-citable[,1],citable[,1]-citable[,2])
return(barp)
}

Heat Maps Auxilary Plots

# the colors of the type of JSW measurement. These are categorical description of the bars 
lType <- c("Original","Derived")
lTypecolors <- c("green","Orange")

feattype <- c("Medial","Lateral","Distance","Location","Joint")
featcolor <- c("green","blue","orange","yellow","red")



# the rowcolors is a global variable that will have the ROI colors of the features
# lJSNBar is a global variable that will have the lateral JSN 
# mJSNBar is a global variable that will have the medial JSN 
# ConData is a global variable that will have the total ConData scores
# PainMedication is a global variable that will have the pain Medication

# the colors for the JSN scores and KL scores. These are ordinal
colorJSNpalette <- c("blue","cyan","red","yellow")

# the colors for the heatmap These are continous
colorpalette <- c("cyan","blue","black","red","yellow")

# the Layoout of the plot will be on 7 rows and four columns
# Color Key     |             |            |Col Dendogram
#               |             |            |ColSideColors
#               |             |            |Lateral JSN
#               |             |            |Medial  JSN
#               |             |            |Total ConData
#               |             |            |Medication
# Row Dendogram |RowSideCOlors|Feature Type|HeatMap
#               |             |            | 

lmat <- rbind( c(6,0, 0,5), c(0,0,0,2),c(0,0,0,8),c(0,0,0,9),c(0,0,0,10),c(0,0,0,11), c(4,1,7,3), c(0,0,0,0) )
# the height is 5.0in
lhei <- c(0.9, 0.15,0.15,0.15,0.20,0.10,3.0,0.20) 
# the width is 6.0in
lwid <- c(1.0, 0.15,0.15,4.2 )

myplot <- function() {
            op <- par(no.readonly = TRUE)
              rindx <- get('rowInd',parent.frame())
              cindx <- get('colInd',parent.frame())
              marg <- get('margins',parent.frame())
              cnames <- rownames(get('orderData',parent.frame(2)))
              nr <-  length(featureColor)
              par(mar = c(marg[1], 0, 0, 0.5),cex = 0.66)
              image(rbind(1:nr), col = featureColor[rindx], axes = FALSE)

             nr <-  length(cnames)
            index <- pmin(floor(nr*(0.9999999*mJSNBar[cnames])/4)+1,nr);
             colcolors <- colorRampPalette(colorJSNpalette)(n = nr)[index];
             par(mar = c(0.5, 0, 0, marg[2]),cex = 0.66)
             image(cbind(1:nr), col = colcolors[cindx], axes = FALSE)
             mtext("mJSN", side=4, line=1, cex=0.4,las=1, col="black")

             index <- pmin(floor(nr*(0.9999999*lJSNBar[cnames])/4)+1,nr);
             colcolors <- colorRampPalette(colorJSNpalette)(n = nr)[index];
             par(mar = c(0.5, 0, 0, marg[2]),cex = 0.66)
             image(cbind(1:nr), col = colcolors[cindx], axes = FALSE)
             mtext("lJSN", side=4, line=1, cex=0.4,las=1, col="black")

             
             par(mar = c(0.5, 0, 0, marg[2]),mgp=c(0,0,0),pty="m",xaxs="i")
             bardata <-  ConData[cnames];
#             barplot(bardata,axes=FALSE,axisnames=FALSE,space=NULL)
             plot(c(1:nr),bardata,type = "s",axes=FALSE,ylab="",xlab="")
             mtext(ConDataName, side=4, line=1, cex=0.4,las=1, col="black")

             index <- pmin(floor(nr*(0.9999999*PainMedication[cnames]))+1,nr);
             colcolors <- colorRampPalette(colorJSNpalette)(n = nr)[index];
             par(mar = c(0.5, 0, 0, marg[2]),cex = 0.66)
             image(cbind(1:nr), col = colcolors[cindx], axes = FALSE)
             mtext("Medication", side=4, line=1, cex=0.4,las=1, col="black")
             
             
              par(op)

}

Loading the data sets.

#Load the JSW Data
OAI_JSW <- read.delim("OAI_JSW_Analysis.txt",na.strings=c("NA","","#N/A","#DIV/0!"))
# just substract 1 from the knee side variable
OAI_JSW$SIDEC = OAI_JSW$SIDE - 1
#Load the Variables to Be processed
OAI_JSWVarDescription <- read.delim("OAI_JSWDescription.txt")
rownames(OAI_JSWVarDescription) <- OAI_JSWVarDescription$Var

Conditioning the Data

JSWtimecontrol <- subset(OAI_JSW, BLKL == 0 & Y2DeltaJSN==0 & BLTOTALWOMAC < 10)

Adjust each site drift


OAI_JSW.adj <- featureAdjustment(OAI_JSWVarDescription, baseModel="1+DaysFomStart", strata = "SITE",data=OAI_JSW,JSWtimecontrol, type = "GLS", pvalue = 0.05,correlationGroup = "IDSIDE")


#lets use the subjects that did not progress to get anthrometric (height) and gender associations
JSWcontrol.BL <- subset(OAI_JSW.adj, VISITID==0 & BLKL == 0 &  Y2DeltaJSN==0 & BLTOTALWOMAC < 10)
#JSWcontrol.BL <- JSWcontrol.BL[complete.cases(JSWcontrol.BL),]



#We adjust for those differences
OAI_JSW.full.adj <- featureAdjustment(OAI_JSWVarDescription, baseModel="1+AGE+HEIGHT+CFWDTH
", strata = "Gender", data=OAI_JSW.adj,JSWcontrol.BL, type = "LM", pvalue = 0.05)

Quantile Normalization


JSWBLadjc <- subset(OAI_JSW.full.adj,SIDE==2 & VISITID==0 & BLKL == 0 & Y2DeltaJSN==0 & BLTOTALWOMAC < 10)
JSWBLadjc <- JSWBLadjc[complete.cases(JSWBLadjc),]
OAI_JSW.full.adj.norm <- rankInverseNormalDataFrame(OAI_JSWVarDescription, OAI_JSW.full.adj, JSWBLadjc,strata="Gender")

I’ll check the longitudinal behavior. y-axis units are standard deviations.

par(mfrow=c(4,4),cex=0.25)
cs1CorCAR1 <- nlme::corCAR1(0.99,form = ~ YearsFromBL | IDSIDE)

OAI_JSW.full.adj.norm$RADOA <- 1*(OAI_JSW.full.adj.norm$BLKL >= 2)
OAI_JSW.full.adj.norm$YearsFromBL <- OAI_JSW.full.adj.norm$DAYSFROMBL/356.25
OAI_JSWVarDescription$FullDescription <- paste("Normalized: ",OAI_JSWVarDescription$Description)
tsa <- timeSerieAnalysis(OAI_JSWVarDescription,"I(YearsFromBL-2)*RADOA",subset(OAI_JSW.full.adj.norm,BLKL<4),timevar="VISITID",contime="YearsFromBL",Outcome="RADOA", cs1CorCAR1,timesign= "+",Ptoshow = c(2,3,4),plegend= c("Visit","ROA","Visit:ROA"),catgo.names=c("Normal", "ROA"),description="FullDescription") 

par(mfrow=c(1,1),cex=1)

tvaluesSiteAgeAdjustedNormalized <- tsa$t.values

Modeling the KL score with Linear Models

Creating a train and a test set


baseline <- subset(OAI_JSW.full.adj.norm,VISITID == 0)

#load(file = "OAIJSW_BASELINE.RDATA")

rownames(baseline) <- baseline$IDSIDE
baseline <- baseline[complete.cases(baseline),c("BLKL","AGE","BMI","WEIGHT",as.character(OAI_JSWVarDescription$Var))]

JSWKL0 <- subset(baseline,BLKL == 0)
JSWKL1 <- subset(baseline,BLKL == 1)
JSWKL2 <- subset(baseline,BLKL == 2)
JSWKL3 <- subset(baseline,BLKL == 3)
JSWKL4 <- subset(baseline,BLKL == 4)
sizePerKL <- nrow(JSWKL4)/2
trainIDs_0 <- sample(nrow(JSWKL0),sizePerKL)
trainIDs_1 <- sample(nrow(JSWKL1),sizePerKL)
trainIDs_2 <- sample(nrow(JSWKL2),sizePerKL)
trainIDs_3 <- sample(nrow(JSWKL3),sizePerKL)
trainIDs_4 <- sample(nrow(JSWKL4),sizePerKL)
trainSubset <- rbind(JSWKL0[trainIDs_0,],JSWKL1[trainIDs_1,],JSWKL2[trainIDs_2,],JSWKL3[trainIDs_3,],JSWKL4[trainIDs_4,])
test2Subset <- rbind(JSWKL0[-trainIDs_0,],JSWKL1[-trainIDs_1,],JSWKL2[-trainIDs_2,],JSWKL3[-trainIDs_3,],JSWKL4[-trainIDs_4,])

testIDs_0 <- sample(c(1:nrow(JSWKL0))[-trainIDs_0],sizePerKL)
testIDs_1 <- sample(c(1:nrow(JSWKL1))[-trainIDs_1],sizePerKL)
testIDs_2 <- sample(c(1:nrow(JSWKL2))[-trainIDs_2],sizePerKL)
testIDs_3 <- sample(c(1:nrow(JSWKL3))[-trainIDs_3],sizePerKL)
testIDs_4 <- c(1:nrow(JSWKL4))[-trainIDs_4]

testSubset <- rbind(JSWKL0[testIDs_0,],JSWKL1[testIDs_1,],JSWKL2[testIDs_2,],JSWKL3[testIDs_3,],JSWKL4[testIDs_4,])

Modeling KL Scores

library("irr")

system.time(OAModelKLOM <- BSWiMS.model(BLKL ~ 1,trainSubset,NumberofRepeats = 5))

bswimpredict <- predict(OAModelKLOM$oridinalModels,test2Subset)[,1]
boxplot(bswimpredict~test2Subset$BLKL,main="B:SWiMS Test Result",ylab="Test Prediction",xlab="KL Score")

cor.test(bswimpredict,test2Subset$BLKL,method="kendall")
kappa2(cbind(as.integer(bswimpredict+0.5),test2Subset$BLKL),"squared")
kappa2(cbind(as.integer(bswimpredict+0.5),test2Subset$BLKL))
table(as.integer(bswimpredict+0.5),test2Subset$BLKL)


bswimpredict <- predict(OAModelKLOM$bagging$bagged.model,test2Subset)
boxplot(bswimpredict~test2Subset$BLKL,main="B:SWiMS Test Result",ylab="Test Prediction",xlab="KL Score")

cor.test(bswimpredict,test2Subset$BLKL,method="kendall")
kappa2(cbind(as.integer(bswimpredict+0.5),test2Subset$BLKL),"squared")
kappa2(cbind(as.integer(bswimpredict+0.5),test2Subset$BLKL))
table(as.integer(bswimpredict+0.5),test2Subset$BLKL)

system.time(OAModelKL <- FRESA.Model(BLKL ~ 1,trainSubset,repeats=5,CVfolds=5,testType="Ftest",equivalent = TRUE))

bswimpredict <- predict(OAModelKL$BSWiMS.model,test2Subset)
bswimpredict <- predict(OAModelKL$bagging$bagged.model,test2Subset)
boxplot(bswimpredict~test2Subset$BLKL,main="B:SWiMS Test Result",ylab="Test Prediction",xlab="KL Score")

cor.test(bswimpredict,test2Subset$BLKL,method="kendall")
kappa2(cbind(as.integer(bswimpredict+0.5),test2Subset$BLKL),"squared")
kappa2(cbind(as.integer(bswimpredict+0.5),test2Subset$BLKL))
table(as.integer(bswimpredict+0.5),test2Subset$BLKL)

KendallTable <- NULL;
KappaTable <- NULL;
RMSETable <- NULL;

plot(OAModelKL$bootstrappedModel)


sumBSWIMS <- summary(OAModelKL$BSWiMS.model)
pander::pander(sumBSWIMS)

sumForward <- summary(OAModelKL$bagging$bagged.model)
pander::pander(sumForward)

Evaluating the Models using the normalized data




boxplot(OAModelKL$cvObject$Models.testPrediction[,"Prediction"]~OAModelKL$cvObject$Models.testPrediction[,"Outcome"],main="CV Test Result",ylab="CV Test Prediction",xlab="KL Score")

KendallTable$CVRawBSWIMS <- cor.test(OAModelKL$cvObject$Models.testPrediction[,"Prediction"],OAModelKL$cvObject$Models.testPrediction[,"Outcome"],method="kendall")
KappaTable$CVRawBSWIMS <- epi.kappa(table(OAModelKL$cvObject$Models.testPrediction[,"Prediction"]>2,OAModelKL$cvObject$Models.testPrediction[,"Outcome"]>2.0))
RMSETable$CVRawBSWIMS <- sqrt(mean((OAModelKL$cvObject$Models.testPrediction[,"Prediction"]-OAModelKL$cvObject$Models.testPrediction[,"Outcome"])^2))

tkrKL.rf <- randomForest(BLKL ~ .,trainSubset)
rf.pred <- predict(tkrKL.rf,testSubset)
boxplot(rf.pred~testSubset$BLKL)

KendallTable$RawRanFor <- cor.test(rf.pred,testSubset$BLKL,method="kendall")
KappaTable$RawRanFor <- epi.kappa(table(rf.pred>2,testSubset$BLKL>2.0))
RMSETable$RawRanFor <- sqrt(mean((rf.pred-testSubset$BLKL)^2))

svm.fit <- svm(BLKL ~ .,trainSubset)
svm.pred <- predict(svm.fit,testSubset)
boxplot(svm.pred~testSubset$BLKL)

KendallTable$RawSVM <- cor.test(svm.pred,testSubset$BLKL,method="kendall")
KappaTable$RawSVM <- epi.kappa(table(svm.pred>2,testSubset$BLKL>2.0))
RMSETable$RawSVM <- sqrt(mean((svm.pred-testSubset$BLKL)^2))

OAModelKLS <- FRESA.Model(BLKL ~ 1,trainSubset)

Doing a Ordinal Fit with zIDI Selection Ordinal Fit will be stored in BSWiMS.models\(oridinalModels Use predict(BSWiMS.models\)oridinalModels,testSet) to get the ordinal prediction on a new dataset Features to test: 38 Adjusted Size: 40.89474

Z: 1.281552 Var Max: 40.89474 FitType: LOGIT Test Type: zIDI >>><……

bswims.pred <- predict(OAModelKLS$BSWiMS.model,testSubset)
boxplot(bswims.pred~testSubset$BLKL)

KendallTable$Rawbeq <- cor.test(bswims.pred,testSubset$BLKL,method="kendall")
KappaTable$Rawbeq <- epi.kappa(table(bswims.pred>2,testSubset$BLKL>2.0))
RMSETable$Rawbeq <- sqrt(mean((bswims.pred-testSubset$BLKL)^2))


bswimpredict <- predict(OAModelKL$BSWiMS.model,testSubset)
boxplot(bswimpredict~testSubset$BLKL,main="B:SWiMS Test Result",ylab="Test Prediction",xlab="KL Score")

KendallTable$RawBSWIMS <- cor.test(bswimpredict,testSubset$BLKL,method="kendall")
KappaTable$RawBSWIMS <- epi.kappa(table(bswimpredict>2,testSubset$BLKL>2))
RMSETable$RawBSWIMS <- sqrt(mean((bswimpredict-testSubset$BLKL)^2))
  
bswimpredict <- predict(OAModelKL$cvObject$Fullenet,as.matrix(testSubset[,OAModelKL$cvObject$LassoFilterVarList]))
boxplot(bswimpredict~testSubset$BLKL,main="LASSO Test Result",ylab="Test Prediction",xlab="KL Score")

KendallTable$LASSO <- cor.test(bswimpredict,testSubset$BLKL,method="kendall")
KappaTable$LASSO <- epi.kappa(table(bswimpredict>2,testSubset$BLKL>2))
RMSETable$LASSO <- sqrt(mean((bswimpredict-testSubset$BLKL)^2))


bswimpredict <- ensemblePredict(OAModelKL$eBSWiMS.model$formula.list,trainSubset,testSubset,Outcome ="BLKL",predictType ="linear",type = "LM",nk=-1 )
boxplot(bswimpredict$ensemblePredict~testSubset$BLKL,main="eB:SWiMS Test Result",ylab="Test Prediction",xlab="KL Score")

KendallTable$raweBSWIMS <- cor.test(bswimpredict$ensemblePredict,testSubset$BLKL,method="kendall")
KappaTable$raweBSWIMS <- epi.kappa(table(bswimpredict$ensemblePredict>2,testSubset$BLKL>2))
RMSETable$raweBSWIMS <- sqrt(mean((bswimpredict$ensemblePredict-testSubset$BLKL)^2))


bswimpredict <- ensemblePredict(OAModelKL$cvObject$formula.list,trainSubset,testSubset,Outcome ="BLKL",predictType ="linear",type = "LM",nk=-1 )
boxplot(bswimpredict$ensemblePredict~testSubset$BLKL,main="Ensemble B:SWiMS Test Result",ylab="Test Prediction",xlab="KL Score")

KendallTable$rawEnBSWIMS <- cor.test(bswimpredict$ensemblePredict,testSubset$BLKL,method="kendall")
KappaTable$rawEnBSWIMS <- epi.kappa(table(bswimpredict$ensemblePredict>2,testSubset$BLKL>2))
RMSETable$rawEnBSWIMS <- sqrt(mean((bswimpredict$ensemblePredict-testSubset$BLKL)^2))

The equivalent matrix

pander::pander(OAModelKL$eBSWiMS.model$equivalentMatrix)
Table continues below
Name Locus Extendend_Name UniPerformance
CVJSW CVJSW CVJSW:CVJSW 0.5559
meanLJSW CVJSW meanLJSW:CVJSW 0.09999
LatIntersept CVJSW LatIntersept:CVJSW 0.0694
LJSW725 CVJSW LJSW725:CVJSW 0.1684
LJSW850 CVJSW LJSW850:CVJSW 0.05151
JointSlope CVJSW JointSlope:CVJSW 0.02647
LJSW800 CVJSW LJSW800:CVJSW 0.07591
LJSW750 CVJSW LJSW750:CVJSW 0.134
MeanJSW CVJSW MeanJSW:CVJSW 0.3426
LateralSlope CVJSW LateralSlope:CVJSW 0.1114
meanLJSW AGE meanLJSW:AGE 0.09999
AGE AGE AGE:AGE 0.02987
LJSW850 AGE LJSW850:AGE 0.05151
JSW150 MCMJSW JSW150:MCMJSW 0.2428
MCMJSW MCMJSW MCMJSW:MCMJSW 0.2472
JointSlope MCMJSW JointSlope:MCMJSW 0.02647
JSW300 MCMJSW JSW300:MCMJSW 0.2587
lCV MCMJSW lCV:MCMJSW 0.0539
meanMJSW MCMJSW meanMJSW:MCMJSW 0.2999
MeanJSW MCMJSW MeanJSW:MCMJSW 0.3426
JSW275 MCMJSW JSW275:MCMJSW 0.2868
mCV MCMJSW mCV:MCMJSW 0.1807
meanLJSW BMI meanLJSW:BMI 0.09999
BMI BMI BMI:BMI 0.01963
WEIGHT BMI WEIGHT:BMI 0.0262
LJSW850 BMI LJSW850:BMI 0.05151
JointSlope BMI JointSlope:BMI 0.02647
MtoLDiff LJSW700 MtoLDiff:LJSW700 0.4776
meanLJSW LJSW700 meanLJSW:LJSW700 0.09999
LatIntersept LJSW700 LatIntersept:LJSW700 0.0694
LJSW725 LJSW700 LJSW725:LJSW700 0.1684
LJSW700 LJSW700 LJSW700:LJSW700 0.201
LJSW850 LJSW700 LJSW850:LJSW700 0.05151
JointSlope LJSW700 JointSlope:LJSW700 0.02647
LJSW800 LJSW700 LJSW800:LJSW700 0.07591
LJSW750 LJSW700 LJSW750:LJSW700 0.134
MeanJSW LJSW700 MeanJSW:LJSW700 0.3426
MaxCV MaxCV MaxCV:MaxCV 0.2537
MtoLDiff STDJSW MtoLDiff:STDJSW 0.4776
STDJSW STDJSW STDJSW:STDJSW 0.4665
FullPerformance DeltaPerformance ImprovementFraction p.value
0.6444 0.005497 -0.01132 0.002186
0.6469 0.008011 0.1057 0.0002842
0.6436 0.004717 0.1283 0.004163
0.6424 0.003548 0.06792 0.01112
0.6481 0.009213 0.1358 0.0001082
0.6469 0.008 0.1321 0.0002865
0.6452 0.006304 0.117 0.001131
0.6431 0.004294 0.1358 0.005924
0.6459 0.007056 0.1057 0.000614
0.6454 0.006529 0.09811 0.0009413
0.639 0.006444 0.1132 0.001104
0.6444 0.01183 0.07547 1.571e-05
0.6392 0.006642 0.1094 0.0009425
0.6396 0.02948 0.1434 5.523e-11
0.6444 0.03428 0.1434 1.473e-12
0.6273 0.01719 0.1094 5.132e-07
0.6258 0.01573 0.07925 1.508e-06
0.6253 0.01521 0.06415 2.21e-06
0.6381 0.02801 0.09811 1.669e-10
0.6315 0.02145 0.1396 2.211e-08
0.6311 0.02102 0.08679 3.025e-08
0.6258 0.01571 0.02642 1.535e-06
0.6418 0.004521 0.1283 0.00498
0.6444 0.007025 -0.003774 0.0006447
0.6418 0.004491 0.04151 0.005106
0.6423 0.004966 0.1358 0.003447
0.6422 0.004876 0.1208 0.003713
0.6274 0.005911 0.09434 0.001946
0.6444 0.02286 0.1623 4.644e-09
0.6348 0.0133 0.1811 6.74e-06
0.6457 0.02423 0.1358 1.629e-09
0.6444 0.02286 0.117 4.648e-09
0.6365 0.01504 0.1283 1.799e-06
0.6408 0.01929 0.1321 7.116e-08
0.6372 0.0157 0.1509 1.09e-06
0.643 0.02152 0.1623 1.291e-08
0.6415 0.02002 0.09434 4.066e-08
0.6444 0.00374 0.04151 0.009297
0.6465 0.0131 0.09434 5.698e-06
0.6444 0.01092 0.117 3.141e-05

Heat map of the features


hmdata <- subset(OAI_JSW.full.adj.norm,VISITID==0)
rownames(hmdata) <- hmdata$IDSIDE
hmdata <- hmdata[complete.cases(hmdata),c("BLKL","mJSN","lJSN","BMI","PainMedication","AGE","WEIGHT",as.character(OAI_JSWVarDescription$Var))]
hmdata <- hmdata[sample(nrow(hmdata),250),]
fnames <- unique(as.character(OAModelKL$eBSWiMS.model$equivalentMatrix$Name))
topFeatures <- OAI_JSWVarDescription[fnames,]
topFeatures <- topFeatures[!is.na(topFeatures$Var),]


featureColor <- c("Black",as.character(topFeatures[,"Type.Color"]))

lJSNBar <- hmdata$lJSN
names(lJSNBar) <- rownames(hmdata)
mJSNBar <- hmdata$mJSN
names(mJSNBar) <- rownames(hmdata)
ConData <- hmdata$BMI
ConDataName <- "BMI"
names(ConData) <- rownames(hmdata)
PainMedication <- hmdata$PainMedication
names(PainMedication) <- rownames(hmdata)
hm <- NULL;
KLprediction <- predict(OAModelKL$BSWiMS.model,hmdata,"linear");
#KLprediction <- NULL;

op <- par(no.readonly = TRUE)
hm <- heatMaps(variableList = topFeatures,Outcome = "BLKL",data = hmdata,hCluster= FALSE,prediction=KLprediction,theFiveColors=colorpalette,outcomeColors=colorJSNpalette,transpose=TRUE,xlab="Baseline",ylab="Features",cexRow=0.50,cexCol=0.10,srtCol=35,reorderfun=function(d, w) reorder(d, w, agglo.FUN = mean),distfun=function(c) dist(c, method = "manhattan"),keysize=1.0,RowSideColors=c("Black",as.character(topFeatures[,"Region.Color"])),lmat=lmat, lhei=lhei, lwid=lwid,extrafun=myplot,key.par=list(cex=0.5))
par(new=TRUE,mar=c(1,2,0,0),omd=c(0.02,0.98,0.02,0.98))
plot(c(0,5), c(0,1), type='n', bty='n', xaxt='n', xlab='', yaxt='n', ylab='',main='')
lg <- legend("bottomleft",legend=feattype,fill=featcolor,cex=0.5)
legend(lg$rect$left+lg$rect$w,lg$rect$top,legend=lType,fill=lTypecolors,cex=0.5)

par(op)

Modeling KL using the Out of Normal Range categorization


system.time(OAModelKL3 <- FRESA.Model(BLKL ~ 1,trainSubset,CVfolds=5,repeats=5,equivalent=TRUE,testType="Ftest"))


#OAModelKL3$eBSWiMS.model$equivalentMatrix
#summary(OAModelKL3$BSWiMS.model)
#am <- modelFitting(OAModelKL3$BSWiMS.model$formula,trainSubset,type="LM",fast=TRUE)
#summary.fitFRESA(am)
pp <- plot(OAModelKL3$bootstrappedModel)

pander::pander(summary(OAModelKL3$BSWiMS.model)$coefficients)
Table continues below
  Estimate lower mean upper u.MSE
MtoLDiff 0.1328 0.1198 0.1328 0.1458 1.047
JSW150 -0.0912 -0.1001 -0.0912 -0.08226 1.517
MCMJSW -0.06576 -0.08422 -0.06576 -0.04729 1.508
LJSW700 -0.1178 -0.158 -0.1178 -0.0776 1.601
meanLJSW -0.2481 -0.3348 -0.2481 -0.1614 1.803
AGE 0.00859 0.004546 0.00859 0.01263 1.944
STDJSW 0.08455 0.04314 0.08455 0.126 1.069
BMI 0.0127 0.004965 0.0127 0.02043 1.964
LatIntersept 0.1385 0.04884 0.1385 0.2282 1.865
CVJSW 0.04705 0.01469 0.04705 0.07941 0.8898
MaxCV -0.01167 -0.02138 -0.01167 -0.001951 1.495
Table continues below
  r.MSE model.MSE NeRI F.pvalue t.pvalue
MtoLDiff 1.165 0.7365 0.3208 0 3.331e-30
JSW150 0.8609 0.7365 0.2 0 1.103e-14
MCMJSW 0.7813 0.7126 0.1434 1.473e-12 4.112e-09
LJSW700 0.7584 0.7126 0.117 4.648e-09 0.0001927
meanLJSW 0.7816 0.7365 0.04906 1.025e-08 0.01188
AGE 0.7364 0.7126 0.07547 1.571e-05 4.543e-05
STDJSW 0.7345 0.7126 0.117 3.141e-05 0.001313
BMI 0.7267 0.7126 -0.003774 0.0006447 0.05474
LatIntersept 0.7494 0.7365 0.003774 0.001233 0.209
CVJSW 0.7237 0.7126 -0.01132 0.002186 0.03716
MaxCV 0.7201 0.7126 0.04151 0.009297 0.005618
  Sign.pvalue Wilcox.pvalue Frequency
MtoLDiff 6.592e-14 4.907e-17 0.5
JSW150 2.378e-06 1.291e-07 0.5
MCMJSW 0.0005517 9.624e-05 0.5
LJSW700 0.003999 0.003125 0.5
meanLJSW 0.1387 0.09839 0.5
AGE 0.04508 0.002814 0.5
STDJSW 0.003999 0.002485 0.5
BMI 0.5 1 0.5
LatIntersept 0.4827 0.4091 0.5
CVJSW 0.5 1 0.5
MaxCV 0.1808 0.04675 0.5

Evaluating the OONR Models


boxplot(OAModelKL3$cvObject$Models.testPrediction[,"Prediction"]~OAModelKL3$cvObject$Models.testPrediction[,"Outcome"],main="CV Test Result",ylab="CV Test Prediction",xlab="KL Score")

KendallTable$CVOONRBSWIMS <- cor.test(OAModelKL3$cvObject$Models.testPrediction[,"Prediction"],OAModelKL3$cvObject$Models.testPrediction[,"Outcome"],method="kendall")
KappaTable$CVOONRBSWIMS <- epi.kappa(table(OAModelKL3$cvObject$Models.testPrediction[,"Prediction"]>2,OAModelKL3$cvObject$Models.testPrediction[,"Outcome"]>2.0))
RMSETable$CVOONRBSWIMS <- sqrt(mean((OAModelKL3$cvObject$Models.testPrediction[,"Prediction"]-OAModelKL3$cvObject$Models.testPrediction[,"Outcome"])^2))

bswimpredict <- predict(OAModelKL3$BSWiMS.model,testSubset)
boxplot(bswimpredict~testSubset$BLKL,main="B:SWiMS Test Result",ylab="Test Prediction",xlab="KL Score")

KendallTable$OONRBSWIMS <- cor.test(bswimpredict,testSubset$BLKL,method="kendall")
KappaTable$OONRBSWIMS <- epi.kappa(table(bswimpredict>2,testSubset$BLKL>2))
RMSETable$OONRBSWIMS <- sqrt(mean((bswimpredict-testSubset$BLKL)^2))


  
bswimpredict <- predict(OAModelKL3$cvObject$Fullenet,as.matrix(testSubset[,OAModelKL3$cvObject$LassoFilterVarList]))
boxplot(bswimpredict~testSubset$BLKL,main="LASSO Test Result",ylab="Test Prediction",xlab="KL Score")

KendallTable$OONRLASSO <- cor.test(bswimpredict,testSubset$BLKL,method="kendall")
KappaTable$OONRLASSO <- epi.kappa(table(bswimpredict>2,testSubset$BLKL>2))
RMSETable$OONRLASSO <- sqrt(mean((bswimpredict-testSubset$BLKL)^2))


bswimpredict <- ensemblePredict(OAModelKL3$eBSWiMS.model$formula.list,trainSubset,testSubset,Outcome ="BLKL",predictType ="linear",type = "LM",nk=-1 )
boxplot(bswimpredict$ensemblePredict~testSubset$BLKL,main="eB:SWiMS Test Result",ylab="Test Prediction",xlab="KL Score")

KendallTable$OONReBSWIMS <- cor.test(bswimpredict$ensemblePredict,testSubset$BLKL,method="kendall")
KappaTable$OONReBSWIMS <- epi.kappa(table(bswimpredict$ensemblePredict>2,testSubset$BLKL>2))
RMSETable$OONReBSWIMS <- sqrt(mean((bswimpredict$ensemblePredict-testSubset$BLKL)^2))


bswimpredict <- ensemblePredict(OAModelKL3$cvObject$formula.list,trainSubset,testSubset,Outcome ="BLKL",predictType ="linear",type = "LM",nk=-1 )
boxplot(bswimpredict$ensemblePredict~testSubset$BLKL,main="Ensemble B:SWiMS Test Result",ylab="Test Prediction",xlab="KL Score")

KendallTable$OONREnBSWIMS <- cor.test(bswimpredict$ensemblePredict,testSubset$BLKL,method="kendall")
KappaTable$OONREnBSWIMS <- epi.kappa(table(bswimpredict$ensemblePredict>2,testSubset$BLKL>2))
RMSETable$OONREnBSWIMS <- sqrt(mean((bswimpredict$ensemblePredict-testSubset$BLKL)^2))


Kendalls <- rbind(c(KendallTable$CVRawBSWIMS$estimate,KendallTable$LASSO$estimate,KendallTable$RawBSWIMS$estimate,KendallTable$raweBSWIMS$estimate,KendallTable$rawEnBSWIMS$estimate),
                  c(KendallTable$CVOONRBSWIMS$estimate,KendallTable$OONRLASSO$estimate,KendallTable$OONRBSWIMS$estimate,KendallTable$OONReBSWIMS$estimate,KendallTable$OONREnBSWIMS$estimate))
colnames(Kendalls) <- c("CV Estimation","LASSO","B:SWiMS","eB:SWiMS","Ensemble B:SWiMS")
barplot(Kendalls,ylab="Correlation",main="Kendall Correlation by Prediction Method",beside=TRUE,legend =c("Raw","Categorical"),ylim=c(0,1.0),args.legend = list(x = "topright"),cex.names=0.7,las=2)


RMSEt <- rbind(c(RMSETable$CVRawBSWIMS,RMSETable$LASSO,RMSETable$RawBSWIMS,RMSETable$raweBSWIMS,RMSETable$rawEnBSWIMS),c(RMSETable$CVOONRBSWIMS,RMSETable$OONRLASSO,RMSETable$OONRBSWIMS,RMSETable$OONReBSWIMS,RMSETable$OONREnBSWIMS))
colnames(RMSEt) <- c("CV Estimation","LASSO","B:SWiMS","eB:SWiMS","Ensemble B:SWiMS")
barplot(RMSEt,ylab="RMSE",main="RMSE by Prediction Method",beside=TRUE,legend =c("Raw","Categorical"),ylim=c(0,1.2),args.legend = list(x = "topright"),cex.names=0.7,las=2)


Kappacis <- as.matrix(rbind(
                   KappaTable$CVRawBSWIMS$kappa[-2],
                   KappaTable$LASSO$kappa[-2],
                   KappaTable$RawBSWIMS$kappa[-2],
                   KappaTable$raweBSWIMS$kappa[-2],
                   KappaTable$rawEnBSWIMS$kappa[-2],
                   KappaTable$CVOONRBSWIMS$kappa[-2],
                   KappaTable$OONRLASSO$kappa[-2],
                   KappaTable$OONRBSWIMS$kappa[-2],
                   KappaTable$OONReBSWIMS$kappa[-2],
                   KappaTable$OONREnBSWIMS$kappa[-2]
                   ))

bplot <- barPlotCiError(Kappacis,"Kappa",c("CV","LASSO","B:SWiMS","eB:SWiMS","En.B:SWiMS"),c("Raw","Categorical"),main="Kappa Agreement",args.legend = list(x = "bottomright"))

Study the longitudinal Beheavior of the predicted KLScores

par(cex=0.75)
suana <- subset(OAI_JSW.full.adj.norm,BLKL<4)
#suana[,-c(1:8)] <- nearestneighborimpute(suana[,-c(1:8)])
suana$KLG0 <- 1*(suana$BLKL>0)

timevar <- rbind(c("CVJSW","Coefficient of Variation JSW"),c("KLPredict","Raw KL Predict"))
colnames(timevar) <- c("Var","FullDescription")
timevar <- rbind(timevar,c("CKLPredict","CKLPredict"),c("LINKLPredict","LINKLPredict"))
timevar <- rbind(timevar,c("LIN2KLPredict","LIN2KLPredict"))



pr <- predict(OAModelKLOM$oridinalModels,suana)
summary(pr)
suana$KLPredict <- pr[,1]
suana$CKLPredict <- pr[,2]
suana$LINKLPredict <- predict(OAModelKLOM$bagging$bagged.model,suana)
suana$LIN2KLPredict <- predict(OAModelKL$bagging$bagged.model,suana)

tsa <- timeSerieAnalysis(timevar,"I(YearsFromBL-2)*KLG0",suana,timevar="VISITID",contime="YearsFromBL",Outcome="KLG0", cs1CorCAR1,timesign= "+",Ptoshow = c(2,3,4),plegend= c("Visit","KLG0","Visit:KLG0"),catgo.names=c("KL=0", "KL>0"),description="FullDescription") 



suana$LASSOKLPredict <- predict(OAModelKL$cvObject$Fullenet,as.matrix(suana[,OAModelKL$cvObject$LassoFilterVarList]))
timevar <- rbind(timevar,c("LASSOKLPredict","LASSO KL"))
suana$CatKLPredict <- predict(OAModelKL3$BSWiMS.model,suana)
timevar <- rbind(timevar,c("CatKLPredict","KL Categorical Model"))

Modeling the TKR event with Logitic Models

Selecting the train and test sets

baseline <- subset(OAI_JSW.full.adj.norm,(VISITID==0) & (BLKL>1) & (TimeToEvent<1825 | TimeToEvent>3200) & (DeathVisit==0))
rownames(baseline) <- baseline$IDSIDE
baseline <- baseline[complete.cases(baseline),c("TKREvent","BLKL","AGE","BMI","Gender","HEIGHT","WEIGHT","AnySurgery","WOMACTOTAL","WOMACPAIN","PainMedication","mJSN","lJSN",as.character(OAI_JSWVarDescription$Var))]
TKR0 <- subset(baseline,TKREvent == 0)
TKR1 <- subset(baseline,TKREvent == 1)
trainIDs <- sample(nrow(TKR1),nrow(TKR1)/2)
trainIDCs <- sample(nrow(TKR0),nrow(TKR0)/10)
trainSubset <- rbind(TKR0[trainIDCs,],TKR1[trainIDs,])
testSubset <- rbind(TKR0[-trainIDCs,],TKR1[-trainIDs,])

Modling with normalized features

system.time(OAModelTKR <- FRESA.Model(TKREvent ~ 1,trainSubset,CVfolds=5,repeats=5,categorizationType="Raw",equivalent=TRUE,usrFitFun=rpart))
BSModelTKR <- BSWiMS.model(TKREvent ~ 1,trainSubset)
OAModelTKR$BSWiMS.model$coefficients
OAModelTKR$BSWiMS.models$formula.list
OAModelTKR$forwardModel$final.model$coefficients
#OAModelTKR$eBSWiMS.model$equivalentMatrix
#summary(OAModelTKR$BSWiMS.model)
#am <- modelFitting(OAModelTKR$BSWiMS.model$formula,trainSubset,type="LOGIT",fast=TRUE)
#summary.fitFRESA(am)
#summary.fitFRESA(am,type="Residual")

Analyzing the CV results

plot(OAModelTKR$bootstrappedModel)



pmr <- plotModels.ROC(OAModelTKR$cvObject$Models.testPrediction,main="FRESA.CAD Models",cex=0.80)

pmr <- plotModels.ROC(OAModelTKR$cvObject$Models.testPrediction,theCVfolds=5,predictor="Ensemble.Forward",main="Ensemble",cex=0.80)

pmr <- plotModels.ROC(OAModelTKR$cvObject$Models.testPrediction,theCVfolds=5,predictor="Forward.Selection.Bagged",main="Bagging",cex=0.80)

pmr <- plotModels.ROC(OAModelTKR$cvObject$Models.testPrediction,theCVfolds=5,predictor="Forward",main="Forward",cex=0.80)

pmr <- plotModels.ROC(OAModelTKR$cvObject$Models.testPrediction,theCVfolds=5,predictor="Backwards",main="Backwards",cex=0.80)

pmr <- plotModels.ROC(OAModelTKR$cvObject$Models.testPrediction,theCVfolds=5,predictor="Prediction",main="B:SWiMS",cex=0.80)

pmr <- plotModels.ROC(OAModelTKR$cvObject$Models.testPrediction,theCVfolds=5,predictor="eB.SWiMS",main="eB:SWIMS",cex=0.80)

pmr <- plotModels.ROC(OAModelTKR$cvObject$Models.testPrediction,theCVfolds=5,predictor="Ensemble.B.SWiMS",main="eB:SWIMS",cex=0.80)

pmr <- plotModels.ROC(OAModelTKR$cvObject$Models.testPrediction,theCVfolds=5,predictor="usrFitFunction",main="rpart",cex=0.80)

Evaluating the Model performance on test



ROCTable <- NULL;
epiTable <- NULL;
ROCMedians <- NULL;


bswimpredict <- as.vector(predict(OAModelTKR$BSWiMS.model,testSubset))
pm<-plotModels.ROC(cbind(testSubset$TKREvent,bswimpredict),main="B:SWiMS")

ROCTable$rawBSWIMS  <- pm$roc.predictor
epiTable$rawBSWIMS <- epi.tests(pm$predictionTable)

rpt <- rpart(TKREvent ~ .,trainSubset)
rpre <- predict(rpt,testSubset)
pm <- plotModels.ROC(cbind(testSubset$TKREvent,rpre),main="rpart",cex=0.90)

ROCTable$rawRPART  <- pm$roc.predictor
epiTable$rawRPART <- epi.tests(pm$predictionTable)


tkr.rf <- randomForest(TKREvent ~ .,trainSubset)
rf.pred <- predict(tkr.rf,testSubset)-0.5
pm <- plotModels.ROC(cbind(testSubset$TKREvent,rf.pred),main="rpart",cex=0.90)

ROCTable$rawRandFor  <- pm$roc.predictor
epiTable$rawRandFor <- epi.tests(pm$predictionTable)

bswimpredict <- as.vector(predict(OAModelTKR$cvObject$Fullenet,as.matrix(testSubset[,OAModelTKR$cvObject$LassoFilterVarList])))
pm <- plotModels.ROC(cbind(testSubset$TKREvent,bswimpredict),main="LASSO",cex=0.90)

ROCTable$rawLASSO  <- pm$roc.predictor
epiTable$rawLASSO <- epi.tests(pm$predictionTable)

bswimpredict <- ensemblePredict(OAModelTKR$eBSWiMS.model$formula.list,trainSubset,testSubset,Outcome ="TKREvent",predictType ="linear",type = "LOGIT")
colnames(bswimpredict$predictions)
ROCMedians$raweBSWIMS  <- plotModels.ROC(bswimpredict$predictions,main="eB:SWiMS",cex=0.90)

pm <- plotModels.ROC(cbind(testSubset$TKREvent,bswimpredict$ensemblePredict),main="eB:SWiMS")

ROCTable$raweBSWIMS  <- pm$roc.predictor
epiTable$raweBSWIMS <- epi.tests(pm$predictionTable)

bswimpredict <- ensemblePredict(OAModelTKR$cvObject$formula.list,trainSubset,testSubset,Outcome ="TKREvent",predictType ="linear",type = "LOGIT")
ROCMedians$rawCVBSWIMS <- plotModels.ROC(bswimpredict$predictions,main="CV B:SWiMS",cex=0.90)

pm <- plotModels.ROC(cbind(testSubset$TKREvent,bswimpredict$ensemblePredict),main="CV B:SWiMS")

ROCTable$rawCVBSWIMS  <- pm$roc.predictor
epiTable$rawCVBSWIMS <- epi.tests(pm$predictionTable)

bswimpredict <- ensemblePredict(OAModelTKR$cvObject$LASSOVariables,trainSubset,testSubset,Outcome ="TKREvent",predictType ="linear",type = "LOGIT")
ROCMedians$rawCVLASSO <- plotModels.ROC(bswimpredict$predictions,main="CV LASSO",cex=0.90)

pm <- plotModels.ROC(cbind(testSubset$TKREvent,bswimpredict$ensemblePredict),main="CV LASSO")

ROCTable$rawCVLASSO  <- pm$roc.predictor
epiTable$rawCVLASSO <- epi.tests(pm$predictionTable)

The equivalent matrix

pander::pander(OAModelTKR$eBSWiMS.model$equivalentMatrix)
Table continues below
Name Locus Extendend_Name
WOMACTOTAL PainMedication WOMACTOTAL:PainMedication
WOMACPAIN PainMedication WOMACPAIN:PainMedication
PainMedication PainMedication PainMedication:PainMedication
STDJSW BLKL STDJSW:BLKL
CVJSW BLKL CVJSW:BLKL
WOMACTOTAL BLKL WOMACTOTAL:BLKL
WOMACPAIN BLKL WOMACPAIN:BLKL
MaxCV BLKL MaxCV:BLKL
MtoLDiff BLKL MtoLDiff:BLKL
DiffSlope BLKL DiffSlope:BLKL
JSW300 BLKL JSW300:BLKL
LatSlope BLKL LatSlope:BLKL
mJSN BLKL mJSN:BLKL
lJSN BLKL lJSN:BLKL
AnySurgery BLKL AnySurgery:BLKL
BLKL BLKL BLKL:BLKL
LJSW700 BLKL LJSW700:BLKL
MeanJSW BLKL MeanJSW:BLKL
lCV BLKL lCV:BLKL
LJSW725 BLKL LJSW725:BLKL
meanMJSW BLKL meanMJSW:BLKL
LJSW750 BLKL LJSW750:BLKL
LateralSlope BLKL LateralSlope:BLKL
mCV BLKL mCV:BLKL
JSW250 BLKL JSW250:BLKL
JSW200 BLKL JSW200:BLKL
JSW175 BLKL JSW175:BLKL
Table continues below
UniPerformance FullPerformance DeltaPerformance ImprovementFraction
0.5342 0.6501 0.01797 0.3385
0.5379 0.6248 -0.007379 0.3531
0.5 0.6916 0.05947 0.3781
0.5977 0.6465 0.1465 0.306
0.6303 0.6483 0.1483 0.3571
0.5342 0.5867 0.08672 0.2732
0.5379 0.5378 0.03781 0.208
0.5597 0.5977 0.09767 0.2195
0.6068 0.6574 0.1574 0.2951
0.5181 0.6067 0.1067 0.1857
0.5561 0.5958 0.0958 0.2075
0.5144 0.6085 0.1085 0.1891
0.596 0.6428 0.1428 0.2517
0.5 0.5488 0.04883 0.04578
0.5 0.5995 0.09945 0.1799
0.6321 0.6916 0.1916 0.3489
0.5199 0.5977 0.09767 0.1963
0.5199 0.6175 0.1175 0.2037
0.5199 0.5507 0.0507 0.05228
0.5145 0.5579 0.05788 0.1385
0.5543 0.6175 0.1175 0.1204
0.5054 0.5397 0.03975 0.1313
0.5127 0.5995 0.09949 0.1565
0.5253 0.5614 0.06145 0.09172
0.5543 0.6175 0.1175 0.1423
0.5398 0.6121 0.1121 0.2003
0.5235 0.6067 0.1067 0.1929
p.value
7.938e-06
5.001e-05
8.716e-08
1.693e-07
6.986e-09
0.001417
0.00283
0.0002893
1.177e-07
0.002946
0.000296
0.001859
0.0004609
0.01979
0.001796
7.859e-10
0.00079
6.299e-05
0.01858
0.00242
0.0004949
0.008086
0.005137
0.01217
0.0003735
0.0006663
0.001532

Heat map of the features


fnames <- unique(as.character(OAModelTKR$eBSWiMS.model$equivalentMatrix$Name))
topFeatures <- OAI_JSWVarDescription[fnames,]
topFeatures <- topFeatures[!is.na(topFeatures$Var),]
if (length(topFeatures[,1])>3)
{
hmdata <- baseline[sample(nrow(baseline),250),]


featureColor <- as.character(topFeatures[,"Type.Color"])
lJSNBar <- hmdata$lJSN
names(lJSNBar) <- rownames(hmdata)
mJSNBar <- hmdata$mJSN
names(mJSNBar) <- rownames(hmdata)
ConData <- hmdata$WOMACTOTAL
ConDataName <- "Total WOMAC"
names(ConData) <- rownames(hmdata)
PainMedication <- hmdata$PainMedication
names(PainMedication) <- rownames(hmdata)
hm <- NULL;


op <- par(no.readonly = TRUE)
hm <- heatMaps(variableList = topFeatures,Outcome = "TKREvent",data = hmdata,hCluster= FALSE,theFiveColors=colorpalette,outcomeColors=colorJSNpalette,transpose=TRUE,xlab="Baseline",ylab="Features",cexRow=0.50,cexCol=0.10,srtCol=35,reorderfun=function(d, w) reorder(d, w, agglo.FUN = mean),distfun=function(c) dist(c, method = "manhattan"),keysize=1.0,RowSideColors=as.character(topFeatures[,"Region.Color"]),lmat=lmat, lhei=lhei, lwid=lwid,extrafun=myplot,key.par=list(cex=0.5))
par(new=TRUE,mar=c(1,2,0,0),omd=c(0.02,0.98,0.02,0.98))
plot(c(0,5), c(0,1), type='n', bty='n', xaxt='n', xlab='', yaxt='n', ylab='',main='')
lg <- legend("bottomleft",legend=feattype,fill=featcolor,cex=0.5)
legend(lg$rect$left+lg$rect$w,lg$rect$top,legend=lType,fill=lTypecolors,cex=0.5)
par(op)
}

Modeling TKR with out-of-normal-range Features

system.time(OAModelTKR2 <- FRESA.Model(TKREvent ~ 1,trainSubset,CVfolds=5,repeats=5,categorizationType="RawZTail",equivalent=TRUE))

Analyzing the with out-of-normal-range Model


bswimpredict <- as.vector(predict(OAModelTKR2$BSWiMS.model,testSubset))
pm <- plotModels.ROC(cbind(testSubset$TKREvent,bswimpredict),main="B:SWiMS")

ROCTable$OONRBSWIMS  <- pm$roc.predictor
epiTable$OONRBSWIMS <- epi.tests(pm$predictionTable)

bswimpredict <- as.vector(predict(OAModelTKR2$cvObject$Fullenet,as.matrix(testSubset[,OAModelTKR2$cvObject$LassoFilterVarList])))
pm <- plotModels.ROC(cbind(testSubset$TKREvent,bswimpredict),main="LASSO")

ROCTable$OONRLASSO  <- pm$roc.predictor
epiTable$OONRLASSO <- epi.tests(pm$predictionTable)

bswimpredict <- ensemblePredict(OAModelTKR2$eBSWiMS.model$formula.list,trainSubset,testSubset,Outcome ="TKREvent",predictType ="linear",type = "LOGIT")
ROCMedians$OONReBSWIMS  <- plotModels.ROC(bswimpredict$predictions,main="eB:SWiMS",cex=0.80)

pm <- plotModels.ROC(cbind(testSubset$TKREvent,bswimpredict$ensemblePredict),main="eB:SWiMS")

ROCTable$OONReBSWIMS  <- pm$roc.predictor
epiTable$OONReBSWIMS <- epi.tests(pm$predictionTable)

bswimpredict <- ensemblePredict(OAModelTKR2$cvObject$formula.list,trainSubset,testSubset,Outcome ="TKREvent",predictType ="linear",type = "LOGIT")
ROCMedians$OONRCVBSWIMS <- plotModels.ROC(bswimpredict$predictions,main="CV B:SWiMS",cex=0.80)

pm <- plotModels.ROC(cbind(testSubset$TKREvent,bswimpredict$ensemblePredict),main="CV B:SWiMS")

ROCTable$OONRCVBSWIMS  <- pm$roc.predictor
epiTable$OONRCVBSWIMS <- epi.tests(pm$predictionTable)

bswimpredict <- ensemblePredict(OAModelTKR2$cvObject$LASSOVariables,trainSubset,testSubset,Outcome ="TKREvent",predictType ="linear",type = "LOGIT")
ROCMedians$OONRCVLASSO <- plotModels.ROC(bswimpredict$predictions,main="CV LASSO",cex=0.80)

pm <-  plotModels.ROC(cbind(testSubset$TKREvent,bswimpredict$ensemblePredict),main="CV LASSO")

ROCTable$OONRCVLASSO  <- pm$roc.predictor
epiTable$OONRCVLASSO <- epi.tests(pm$predictionTable)

pander::pander(OAModelTKR2$eBSWiMS.model$equivalentMatrix)
Table continues below
Name Locus
WOMACPAIN PainMedication
WOMACTOTAL PainMedication
PainMedication PainMedication
I(MtoLDiff * (MtoLDiff > 1.282)) BLKL
I(LateralSlope * (LateralSlope > 1.282)) BLKL
WOMACPAIN BLKL
WOMACTOTAL BLKL
I(CVJSW * (CVJSW > 1.282)) BLKL
CVJSW BLKL
I(LJSW725 * (LJSW725 < -1.282)) BLKL
I(MeanJSW * (MeanJSW < -1.282)) BLKL
I(STDJSW * (STDJSW > 1.282)) BLKL
STDJSW BLKL
BLKL BLKL
I(LatSlope * (LatSlope > 1.282)) BLKL
MtoLDiff BLKL
I(JSW275 * (JSW275 < -1.282)) BLKL
I(MaxCV * (MaxCV > 1.282)) BLKL
I(DiffSlope * (DiffSlope < -1.282)) BLKL
I(LJSW700 * (LJSW700 < -1.282)) BLKL
I(JSW300 * (JSW300 < -1.282)) BLKL
I(meanLJSW * (meanLJSW < -1.282)) BLKL
I(MCMJSW * (MCMJSW < -1.282)) BLKL
LJSW700 BLKL
I(JSW225 * (JSW225 < -1.282)) BLKL
I(LJSW775 * (LJSW775 < -1.282)) BLKL
LJSW725 BLKL
MeanJSW BLKL
I(MedIntersept * (MedIntersept < -1.282)) BLKL
I(lCV * (lCV > 1.282)) BLKL
JSW300 BLKL
mJSN BLKL
MCMJSW BLKL
I(LJSW850 * (LJSW850 < -1.282)) BLKL
DiffSlope BLKL
AnySurgery BLKL
LatSlope BLKL
I(JSW150 * (JSW150 < -1.282)) BLKL
MaxCV BLKL
mCV BLKL
Table continues below
Extendend_Name UniPerformance FullPerformance
WOMACPAIN:PainMedication 0.5379 0.6248
WOMACTOTAL:PainMedication 0.5342 0.6501
PainMedication:PainMedication 0.5 0.6916
I(MtoLDiff * (MtoLDiff > 1.282)):BLKL 0.6049 0.6574
I(LateralSlope * (LateralSlope > 1.282)):BLKL 0.5434 0.6175
WOMACPAIN:BLKL 0.5379 0.5378
WOMACTOTAL:BLKL 0.5342 0.5867
I(CVJSW * (CVJSW > 1.282)):BLKL 0.6321 0.6465
CVJSW:BLKL 0.6303 0.6483
I(LJSW725 * (LJSW725 < -1.282)):BLKL 0.5525 0.6031
I(MeanJSW * (MeanJSW < -1.282)):BLKL 0.5543 0.6211
I(STDJSW * (STDJSW > 1.282)):BLKL 0.5977 0.6447
STDJSW:BLKL 0.5977 0.6465
BLKL:BLKL 0.6321 0.6916
I(LatSlope * (LatSlope > 1.282)):BLKL 0.5561 0.6175
MtoLDiff:BLKL 0.6068 0.6574
I(JSW275 * (JSW275 < -1.282)):BLKL 0.5996 0.6175
I(MaxCV * (MaxCV > 1.282)):BLKL 0.5706 0.5977
I(DiffSlope * (DiffSlope < -1.282)):BLKL 0.5307 0.6067
I(LJSW700 * (LJSW700 < -1.282)):BLKL 0.5416 0.6393
I(JSW300 * (JSW300 < -1.282)):BLKL 0.5814 0.6193
I(meanLJSW * (meanLJSW < -1.282)):BLKL 0.5235 0.5632
I(MCMJSW * (MCMJSW < -1.282)):BLKL 0.5742 0.6139
LJSW700:BLKL 0.5199 0.5977
I(JSW225 * (JSW225 < -1.282)):BLKL 0.596 0.6229
I(LJSW775 * (LJSW775 < -1.282)):BLKL 0.5127 0.5651
LJSW725:BLKL 0.5145 0.5579
MeanJSW:BLKL 0.5199 0.6175
I(MedIntersept * (MedIntersept < -1.282)):BLKL 0.596 0.6102
I(lCV * (lCV > 1.282)):BLKL 0.5235 0.5597
JSW300:BLKL 0.5561 0.5958
mJSN:BLKL 0.596 0.6428
MCMJSW:BLKL 0.5271 0.5958
I(LJSW850 * (LJSW850 < -1.282)):BLKL 0.5217 0.5633
DiffSlope:BLKL 0.5181 0.6067
AnySurgery:BLKL 0.5 0.5995
LatSlope:BLKL 0.5144 0.6085
I(JSW150 * (JSW150 < -1.282)):BLKL 0.5724 0.6139
MaxCV:BLKL 0.5597 0.5977
mCV:BLKL 0.5253 0.5614
DeltaPerformance ImprovementFraction p.value
-0.007379 0.3531 5.001e-05
0.01797 0.3385 7.938e-06
0.05947 0.3781 8.716e-08
0.1574 0.2987 1.272e-07
0.1175 0.2487 8.82e-05
0.03781 0.208 0.00283
0.08672 0.2732 0.001417
0.1465 0.3535 1.224e-08
0.1483 0.3571 6.986e-09
0.1031 0.1795 0.0002208
0.1211 0.2333 1.176e-05
0.1447 0.3097 3.468e-07
0.1465 0.306 1.693e-07
0.1916 0.3489 7.859e-10
0.1175 0.2449 0.0001278
0.1574 0.2951 1.177e-07
0.1175 0.215 2.306e-05
0.09767 0.2085 0.000271
0.1067 0.2126 0.0002918
0.1393 0.2375 3.505e-05
0.1193 0.2299 4.158e-05
0.06323 0.1105 0.005045
0.1139 0.1857 0.0001472
0.09767 0.1963 0.00079
0.1229 0.1968 4.704e-05
0.06506 0.1106 0.004372
0.05788 0.1385 0.00242
0.1175 0.2037 6.299e-05
0.1102 0.1788 4.747e-05
0.05974 0.1078 0.01011
0.0958 0.2075 0.000296
0.1428 0.2517 0.0004609
0.0958 0.1818 0.001219
0.06327 0.09616 0.0123
0.1067 0.1857 0.002946
0.09945 0.1799 0.001796
0.1085 0.1891 0.001859
0.1139 0.1857 0.0001863
0.09767 0.2195 0.0002893
0.06145 0.09172 0.01217 

acctables <- as.matrix(rbind(
                   epiTable$rawLASSO$elements$diag.acc,
                   epiTable$rawCVLASSO$elements$diag.acc,
                   epiTable$rawBSWIMS$elements$diag.acc,
                   epiTable$raweBSWIMS$elements$diag.acc,
                   epiTable$rawCVBSWIMS$elements$diag.acc,
                   epiTable$OONRLASSO$elements$diag.acc,
                   epiTable$OONRCVLASSO$elements$diag.acc,
                   epiTable$OONRBSWIMS$elements$diag.acc,
                   epiTable$OONReBSWIMS$elements$diag.acc,
                   epiTable$OONRCVBSWIMS$elements$diag.acc
                   ))

bplot <- barPlotCiError(acctables,"Accuracy",c("LASSO","CV LASSO","B:SWiMS","eB:SWiMS","CV B:SWiMS"),c("Standard","Categorical"),main="Prediction Accuracy",args.legend = list(x = "bottomright"))


sentables <- as.matrix(rbind(
                   epiTable$rawLASSO$elements$sensitivity,
                   epiTable$rawCVLASSO$elements$sensitivity,
                   epiTable$rawBSWIMS$elements$sensitivity,
                   epiTable$raweBSWIMS$elements$sensitivity,
                   epiTable$rawCVBSWIMS$elements$sensitivity,
                   epiTable$OONRLASSO$elements$sensitivity,
                   epiTable$OONRCVLASSO$elements$sensitivity,
                   epiTable$OONRBSWIMS$elements$sensitivity,
                   epiTable$OONReBSWIMS$elements$sensitivity,
                   epiTable$OONRCVBSWIMS$elements$sensitivity
                   ))

bplot <- barPlotCiError(sentables,"Sensitivity",c("LASSO","CV LASSO","B:SWiMS","eB:SWiMS","CV B:SWiMS"),c("Standard","Categorical"),main="Prediction Sensitivity",args.legend = list(x = "bottomright"))



errtables <- 1.0-as.matrix(0.5*rbind(
                   epiTable$rawLASSO$elements$sensitivity+epiTable$rawLASSO$elements$specificity,
                   epiTable$rawCVLASSO$elements$sensitivity+epiTable$rawCVLASSO$elements$specificity,
                   epiTable$rawBSWIMS$elements$sensitivity+epiTable$rawBSWIMS$elements$specificity,
                   epiTable$raweBSWIMS$elements$sensitivity+epiTable$raweBSWIMS$elements$specificity,
                   epiTable$rawCVBSWIMS$elements$sensitivity+epiTable$rawCVBSWIMS$elements$specificity,
                   epiTable$OONRLASSO$elements$sensitivity+epiTable$OONRLASSO$elements$specificity,
                   epiTable$OONRCVLASSO$elements$sensitivity+epiTable$OONRCVLASSO$elements$specificity,
                   epiTable$OONRBSWIMS$elements$sensitivity+epiTable$OONRBSWIMS$elements$specificity,
                   epiTable$OONReBSWIMS$elements$sensitivity+epiTable$OONReBSWIMS$elements$specificity,
                   epiTable$OONRCVBSWIMS$elements$sensitivity+epiTable$OONRCVBSWIMS$elements$specificity
                   ))

bplot <- barPlotCiError(errtables,"Balanced Error",c("LASSO","CV LASSO","B:SWiMS","eB:SWiMS","CV B:SWiMS"),c("Standard","Categorical"),main="Prediction Error",args.legend = list(x = "topright"))



auctables <- rbind(
                   ci(ROCTable$rawLASSO,of="auc"),
                   ci(ROCTable$rawCVLASSO,of="auc"),
                   ci(ROCTable$rawBSWIMS,of="auc"),
                   ci(ROCTable$raweBSWIMS,of="auc"),
                   ci(ROCTable$rawCVBSWIMS,of="auc"),
                   ci(ROCTable$OONRLASSO,of="auc"),
                   ci(ROCTable$OONRCVLASSO,of="auc"),
                   ci(ROCTable$OONRBSWIMS,of="auc"),
                   ci(ROCTable$OONReBSWIMS,of="auc"),
                   ci(ROCTable$OONRCVBSWIMS,of="auc")
                   )
low <- auctables[,1]; auctables[,1] <- auctables[,2]; auctables[,2] <- low;
bplot <- barPlotCiError(auctables,"ROC AUC",c("LASSO","CV LASSO","B:SWiMS","eB:SWiMS","CV B:SWiMS"),c("Standard","Categorical"),main="Prediction AUC",args.legend = list(x = "topright"))