Homework #3: Case Studies
Medical Physics (Biophys430)
Kyle Garcia-Rogers
March 2024
Problem #1
# Download the 4D fMRI hyper-volume (3D space * 1D time)
# library(brainR) # for NIFTI file input
library(fmri) # for NIFTI file input
library(RNifti)
library(oro.nifti)## oro.nifti 0.11.4##
## Attaching package: 'oro.nifti'## The following object is masked from 'package:RNifti':
##
## originfMRIURL <- "https://socr.umich.edu/HTML5/BrainViewer/data/fMRI_FilteredData_4D.nii.gz"
fMRIFile <- file.path(tempdir(), "fMRI_FilteredData_4D.nii.gz")
download.file(fMRIURL, dest=fMRIFile, quiet=TRUE)
(fMRIVolume <- readNIfTI(fMRIFile, reorient=FALSE))## NIfTI-1 format
## Type : nifti
## Data Type : 4 (INT16)
## Bits per Pixel : 16
## Slice Code : 0 (Unknown)
## Intent Code : 0 (None)
## Qform Code : 1 (Scanner_Anat)
## Sform Code : 0 (Unknown)
## Dimension : 64 x 64 x 21 x 180
## Pixel Dimension : 4 x 4 x 6 x 3
## Voxel Units : mm
## Time Units : sec## Loading required package: ggplot2##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:oro.nifti':
##
## slice## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutx1 <- c(1:180)
y1 <- loess(fMRIVolume[30, 30, 10,]~ x1, family = "gaussian")$y
y2 <- ksmooth(x1, fMRIVolume[30, 30, 10,], kernel = "normal", bandwidth = 5)$y
plot_ly(x=~x1, y=~fMRIVolume[30, 30, 10,], type="scatter", mode="lines",
name="Raw fMRI[30, 30, 10,t]") %>%
add_trace(y=~smooth(fMRIVolume[30, 30, 10,]), name="Smooth(fMRI)") %>%
add_trace(y=~y1, name="LOESS(fMRI)") %>%
add_trace(y=~y2, name="ksmooth(fMRI)") %>%
layout(xaxis=list(title="time"), yaxis=list(title="intensity"), legend = list(orientation='h'))# Define a pair of fMRI timeseries at 2 separate voxel locations
y_30_30_10 <- fMRIVolume[30, 30, 10,]
y_20_30_15 <- fMRIVolume[20, 30, 15,]
# Examine the top-6 of the pair of intensities
head(cbind(y_30_30_10, y_20_30_15))## y_30_30_10 y_20_30_15
## [1,] 14140 6278
## [2,] 14038 6294
## [3,] 13888 6326
## [4,] 13989 6355
## [5,] 14015 6341
## [6,] 14035 6288# Plot y_30_30_10 (x-axis) against y_20_30_15 (y-axis)
#plot(y_30_30_10, y_20_30_15)
plot_ly(x=~y_30_30_10, y=~y_20_30_15, type="scatter", mode="markers")## [1] 0.01868463# Compute the Auto-correlation Function for the first time-series (y_30_30_10)
acf(y_30_30_10, lag.max = 20, plot = T)
##
## Call:
## lm(formula = y_30_30_10 ~ y_20_30_15)
##
## Residuals:
## Min 1Q Median 3Q Max
## -202.850 -35.762 0.349 39.126 132.700
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.392e+04 4.223e+02 32.961 <2e-16 ***
## y_20_30_15 1.665e-02 6.678e-02 0.249 0.803
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 60.77 on 178 degrees of freedom
## Multiple R-squared: 0.0003491, Adjusted R-squared: -0.005267
## F-statistic: 0.06216 on 1 and 178 DF, p-value: 0.8034plot_ly(x=~y_20_30_15, y=~y_30_30_10, type="scatter", mode="markers") %>%
add_trace(y=~(lm1$coefficients[1] + lm1$coefficients[2]*y_20_30_15), mode="lines", name="Linear Model") %>%
layout(legend = list(orientation='h'))Determining if a time series is auto-correlated is essential because it indicates how much past values in the series influence future ones, which affects forecasting and model accuracy. Testing for auto-correlation is also a key step in validating the assumptions of many statistical tests that require independence of observations. As for the correlation between time series at two different locations, assessing the degree of correlation is vital for uncovering potential relationships or interactions between them, which could provide insights into the mechanisms that underpin the observed data, reveal synchronized behaviors, or inform predictions about one series based on the other.
Problem #2
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ lubridate 1.9.3 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks plotly::filter(), stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ lubridate::origin() masks oro.nifti::origin(), RNifti::origin()
## ✖ dplyr::slice() masks plotly::slice(), oro.nifti::slice()
## ✖ tibble::view() masks RNifti::view()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(plotly)
Kidney_dataset <- read.csv("https://umich.instructure.com/files/29594707/download?download_frd=1",
header = TRUE)
library(psych)##
## Attaching package: 'psych'
##
## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha## Kidney_dataset
##
## 63 Variables 210 Observations
## --------------------------------------------------------------------------------
## case_id
## n missing distinct
## 210 0 210
##
## lowest : case_00000 case_00001 case_00002 case_00003 case_00004
## highest: case_00205 case_00206 case_00207 case_00208 case_00209
## --------------------------------------------------------------------------------
## age_at_nephrectomy
## n missing distinct Info Mean Gmd .05 .10
## 210 0 59 0.999 58.35 15.77 35.00 39.80
## .25 .50 .75 .90 .95
## 50.00 61.00 68.00 74.00 77.55
##
## lowest : 1 11 12 17 21, highest: 82 83 84 86 90
## --------------------------------------------------------------------------------
## gender
## n missing distinct
## 210 0 2
##
## Value female male
## Frequency 87 123
## Proportion 0.414 0.586
## --------------------------------------------------------------------------------
## body_mass_index
## n missing distinct Info Mean Gmd .05 .10
## 210 0 198 1 31.17 7.293 22.05 23.72
## .25 .50 .75 .90 .95
## 26.47 30.07 35.61 40.07 42.87
##
## lowest : 16.2 18.92 19.96 20.13 20.16, highest: 45.6 46.62 48.24 48.42 49.61
## --------------------------------------------------------------------------------
## comorbidities.myocardial_infarction
## n missing distinct
## 210 0 2
##
## Value false true
## Frequency 201 9
## Proportion 0.957 0.043
## --------------------------------------------------------------------------------
## comorbidities.congestive_heart_failure
## n missing distinct
## 210 0 2
##
## Value false true
## Frequency 204 6
## Proportion 0.971 0.029
## --------------------------------------------------------------------------------
## comorbidities.peripheral_vascular_disease
## n missing distinct
## 210 0 2
##
## Value false true
## Frequency 201 9
## Proportion 0.957 0.043
## --------------------------------------------------------------------------------
## comorbidities.cerebrovascular_disease
## n missing distinct
## 210 0 2
##
## Value false true
## Frequency 206 4
## Proportion 0.981 0.019
## --------------------------------------------------------------------------------
## comorbidities.dementia
## n missing distinct value
## 210 0 1 false
##
## Value false
## Frequency 210
## Proportion 1
## --------------------------------------------------------------------------------
## comorbidities.copd
## n missing distinct
## 210 0 2
##
## Value false true
## Frequency 199 11
## Proportion 0.948 0.052
## --------------------------------------------------------------------------------
## comorbidities.connective_tissue_disease
## n missing distinct
## 210 0 2
##
## Value false true
## Frequency 209 1
## Proportion 0.995 0.005
## --------------------------------------------------------------------------------
## comorbidities.peptic_ulcer_disease
## n missing distinct
## 210 0 2
##
## Value false true
## Frequency 209 1
## Proportion 0.995 0.005
## --------------------------------------------------------------------------------
## comorbidities.uncomplicated_diabetes_mellitus
## n missing distinct
## 210 0 2
##
## Value false true
## Frequency 172 38
## Proportion 0.819 0.181
## --------------------------------------------------------------------------------
## comorbidities.diabetes_mellitus_with_end_organ_damage
## n missing distinct
## 210 0 2
##
## Value false true
## Frequency 208 2
## Proportion 0.99 0.01
## --------------------------------------------------------------------------------
## comorbidities.chronic_kidney_disease
## n missing distinct
## 210 0 2
##
## Value false true
## Frequency 194 16
## Proportion 0.924 0.076
## --------------------------------------------------------------------------------
## comorbidities.hemiplegia_from_stroke
## n missing distinct
## 210 0 2
##
## Value false true
## Frequency 209 1
## Proportion 0.995 0.005
## --------------------------------------------------------------------------------
## comorbidities.leukemia
## n missing distinct
## 210 0 2
##
## Value false true
## Frequency 209 1
## Proportion 0.995 0.005
## --------------------------------------------------------------------------------
## comorbidities.malignant_lymphoma
## n missing distinct
## 210 0 2
##
## Value false true
## Frequency 208 2
## Proportion 0.99 0.01
## --------------------------------------------------------------------------------
## comorbidities.localized_solid_tumor
## n missing distinct
## 210 0 2
##
## Value false true
## Frequency 181 29
## Proportion 0.862 0.138
## --------------------------------------------------------------------------------
## comorbidities.metastatic_solid_tumor
## n missing distinct
## 210 0 2
##
## Value false true
## Frequency 198 12
## Proportion 0.943 0.057
## --------------------------------------------------------------------------------
## comorbidities.mild_liver_disease
## n missing distinct
## 210 0 2
##
## Value false true
## Frequency 206 4
## Proportion 0.981 0.019
## --------------------------------------------------------------------------------
## comorbidities.moderate_to_severe_liver_disease
## n missing distinct
## 210 0 2
##
## Value false true
## Frequency 209 1
## Proportion 0.995 0.005
## --------------------------------------------------------------------------------
## comorbidities.aids
## n missing distinct value
## 210 0 1 false
##
## Value false
## Frequency 210
## Proportion 1
## --------------------------------------------------------------------------------
## smoking_history
## n missing distinct
## 210 0 3
##
## Value current_smoker never_smoked previous_smoker
## Frequency 34 100 76
## Proportion 0.162 0.476 0.362
## --------------------------------------------------------------------------------
## age_when_quit_smoking
## n missing distinct
## 188 22 36
##
## lowest : 22 24 25 26 29
## highest: 60 64 71 73 not_applicable
## --------------------------------------------------------------------------------
## pack_years
## n missing distinct Info Mean Gmd .05 .10
## 143 67 23 0.658 9.909 16.58 0 0
## .25 .50 .75 .90 .95
## 0 0 8 34 45
##
## lowest : 0 2 5 6 7, highest: 50 64 75 115 210
## --------------------------------------------------------------------------------
## chewing_tobacco_use
## n missing distinct
## 210 0 2
##
## Value never_or_not_in_last_3mo quit_in_last_3mo
## Frequency 209 1
## Proportion 0.995 0.005
## --------------------------------------------------------------------------------
## alcohol_use
## n missing distinct
## 210 0 3
##
## Value more_than_two_daily never_or_not_in_last_3mo
## Frequency 13 95
## Proportion 0.062 0.452
##
## Value two_or_less_daily
## Frequency 102
## Proportion 0.486
## --------------------------------------------------------------------------------
## intraoperative_complications.blood_transfusion
## n missing distinct
## 210 0 2
##
## Value false true
## Frequency 207 3
## Proportion 0.986 0.014
## --------------------------------------------------------------------------------
## intraoperative_complications.injury_to_surrounding_organ
## n missing distinct
## 210 0 2
##
## Value false true
## Frequency 207 3
## Proportion 0.986 0.014
## --------------------------------------------------------------------------------
## intraoperative_complications.cardiac_event
## n missing distinct value
## 210 0 1 false
##
## Value false
## Frequency 210
## Proportion 1
## --------------------------------------------------------------------------------
## hospitalization
## n missing distinct
## 210 0 15
##
## lowest : 0 1 10 11 2
## highest: 6 7 8 9 died_before_discharge
## --------------------------------------------------------------------------------
## ischemia_time
## n missing distinct
## 200 10 39
##
## lowest : 0 10 11 12 13
## highest: 5 50 54 8 not_applicable
## --------------------------------------------------------------------------------
## radiographic_size
## n missing distinct Info Mean Gmd .05 .10
## 210 0 84 0.999 4.772 3.211 1.678 2.000
## .25 .50 .75 .90 .95
## 2.500 4.000 5.600 9.611 11.455
##
## lowest : 1.2 1.4 1.5 1.6 1.66, highest: 13.2 13.9 14.1 15 16.2
## --------------------------------------------------------------------------------
## pathologic_size
## n missing distinct Info Mean Gmd .05 .10
## 210 0 77 0.999 4.716 3.217 1.500 1.990
## .25 .50 .75 .90 .95
## 2.400 3.900 5.975 9.200 11.065
##
## lowest : 0.8 1 1.2 1.3 1.4 , highest: 13 13.6 15 16 22.5
## --------------------------------------------------------------------------------
## malignant
## n missing distinct
## 210 0 2
##
## Value false true
## Frequency 18 192
## Proportion 0.086 0.914
## --------------------------------------------------------------------------------
## pathology_t_stage
## n missing distinct
## 210 0 7
##
## Value 0 1a 1b 2a 2b 3 4
## Frequency 18 92 41 9 3 44 3
## Proportion 0.086 0.438 0.195 0.043 0.014 0.210 0.014
## --------------------------------------------------------------------------------
## pathology_n_stage
## n missing distinct
## 210 0 3
##
## Value 0 1 X
## Frequency 90 6 114
## Proportion 0.429 0.029 0.543
## --------------------------------------------------------------------------------
## pathology_m_stage
## n missing distinct
## 210 0 3
##
## Value 0 1 X
## Frequency 93 17 100
## Proportion 0.443 0.081 0.476
## --------------------------------------------------------------------------------
## tumor_histologic_subtype
## n missing distinct
## 210 0 12
##
## lowest : angiomyolipoma chromophobe clear_cell_papillary_rcc clear_cell_rcc mest
## highest: papillary rcc_unclassified spindle_cell_neoplasm urothelial wilms
## --------------------------------------------------------------------------------
## tumor_necrosis
## n missing distinct
## 187 23 2
##
## Value false true
## Frequency 144 43
## Proportion 0.77 0.23
## --------------------------------------------------------------------------------
## tumor_isup_grade
## n missing distinct Info Mean Gmd
## 172 38 4 0.826 2.238 0.8374
##
## Value 1 2 3 4
## Frequency 26 93 39 14
## Proportion 0.151 0.541 0.227 0.081
##
## For the frequency table, variable is rounded to the nearest 0
## --------------------------------------------------------------------------------
## clavien_surgical_complications
## n missing distinct
## 209 1 7
##
## Value 0 1 2 3a 3b 4 5
## Frequency 149 27 15 6 7 4 1
## Proportion 0.713 0.129 0.072 0.029 0.033 0.019 0.005
## --------------------------------------------------------------------------------
## er_visit
## n missing distinct
## 208 2 2
##
## Value false true
## Frequency 191 17
## Proportion 0.918 0.082
## --------------------------------------------------------------------------------
## readmission
## n missing distinct
## 208 2 2
##
## Value false true
## Frequency 187 21
## Proportion 0.899 0.101
## --------------------------------------------------------------------------------
## estimated_blood_loss
## n missing distinct Info Mean Gmd .05 .10
## 209 1 43 0.994 314 339.9 25 50
## .25 .50 .75 .90 .95
## 75 200 400 800 1000
##
## lowest : 0 5 10 20 25, highest: 1200 1300 1500 1800 2000
## --------------------------------------------------------------------------------
## surgery_type
## n missing distinct
## 210 0 3
##
## Value laparoscopic open robotic
## Frequency 28 60 122
## Proportion 0.133 0.286 0.581
## --------------------------------------------------------------------------------
## surgical_procedure
## n missing distinct
## 210 0 2
##
## Value partial_nephrectomy radical_nephrectomy
## Frequency 140 70
## Proportion 0.667 0.333
## --------------------------------------------------------------------------------
## surgical_approach
## n missing distinct
## 210 0 2
##
## Value Retroperitoneal Transperitoneal
## Frequency 39 171
## Proportion 0.186 0.814
## --------------------------------------------------------------------------------
## operative_time
## n missing distinct Info Mean Gmd .05 .10
## 208 2 151 1 243 106.6 120.4 139.4
## .25 .50 .75 .90 .95
## 174.8 229.0 283.2 372.3 437.8
##
## lowest : 44 47 64 91 97, highest: 530 553 586 592 613
## --------------------------------------------------------------------------------
## cytoreductive
## n missing distinct
## 210 0 2
##
## Value false true
## Frequency 193 17
## Proportion 0.919 0.081
## --------------------------------------------------------------------------------
## positive_resection_margins
## n missing distinct
## 210 0 2
##
## Value false true
## Frequency 199 11
## Proportion 0.948 0.052
## --------------------------------------------------------------------------------
## last_preop_egfr.value
## n missing distinct
## 153 57 49
##
## lowest : >=90 17 38 41 42 , highest: 86 87 88 89 90
## --------------------------------------------------------------------------------
## last_preop_egfr.days_before_nephrectomy
## n missing distinct Info Mean Gmd .05 .10
## 153 57 54 0.998 23.39 23.93 1.0 2.0
## .25 .50 .75 .90 .95
## 6.0 15.0 36.0 58.2 69.0
##
## lowest : 1 2 3 4 5, highest: 75 76 77 83 90
## --------------------------------------------------------------------------------
## first_postop_egfr.value
## n missing distinct
## 157 53 59
##
## lowest : >=90 25 27 28 29 , highest: 88 89 9 90 age<16
## --------------------------------------------------------------------------------
## first_postop_egfr.days_before_nephrectomy
## n missing distinct Info Mean Gmd .05 .10
## 157 53 112 1 437.6 526.5 78.0 82.2
## .25 .50 .75 .90 .95
## 103.0 154.0 500.0 1162.2 2078.4
##
## lowest : 76 77 78 79 80, highest: 2317 2380 2393 2473 2503
## --------------------------------------------------------------------------------
## last_postop_egfr.value
## n missing distinct
## 157 53 60
##
## lowest : >=90 27 28 30 32 , highest: 88 89 9 90 age<16
## --------------------------------------------------------------------------------
## last_postop_egfr.days_before_nephrectomy
## n missing distinct Info Mean Gmd .05 .10
## 157 53 147 1 1048 818.1 149.8 197.6
## .25 .50 .75 .90 .95
## 442.0 932.0 1487.0 2256.8 2426.0
##
## lowest : 76 77 93 100 114, highest: 2499 2503 2584 2604 3070
## --------------------------------------------------------------------------------
## vital_status
## n missing distinct
## 210 0 2
##
## Value censored dead
## Frequency 189 21
## Proportion 0.9 0.1
## --------------------------------------------------------------------------------
## vital_days_after_surgery
## n missing distinct Info Mean Gmd .05 .10
## 210 0 185 1 956.8 875.2 15.45 23.80
## .25 .50 .75 .90 .95
## 271.25 790.50 1427.50 2288.10 2431.60
##
## lowest : 0 2 8 9 11, highest: 2539 2607 2611 2900 3070
## --------------------------------------------------------------------------------
## voxel_spacing.x_spacing
## n missing distinct Info Mean Gmd .05 .10
## 210 0 125 0.999 0.7966 0.1216 0.6565 0.6832
## .25 .50 .75 .90 .95
## 0.7183 0.7812 0.8682 0.9631 0.9766
##
## lowest : 0.4375 0.451172 0.518 0.546875 0.550781
## highest: 0.970703 0.976562 0.976562 0.984 1.041
## --------------------------------------------------------------------------------
## voxel_spacing.y_spacing
## n missing distinct Info Mean Gmd .05 .10
## 210 0 125 0.999 0.7966 0.1216 0.6565 0.6832
## .25 .50 .75 .90 .95
## 0.7183 0.7812 0.8682 0.9631 0.9766
##
## lowest : 0.4375 0.451172 0.518 0.546875 0.550781
## highest: 0.970703 0.976562 0.976562 0.984 1.041
## --------------------------------------------------------------------------------
## voxel_spacing.z_spacing
## n missing distinct Info Mean Gmd .05 .10
## 210 0 10 0.92 3.162 1.93 0.50 0.50
## .25 .50 .75 .90 .95
## 1.25 3.00 5.00 5.00 5.00
##
## Value 0.50 1.00 1.25 1.50 2.00 2.50 3.00 3.75 4.00 5.00
## Frequency 31 18 5 2 4 23 38 1 2 86
## Proportion 0.148 0.086 0.024 0.010 0.019 0.110 0.181 0.005 0.010 0.410
##
## For the frequency table, variable is rounded to the nearest 0
## --------------------------------------------------------------------------------## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2Kidney_DataDescr <- ExpData(Kidney_dataset, type=2, fun=c("mean", "median", "var")) %>% datatable()
library(dplyr)
Kidney_dataset <- select(Kidney_dataset, -starts_with("comorbidities."))
Kidney_DataDescr <- ExpData(Kidney_dataset, type=2, fun=c("mean", "median", "var")) %>% datatable()
plot_ly(Kidney_dataset, x = ~surgery_type, y = ~radiographic_size, split = ~surgery_type,
type = 'violin', box = list(visible = T), meanline = list(visible = T)) %>%
layout(title="Distribution of radiographic_size by surgery_type",
xaxis=list(title="Surgery Type"), yaxis=list(title="Tumor Size"),
legend = list(title="Surgery Type", orientation = 'h'))plot_ly(Kidney_dataset, type = 'violin') %>%
add_trace(x = ~surgical_approach[Kidney_dataset$gender == 'female'],
y = ~vital_days_after_surgery[Kidney_dataset$gender == 'female'],
legendgroup = 'Yes', scalegroup = 'Yes', name = 'Female',
side='negative', box=list(visible=T), meanline=list(visible=T),
color = I("blue"), points = 'all', jitter = 0.5, pointpos = -1) %>%
add_trace(x = ~surgical_approach[Kidney_dataset$gender == 'male'],
y = ~vital_days_after_surgery[Kidney_dataset$gender == 'male'],
legendgroup = 'Yes', scalegroup = 'Yes', name = 'Male',
side='positive', box=list(visible=T), meanline=list(visible=T),
color = I("yellow"), points = 'all', jitter = 0.5, pointpos = +1) %>%
layout(title="Distribution of vital_days_after_surgery by surgical_approach and `ender",
xaxis=list(title="Surgical Approach"), yaxis=list(title="Vital days after surgery"),
legend = list(title="Surgical Approach", orientation = 'h')) #, ## Warning: Can't display both discrete & non-discrete data on same axis##
## Attaching package: 'extraDistr'
##
## The following object is masked from 'package:purrr':
##
## rdunif##
## Attaching package: 'pracma'
##
## The following objects are masked from 'package:psych':
##
## logit, polar
##
## The following object is masked from 'package:purrr':
##
## cross
##
## The following object is masked from 'package:oro.nifti':
##
## magic# marginal domains, # summary(Kidney_dataset$vital_days_after_surgery)
len <- 300
t <- 1:len
x_dom <- seq(0, 3070, length=len) # range of vital_days_after_surgery
y_dom <- seq(16, 50, length=len) # range of body_mass_index
rho <- 0.3 # initial dependence (correlation), -1 < rho = Corr(t,phi) < 1
# Base function definitions
# f(t), Kidney_dataset$vital_days_after_surgery
pdf_x <- density(Kidney_dataset$vital_days_after_surgery, bw = "sj")
f <- approxfun(pdf_x, rule=1:2) # density as an abstract function
fMax <- max(abs(pdf_x$y))
F1 <- function (x) { # CDF as an abstract function
# message(paste('input to F1 =',x))
if (x < 1) { val = 0 }
else if (x > len) { val = 1 }
else { val=integral(f, xmin=x_dom[1]-200, xmax=x_dom[x],
reltol=1e-12, method="Kronrod")
}
# val=integral(f, xmin=x_dom[1]-200, xmax=x, reltol=1e-12, method="Kronrod")
return(val)
} # F1(100)
# g(phi), Kidney_dataset$body_mass_index
pdf_y <- density(Kidney_dataset$body_mass_index)
g <- approxfun(pdf_y, rule=1:2) # density as an abstract function
gMax <- max(abs(pdf_y$y))
G1 <- function (y) {
if (y < 1) { val = 0 }
else if (y > len) { val = 1 }
else { val = integral(g, y_dom[1], y_dom[y], method = "Kronrod") }
return(val)
} # G1(100)
# Gauss-copula
# C_{XY}^{Gauss} (x,y)=\frac{1}{\sqrt{1-ρ^2}} \times \exp \left (-\frac{ρ^2 (a^2+b^2 )-2abρ)}{2(1-ρ^2 )} \right ),
# where the pairwise correlation $-1 < \rho < 1$ and
# a(x)= \sqrt{2}\times erf^{-1} (2\times F_X (x)-1),
# b(y)= \sqrt{2}\times erf^{-1} (2\times F_Y (y)-1),
# and the error function is
# $$erf^{-1}(s)= \frac{2}{\sqrt{\pi}}\int_{0}^{\infty}{e^{-t^2}dt}.$$
CGauss <- function (t=len/2, s=len/2) {
if (t < 1) { A= -1e+9 }
else if (t >= len) { A= 1e+9 }
else { A = sqrt(2) * erfinv(2*(F1(t))-1) }
if (s < 1) { B= -1e+9 }
else if (s > len) { B= 1e+9 }
else { B = sqrt(2) * erfinv(2*(G1(s))-1) }
# message(paste('A=', A, 'B=', B))
val = (1/(sqrt(1-rho^2))) *
exp(-(rho^2*(A^2+B^2) - 2*A*B*rho)/(2*(1-rho^2)))
if (is.finite(val)) return (val)
else return (0)
} # CGauss(1, 17)
# F(t)
fJoint <- function (t=len/2, s=len/2, rho = rho) {
return ( f(x_dom[t]) * g(y_dom[s]) * CGauss(t, s) )
} # fJoint(1,18)
### Part 1: Classical Time Dynamics: $\phi=0$.
# plot(s, f(s)/fMax)
# plot(s, g(s)/gMax)
# # plot(s, 2*g(rep(0, length(s)))/gMax)
# # plot(s, log(F(t=s, phi=0, rho=0.1)))
# plot(s, log(F(t=s, phi=rep(0, length(s)), rho=0.1)))
# plot(s, CGauss(s))
library(plotly)
y1 <- 1000*f(x_dom) # f: first marginal
y2 <- 10* g(y_dom) # g: second marginal
y3 <- c() # Gaussian Copula (. , s=sMax/2)
for (i in 1:len) y3[i]=CGauss(i, 55) # or y_dom[len/2]
y4 <- c() # joint density f(t,t)
for (i in 1:len) y4[i] <- 10000*fJoint(t=i, s=55, rho)
plot_ly() %>%
add_trace(x=~t, y=~y1, type='scatter', mode='markers+lines',
name="f: first marginal") %>%
add_trace(x=~t, y=~y2, type='scatter', mode='markers+lines',
name="g: second marginal", opacity=1) %>%
add_trace(x=~t, y=~y3, type='scatter', mode='markers+lines',
name="C_Gauss(t,s): Gaussian Copula", opacity=0.5) %>%
add_trace(x=~t, y=~y4, type='scatter',
mode='markers+lines',
name="f(t,s, rho)): Joint Density") %>%
layout(title = "Copula-based Modeling of Joint Distribution",
xaxis = list(title = "time"),
yaxis = list(title = "Functional Value"), # showlegend = FALSE,
legend = list(orientation = 'h'),
scene = list(aspectmode = "manual",
aspectratio = list(x=1, y=1)))#### Path Kynamics over curves \gamma:\mathbb{R}\to \mathbb{C}
### Part 2: Path Kynamics over curves \gamma:\mathbb{R}\to \mathbb{C}
### Kime tube gamma(k1,k2) on torus ########## CARTESIAN
#### Parametric curve on torus
# define kime-surface
#
#xsurf <- x_dom[t] # longitudinal trajectory
#ysurf <- y_dom[t] # phase trajectory
#zsurf <- mat<-matrix(0, nrow=len, ncol=len) # joint density f(t,s)
#for (i in 1:len) {
# for (j in 1:len) {
# zsurf[i,j] <- 100000*fJoint(t=i, s=j, rho)
# }
# if (i %% 30 ==0) message(paste0('Progress: ... ', round((100*i)/len, 0), '%'))
#}
# image(zsurf)
#plot_ly(z= zsurf, type="heatmap") # %>% hide_colorbar()
# 3D Kime-surface Plot
#f <- list(family = "Courier New, monospace", size = 18, color = "black")
#xlab <- list(title = "post surgery vitality", titlefont = f)
#ylab <- list(title = "BMI", titlefont = f)
#zlab <- list(title = "Density surface", titlefont = f)
#
#plot_ly(x=~xsurf, y=~ysurf, z=~t(zsurf), type="surface", opacity=0.95,
# showlegend=FALSE, name="Joint Density-Surface") %>%
# # add vital days after surgery projection curve, for a fixed BMI=22
# add_trace(x=~xsurf, y=ysurf[len/2], z=~zsurf[, len/2], type="scatter3d",
# mode="lines", line=list(width=20), name="BMI=22") %>%
# # add BMI projection curve vital days after surgery=1000
# add_trace(x=xsurf[100], y=~ysurf, z=~zsurf[100,], type="scatter3d",
# mode="lines", line=list(width=20), name="1000 days post surgery") %>%
# hide_colorbar() %>%
# layout(title = "Joint Bivariate Density (Gaussian copula, Rho=0.3)",
# scene = list(aspectmode = "manual",
# xaxis = xlab, yaxis = ylab, zaxis = zlab,
# aspectratio = list(x=1, y=1, z=1))) # 1:1:1 aspect ratio
# Model 1 -- may need to split 80:20 the data into training set (for model fitting) and
# testing set (for independent valuation)
Kidney_dataset$vital_status <- as.factor(Kidney_dataset$vital_status)
table(Kidney_dataset$vital_status)##
## censored dead
## 189 21# 80% training + 20% testing
set.seed(1234)
subset_int <- sample(nrow(Kidney_dataset), floor(nrow(Kidney_dataset)*0.8))
Kidney_dataset_train <- Kidney_dataset[subset_int, ]
Kidney_dataset_test <- Kidney_dataset[-subset_int, ]
# Build logistic regression model (based on training set)
Log.Reg.Model <- glm(vital_status ~ age_at_nephrectomy+tumor_isup_grade+estimated_blood_loss+
operative_time+pathologic_size+body_mass_index,
data = Kidney_dataset_train, family = "binomial")
summary(Log.Reg.Model)##
## Call:
## glm(formula = vital_status ~ age_at_nephrectomy + tumor_isup_grade +
## estimated_blood_loss + operative_time + pathologic_size +
## body_mass_index, family = "binomial", data = Kidney_dataset_train)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.925453 2.573797 -2.302 0.02132 *
## age_at_nephrectomy 0.015919 0.023118 0.689 0.49108
## tumor_isup_grade 1.154501 0.396964 2.908 0.00363 **
## estimated_blood_loss -0.001089 0.001182 -0.921 0.35680
## operative_time 0.001753 0.003293 0.532 0.59453
## pathologic_size 0.081326 0.081302 1.000 0.31717
## body_mass_index -0.019959 0.052095 -0.383 0.70163
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 98.771 on 136 degrees of freedom
## Residual deviance: 77.690 on 130 degrees of freedom
## (31 observations deleted due to missingness)
## AIC: 91.69
##
## Number of Fisher Scoring iterations: 6## Classes and Methods for R originally developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University (2002-2015),
## by and under the direction of Simon Jackman.
## hurdle and zeroinfl functions by Achim Zeileis.## $N
## [1] 137
##
## $R2
## [1] 0.2776233## fitting null model for pseudo-r2## llh llhNull G2 McFadden r2ML r2CU
## -38.8450735 -49.3853113 21.0804755 0.2134286 0.1426183 0.2776233# Prediction (using testing data)
Log.Reg.Predict <- predict(Log.Reg.Model, Kidney_dataset_test, type = "response")
# Check Accuracy
library(caret)## Loading required package: lattice
## Registered S3 methods overwritten by 'pROC':
## method from
## print.roc fmsb
## plot.roc fmsb
##
## Attaching package: 'caret'
##
## The following object is masked from 'package:purrr':
##
## lift# Confusion Matrix
confusionMatrix(data = factor(as.factor(Log.Reg.Predict <= 0.02), labels=c("censored", "dead")),
reference = Kidney_dataset_test$vital_status)## Confusion Matrix and Statistics
##
## Reference
## Prediction censored dead
## censored 27 2
## dead 2 1
##
## Accuracy : 0.875
## 95% CI : (0.7101, 0.9649)
## No Information Rate : 0.9062
## P-Value [Acc > NIR] : 0.8239
##
## Kappa : 0.2644
##
## Mcnemar's Test P-Value : 1.0000
##
## Sensitivity : 0.9310
## Specificity : 0.3333
## Pos Pred Value : 0.9310
## Neg Pred Value : 0.3333
## Prevalence : 0.9062
## Detection Rate : 0.8438
## Detection Prevalence : 0.9062
## Balanced Accuracy : 0.6322
##
## 'Positive' Class : censored
## # print(paste('Accuracy is ', 1 - MisClassificError))
# Area Under Curve
library(ROCR)
# replace NA predictions to 0 (trivial) probabilities
Log.Reg.Predict <- Log.Reg.Predict %>% replace(is.na(.), 0)
TestDataDependentVar <- Kidney_dataset_test$vital_status
pred <- ROCR::prediction(Log.Reg.Predict, TestDataDependentVar)
perf <- ROCR::performance(prediction.obj = pred, measure = "tpr", x.measure = "fpr")
plot(perf, col = "black", lty = 2, lwd = 2)
abline(a = 0, b = 1, col = "red")
# Calculating AUC
auc <- ROCR::performance(prediction.obj = pred, measure = "auc")
auc <- as.numeric(auc@y.values)
minauc <- (round(auc, digits = 2))
minauct <- paste(c("min (AUC) = "), minauc, sep = "")
legend(0.52, 0.5, c(minauct, "\n"), cex = 1, box.col = "Red")
##
## Attaching package: 'e1071'
##
## The following object is masked from 'package:pracma':
##
## sigmoid#nb_model <- naiveBayes(vital_status ~ age_at_nephrectomy+tumor_isup_grade+estimated_blood_loss+
# operative_time+pathologic_size+body_mass_index,
# data = Kidney_dataset_train, family = "binomial")
nb_model <- naiveBayes(surgical_approach ~ vital_status, data = Kidney_dataset_train)
print(nb_model)##
## Naive Bayes Classifier for Discrete Predictors
##
## Call:
## naiveBayes.default(x = X, y = Y, laplace = laplace)
##
## A-priori probabilities:
## Y
## Retroperitoneal Transperitoneal
## 0.2142857 0.7857143
##
## Conditional probabilities:
## vital_status
## Y censored dead
## Retroperitoneal 0.83333333 0.16666667
## Transperitoneal 0.90909091 0.09090909train_predictions <- predict(nb_model, Kidney_dataset_train)
#table(train_predictions, Kidney_dataset_train$Outcome)
accuracy <- sum(train_predictions == Kidney_dataset_train$Outcome) / nrow(Kidney_dataset_train)
print(accuracy)## [1] 0statistical significance of each model parameter (beta effects) corresponding to each explanatory predictor variable:
age_at_nephrectomy 0.015919
tumor_isup_grade 1.154501
estimated_blood_loss -0.001089
operative_time 0.001753
pathologic_size 0.081326
body_mass_index -0.019959