HW3

#Question 1 Why is it important to find whether a time-series is auto-correlated, and whether the time series at two voxels are correlated?

Time-series autocorrelation involves analyzing how consecutive measurements of a single variable evolve over time. For fMRI data, this means observing changes in signal intensity across time points, which helps to determine if a particular data point or brain region is being engaged in the process being studied.

Examining the correlation between time series data from two distinct voxels can help distinguish whether an observed correlation at a single voxel represents background noise that is common across all tissues or a genuine effect. This distinction is important for accurately interpreting neural activity and its relevance to physiological processes.

# Download the 4D fMRI hyper-volume (3D space * 1D time)
# library(brainR)  # for NIFTI file input
library(oro.nifti)

## oro.nifti 0.11.4

library(fmri)      # for NIFTI file input
fMRIURL <- "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

# Extract the 1D time course of a specific 3D voxel (say the one at x=30, y=30, z=10):
 plot(fMRIVolume[30, 30, 10,], type='l', main="Time Series of 3D Voxel \n (x=30, y=30, z=10)", col="blue")
# Smooth the noisy fMRI time-series
 x1 <- c(1:180)
 y1 <- loess(fMRIVolume[30, 30, 10,]~ x1, family = "gaussian")
 lines(x1, smooth(fMRIVolume[30, 30, 10,]), col = "red", lwd = 2);
 lines(ksmooth(x1, fMRIVolume[30, 30, 10,], kernel = "normal", bandwidth = 5), col = "green", lwd = 3);
 legend("bottomright", legend=c("(raw) fMRI", "smooth(fMRI)", "ksmooth(fMRI"),
        col=c("blue", "red", "green"), lty=1, cex=0.8, y.intersp=0.8)

library(plotly)

## 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':
## 
##     layout

x1 <- 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")

# Compute their correlation
cor(y_30_30_10, y_20_30_15)

## [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)

# fit as linear model
lm1 <- lm(y_30_30_10 ~ y_20_30_15); summary(lm1)

## 
## 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.8034

plot_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'))

#Question 2 For predicting patient vital status, there were five potential variables that were considered, which included the age at surgery, duration of the surgery, blood loss, tumor grade, and body mass index (BMI). Within the analysis, only the tumor ISUP grade showed statistical significance, with a p-value of 0.00363 which is less than 0.05. The model demonstrated an accuracy of 87.5% in predicting the patient’s survival based on these variables when applied to the test data.

library("tidyverse")

## ── 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()
## ✖ dplyr::slice()  masks plotly::slice(), oro.nifti::slice()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(plotly)
library(nnet) # For multinomial logistic regression
library(lmtest)

## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

Kidney_dataset <- read.csv("https://umich.instructure.com/files/29594707/download?download_frd=1", 
                    header = TRUE)

Kidney_dataset$vital_status <- as.factor(Kidney_dataset$vital_status)
table(Kidney_dataset$vital_status)

## 
## censored     dead 
##      189       21

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,  ] 
Kidney_dataset$vital_status <- ifelse(Kidney_dataset$vital_status == "censored", 0, 1)

# 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

library(fmsb)
library(pscl)

## 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.

NagelkerkeR2(Log.Reg.Model)

## $N
## [1] 137
## 
## $R2
## [1] 0.2776233

# Other pseudo r-squared
pR2(Log.Reg.Model)

## 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        
## 

Using the K-nearest neighbor method:

library(class)
library(caret)
library(gmodels)


Kidney_dataset$malignant <- as.numeric(as.factor(Kidney_dataset$malignant))

# Split the dataset first to avoid data leakage
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,]

# Identify numeric predictors for imputation
numericPredictors <- sapply(Kidney_dataset_train, is.numeric)
numericPredictors["vital_status"] <- FALSE  # Exclude outcome variable

# Impute missing values for numeric predictors in both training and test datasets
for(col in names(numericPredictors[numericPredictors])) {
  Kidney_dataset_train[[col]] <- ifelse(is.na(Kidney_dataset_train[[col]]), mean(Kidney_dataset_train[[col]], na.rm = TRUE), Kidney_dataset_train[[col]])
  Kidney_dataset_test[[col]] <- ifelse(is.na(Kidney_dataset_test[[col]]), mean(Kidney_dataset_test[[col]], na.rm = TRUE), Kidney_dataset_test[[col]])
}

# Scale numeric predictors
Kidney_dataset_train[numericPredictors] <- scale(Kidney_dataset_train[numericPredictors])
Kidney_dataset_test[numericPredictors] <- scale(Kidney_dataset_test[numericPredictors])

# Ensure no NA, NaN, or Inf values remain
Kidney_dataset_train <- na.omit(Kidney_dataset_train)
Kidney_dataset_test <- na.omit(Kidney_dataset_test)

# Prepare for KNN
X_train <- Kidney_dataset_train[, numericPredictors]
Y_train <- Kidney_dataset_train$vital_status
X_test <- Kidney_dataset_test[, numericPredictors]
Y_test <- Kidney_dataset_test$vital_status

# Train KNN
k <- 5
knn_fit <- knn(train = X_train, test = X_test, cl = Y_train, k = k)
Y_test <- factor(Y_test, levels = c(0, 1))
knn_fit <- factor(knn_fit, levels = c(0, 1))

# Calculate accuracy
accuracy <- mean(Y_test == knn_fit)
print(paste("Accuracy:", accuracy))

## [1] "Accuracy: 0.904761904761905"

# Confusion Matrix
confusionMatrix <- table(Predicted = knn_fit, Actual = Y_test)
print(confusionMatrix)

##          Actual
## Predicted  0  1
##         0 38  3
##         1  1  0

# Using caret's confusionMatrix function for detailed metrics
confMatrix <- confusionMatrix(confusionMatrix)
print(confMatrix)

## Confusion Matrix and Statistics
## 
##          Actual
## Predicted  0  1
##         0 38  3
##         1  1  0
##                                           
##                Accuracy : 0.9048          
##                  95% CI : (0.7738, 0.9734)
##     No Information Rate : 0.9286          
##     P-Value [Acc > NIR] : 0.8217          
##                                           
##                   Kappa : -0.037          
##                                           
##  Mcnemar's Test P-Value : 0.6171          
##                                           
##             Sensitivity : 0.9744          
##             Specificity : 0.0000          
##          Pos Pred Value : 0.9268          
##          Neg Pred Value : 0.0000          
##              Prevalence : 0.9286          
##          Detection Rate : 0.9048          
##    Detection Prevalence : 0.9762          
##       Balanced Accuracy : 0.4872          
##                                           
##        'Positive' Class : 0               
## 

# Ensure both Y_test and predictions are factors with explicit levels
Y_test <- factor(Y_test, levels = c(0, 1))
knn_fit <- factor(knn_fit, levels = c(0, 1))

# Calculate accuracy for the testing set
accuracy <- mean(Y_test == knn_fit)
cat("Accuracy on the Testing Set:", accuracy, "\n\n")

## Accuracy on the Testing Set: 0.9047619

# Confusion Matrix for the testing set
confusionMatrix <- table(Predicted = knn_fit, Actual = Y_test)
print("Confusion Matrix:")

## [1] "Confusion Matrix:"

print(confusionMatrix)

##          Actual
## Predicted  0  1
##         0 38  3
##         1  1  0

# Detailed metrics using caret's confusionMatrix
# Recompute confusion matrix as a 'table' might be required by confusionMatrix function from caret
confusionMatrix <- as.table(confusionMatrix)
confMatrixDetailed <- confusionMatrix(data = knn_fit, reference = Y_test)
print("Detailed Evaluation Metrics:")

## [1] "Detailed Evaluation Metrics:"

print(confMatrixDetailed)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  0  1
##          0 38  3
##          1  1  0
##                                           
##                Accuracy : 0.9048          
##                  95% CI : (0.7738, 0.9734)
##     No Information Rate : 0.9286          
##     P-Value [Acc > NIR] : 0.8217          
##                                           
##                   Kappa : -0.037          
##                                           
##  Mcnemar's Test P-Value : 0.6171          
##                                           
##             Sensitivity : 0.9744          
##             Specificity : 0.0000          
##          Pos Pred Value : 0.9268          
##          Neg Pred Value : 0.0000          
##              Prevalence : 0.9286          
##          Detection Rate : 0.9048          
##    Detection Prevalence : 0.9762          
##       Balanced Accuracy : 0.4872          
##                                           
##        'Positive' Class : 0               
##