
library(BART)

dir.create('temp')
setwd('temp')
download.file('https://www.mcw.edu/-/media/MCW/Departments/Biostatistics/tr067.zip',
              'tr067.zip')
unzip('tr067.zip')
##unzip('../tr067.zip')

table.=function(y, x, dnn) {
    a=table(y, x, dnn=dnn)    
    print(a)
    spec=round(100*a[1, 1]/sum(a[1, ]), 1) 
    sens=round(100*a[2, 2]/sum(a[2, ]), 1)
    print(c('Spec.'=spec, 'Sens.'=sens))
}

train=readRDS('mesa.train.rds')
test=readRDS('mesa.test.rds')

bart26=readRDS('bart26.rds')
bart26.names=names(bart26$varcount.mean)
print(bart26.names)

x.train=as.matrix(train)[ , bart26.names]
str(x.train)

x.test=as.matrix(test)[ , bart26.names]
str(x.test)

yhat.train=predict(bart26, x.train)
yhat.train.mean=apply(yhat.train, 2, mean)

print(c('R square: Normalized LVM by cMRI with BART26'=
            cor(train$lvh_lvm, yhat.train.mean)^2))

plot(train$lvh_lvm, yhat.train.mean, asp=1,
     xlab='Normalized LVM measured by cMRI', pch='.',
     ylab='Normalized LVM predicted by BART26',
     sub='training data points',
     col=2*train$gender1+2)
dev.copy2pdf(file='NormLVM-train-BART26.pdf')

print(c('R square: LVM by cMRI with BART26'=cor(train$olvedm1,
                          yhat.train.mean*train$exp_lvm)^2))

plot(train$olvedm1, yhat.train.mean*train$exp_lvm, asp=1,
     xlab='LVM measured by cMRI', pch='.',
     ylab='LVM predicted by BART26',
     sub='training data points',
     col=2*train$gender1+2)
dev.copy2pdf(file='LVM-train-BART26.pdf')

table.(train$lvh_lvm>=1.31, yhat.train.mean>=1.19,
       dnn=list('MRI LVH', 'BART26 training'))

yhat.test=predict(bart26, x.test)
yhat.test.mean=apply(yhat.test, 2, mean)

print(c('R square: Normalized LVM by cMRI with BART26'=
            cor(test$lvh_lvm, yhat.test.mean)^2))

plot(test$lvh_lvm, yhat.test.mean, asp=1,
     xlab='Normalized LVM measured by cMRI', pch='.',
     ylab='Normalized LVM predicted by BART26',
     sub='validation data points',
     col=2*test$gender1+2)
dev.copy2pdf(file='NormLVM-test-BART26.pdf')

print(c('R square: LVM by cMRI with BART26'=cor(test$olvedm1,
                          yhat.test.mean*test$exp_lvm)^2))

plot(test$olvedm1, yhat.test.mean*test$exp_lvm, asp=1,
     xlab='LVM measured by cMRI', pch='.',
     ylab='LVM predicted by BART26',
     sub='validation data points',
     col=2*test$gender1+2)
dev.copy2pdf(file='LVM-test-BART26.pdf')

table.(test$lvh_lvm>=1.31, yhat.test.mean>=1.19,
       dnn=list('MRI LVH', 'BART26 validation'))

bart22=readRDS('bart22.rds')
bart22.names=names(bart22$varcount.mean)
print(bart22.names)

x.train=as.matrix(train)[ , bart22.names]
str(x.train)

x.test=as.matrix(test)[ , bart22.names]
str(x.test)

yhat.train=predict(bart22, x.train)
yhat.train.mean=apply(yhat.train, 2, mean)

print(c('R square: Normalized LVM by cMRI with BART22'=
            cor(train$lvh_lvm, yhat.train.mean)^2))

plot(train$lvh_lvm, yhat.train.mean, asp=1,
     xlab='Normalized LVM measured by cMRI', pch='.',
     ylab='Normalized LVM predicted by BART22',
     sub='training data points',
     col=2*train$gender1+2)
dev.copy2pdf(file='NormLVM-train-BART22.pdf')

print(c('R square: LVM by cMRI with BART22'=cor(train$olvedm1,
                          yhat.train.mean*train$exp_lvm)^2))

plot(train$olvedm1, yhat.train.mean*train$exp_lvm, asp=1,
     xlab='LVM measured by cMRI', pch='.',
     ylab='LVM predicted by BART22',
     sub='training data points',
     col=2*train$gender1+2)
dev.copy2pdf(file='LVM-train-BART22.pdf')

table.(train$lvh_lvm>=1.31, yhat.train.mean>=1.19,
       dnn=list('MRI LVH', 'BART22 training'))

yhat.test=predict(bart22, x.test)
yhat.test.mean=apply(yhat.test, 2, mean)

print(c('R square: Normalized LVM by cMRI with BART22'=
            cor(test$lvh_lvm, yhat.test.mean)^2))

plot(test$lvh_lvm, yhat.test.mean, asp=1,
     xlab='Normalized LVM measured by cMRI', pch='.',
     ylab='Normalized LVM predicted by BART22',
     sub='validation data points',
     col=2*test$gender1+2)
dev.copy2pdf(file='NormLVM-test-BART22.pdf')

print(c('R square: LVM by cMRI with BART22'=cor(test$olvedm1,
                          yhat.test.mean*test$exp_lvm)^2))

plot(test$olvedm1, yhat.test.mean*test$exp_lvm, asp=1,
     xlab='LVM measured by cMRI', pch='.',
     ylab='LVM predicted by BART22',
     sub='validation data points',
     col=2*test$gender1+2)
dev.copy2pdf(file='LVM-test-BART22.pdf')

table.(test$lvh_lvm>=1.31, yhat.test.mean>=1.19,
       dnn=list('MRI LVH', 'BART22 training'))
