Introduction

Motivation

Parkinson’s disease is a progressive neurodegenerative disorder that affects millions of people worldwide, leading to a gradual decline in motor and cognitive function. Early detection is crucial for managing symptoms and slowing disease progression. In this project, we aim to analyze a relevant dataset to train a machine learning model capable of accurately predicting early signs of Parkinson’s disease. Our goal is to enable timely intervention, allowing patients to take preventive measures and improve their quality of life.

The Core Research Question

While this project addresses several related questions, our primary objective is to identify the most effective model for detecting Parkinson’s disease. Using a variety of statistical and machine learning techniques covered in our data mining course, we aim to compare different models and select the one that best generalizes the data with minimal misclassification.

Data Description

The Parkinson’s dataset used in this project is sourced from the UC Irvine Machine Learning Repository. It contains 195 voice recordings (the UCI page lists 197 instances, but the file provided has 195) and 22 features, representing voice measurements collected from 31 individuals—23 of whom have been diagnosed with Parkinson’s Disease (the recording names in the file actually contain 32 distinct subject IDs, 8 healthy and 24 with Parkinson’s; the analysis below uses the IDs found in the file). The features are designed to quantify changes in vocal attributes, which are often affected by the motor control impairments associated with Parkinson’s. These measurements capture variations in frequency, amplitude, and other vocal signal characteristics, providing a basis for distinguishing between healthy and affected individuals.

Methodology

Exploring the data.

parkinsons_data <- read.csv("parkinsons.csv")


glue("The dimensions of the given dataset is {dim(parkinsons_data)[1]} X {dim(parkinsons_data)[2]}")

## The dimensions of the given dataset is 195 X 24

glue("The number of empty values in the dataset are {sum(is.na(parkinsons_data))}")

## The number of empty values in the dataset are 0

# Basic stats of all the feature variables

#summary(parkinsons_data)

cor_matrix <- cor(parkinsons_data[ , sapply(parkinsons_data, is.numeric)])


corrplot(cor_matrix, method = 'color', number.cex = 0.7, type = 'lower', tl.cex = 0.7)

As it is evident from the heatmap above, the data has multi-collinearity. The blue color indicates a positive correlation and the red explains a negeative correlation anything in between indicates a low correlation.

Exploring the questions in the proposal by generating various models and generating relevant metrics for analysis

Question 1.

logistic_model <- glm(status ~ ., data = parkinsons_data[, -1], family = "binomial")

summary(logistic_model)

## 
## Call:
## glm(formula = status ~ ., family = "binomial", data = parkinsons_data[, 
##     -1])
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)  
## (Intercept)      -1.224e+01  1.730e+01  -0.707   0.4794  
## MDVP.Fo.Hz.      -1.187e-02  2.005e-02  -0.592   0.5538  
## MDVP.Fhi.Hz.     -2.448e-03  3.908e-03  -0.626   0.5312  
## MDVP.Flo.Hz.      3.359e-03  1.116e-02   0.301   0.7635  
## MDVP.Jitter...   -1.413e+03  1.131e+03  -1.249   0.2117  
## MDVP.Jitter.Abs. -4.242e+04  8.770e+04  -0.484   0.6286  
## MDVP.RAP          1.813e+03  1.240e+05   0.015   0.9883  
## MDVP.PPQ         -1.916e+03  1.830e+03  -1.047   0.2951  
## Jitter.DDP        6.729e+02  4.136e+04   0.016   0.9870  
## MDVP.Shimmer      3.480e+02  9.340e+02   0.373   0.7095  
## MDVP.Shimmer.dB.  2.354e+01  3.072e+01   0.766   0.4435  
## Shimmer.APQ3      1.286e+04  1.083e+05   0.119   0.9054  
## Shimmer.APQ5     -2.627e+02  4.159e+02  -0.632   0.5277  
## MDVP.APQ          2.181e+02  3.667e+02   0.595   0.5519  
## Shimmer.DDA      -4.563e+03  3.612e+04  -0.126   0.8995  
## NHR               2.059e+01  5.085e+01   0.405   0.6856  
## HNR               5.209e-02  2.052e-01   0.254   0.7997  
## RPDE             -3.213e+00  4.760e+00  -0.675   0.4996  
## DFA               7.563e+00  8.721e+00   0.867   0.3858  
## spread1           4.647e-02  1.674e+00   0.028   0.9778  
## spread2           1.047e+01  6.051e+00   1.731   0.0835 .
## D2                1.224e+00  1.410e+00   0.868   0.3856  
## PPE               3.658e+01  2.458e+01   1.488   0.1367  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 217.647  on 194  degrees of freedom
## Residual deviance:  91.027  on 172  degrees of freedom
## AIC: 137.03
## 
## Number of Fisher Scoring iterations: 9

vif(logistic_model)

##      MDVP.Fo.Hz.     MDVP.Fhi.Hz.     MDVP.Flo.Hz.   MDVP.Jitter... 
##     8.667982e+00     1.563976e+00     2.984282e+00     7.108522e+01 
## MDVP.Jitter.Abs.         MDVP.RAP         MDVP.PPQ       Jitter.DDP 
##     3.337170e+01     2.770177e+05     7.542811e+01     2.772481e+05 
##     MDVP.Shimmer MDVP.Shimmer.dB.     Shimmer.APQ3     Shimmer.APQ5 
##     1.011660e+03     9.703918e+01     4.778157e+06     8.884253e+01 
##         MDVP.APQ      Shimmer.DDA              NHR              HNR 
##     7.167481e+01     4.785237e+06     3.816464e+00     6.421831e+00 
##             RPDE              DFA          spread1          spread2 
##     3.039954e+00     2.836295e+00     2.157266e+01     2.057735e+00 
##               D2              PPE 
##     2.099320e+00     3.087806e+01

As you can see from the above summary results none of the predictors are significant this is because of high multi-collinearity between the predictors.

As you can see from the output of the vif function from cars library out 22 features only 6 features have the coefficient as less than 5. This means all the other features bring in multicollinearity.

But from the Residual deviance we can tell that the model fit is very good as their is a significant drop when compared to the null deviance.

So, this model can still be considered.

Now lets evaluate the performance of this model using test and train set.

Data Split

set.seed(32)

# The 195 recordings come from 32 subject IDs (6-7 recordings each); the person is encoded in the
# `name` column (e.g. phon_R01_S01_1 -> subject S01). A random split by recording would put
# recordings of the same person in both train and test, so whole subjects are held out instead,
# stratified by diagnosis (2 healthy and 7 Parkinson's subjects go to the test set).
subject <- sub("_[0-9]+$", "", parkinsons_data$name)
subject_status <- tapply(parkinsons_data$status, subject, max)

test_subjects <- c(sample(names(subject_status)[subject_status == 0], 2),
                   sample(names(subject_status)[subject_status == 1], 7))

train <- which(!subject %in% test_subjects)

glue("Training set: {length(train)} recordings from {length(unique(subject[train]))} people. Test set: {nrow(parkinsons_data) - length(train)} recordings from {length(test_subjects)} people ({sum(parkinsons_data$status[-train] == 0)} healthy, {sum(parkinsons_data$status[-train] == 1)} Parkinson's).")

## Training set: 140 recordings from 23 people. Test set: 55 recordings from 9 people (12 healthy, 43 Parkinson's).

logistic_model <- glm(status ~ ., data = parkinsons_data[train, -1], family = "binomial")

prob <- predict(logistic_model, newdata = parkinsons_data[-train, -1], type = "response")

pred <- ifelse(prob > 0.5, 1, 0)
#pred

conf_mat <- conf_table(parkinsons_data[-train, ]$status, pred)

accuracy_log <- round(mean(parkinsons_data[-train, ]$status == pred), 5) *100

#conf_mat <- confusionMatrix(factor(pred), factor(parkinsons_data[-train, ]$status))

precision_log <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['0', '1']), 5) * 100

recall_log <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['1', '0']), 5) * 100

glue("After the evaluation of logistic regression model the accuracy, precision and recall metrics are the following respectively {accuracy_log}%, {precision_log}% and {recall_log}%")

## After the evaluation of logistic regression model the accuracy, precision and recall metrics are the following respectively 32.727%, 60% and 41.86%

models <- c("Logistic regression")
accuracy <- c(accuracy_log)
precision <- c(precision_log)
recall <- c(recall_log)

Evaluating on people the model has never seen

An earlier version of this analysis split the data by recording (sample(1:195, 140)). Because each person contributes about 6 recordings that sound alike, that split leaks: the same person appears in both the training and test sets, and the models are partly rewarded for recognising a voice rather than the disease. Under that random split, plain logistic regression reached 85.5% accuracy and the decision tree 94.5%. With whole people held out, as above, the picture is very different, and this is the realistic scenario for a screening tool.

Two things are worth keeping in mind for every result below:

Baseline. 43 of the 55 test recordings (78.2%) are from people with Parkinson’s, so labelling every recording as Parkinson’s already gives 78.2% accuracy, 78.2% precision and 100% recall. A model has to beat that to be useful.
Small test set. The test set contains only 9 people (2 of them healthy), so a single split is very noisy. The Robustness check section repeats the split 30 times.

For plain logistic regression the held-out accuracy is only 32.7%, well below the baseline. With 22 strongly collinear predictors and only 23 training people, the unregularised model overfits, and the significance and VIF results above already warned that its coefficients are unstable.

Question 2:

Parametric

The parametric model we are going to evaluate are: 1) Logistic regression (already covered in Question 1) 2) linear and Quadratic discriminant analysis. 3) Naive Bayes 4) Lasso and Ridge regularization.

set.seed(32)

#parkinsons_data$status <- as.factor(parkinsons_data$status)

#lda_model <- lda(status ~ ., data = parkinsons_data[train, -1])

#parkinsons_data$MDVP.Jitter... <- scale(parkinsons_data$MDVP.Jitter...)

lda_model <- lda(status ~ ., data = parkinsons_data[train, -c(1,5, 4, 3, 2, 6)])

pred <- predict(lda_model, newdata = parkinsons_data[-train, -1])

conf_mat <- conf_table(parkinsons_data[-train, 18], pred$class)

conf_mat

##    
##      0  1
##   0  2 10
##   1  6 37

accuracy_lda <- round(mean(parkinsons_data[-train, 18] == pred$class), 5) *100

precision_lda <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['0', '1']), 5) * 100

recall_lda <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['1', '0']), 5) * 100

glue("After the evaluation of linear Discriminant analysis regression model the accuracy, precision and recall metrics are the following respectively {accuracy_lda}%, {precision_lda}% and {recall_lda}%")

## After the evaluation of linear Discriminant analysis regression model the accuracy, precision and recall metrics are the following respectively 70.909%, 78.723% and 86.047%

models <- c(models, "LDA model")
accuracy <- c(accuracy, accuracy_lda)
precision <- c(precision, precision_lda)
recall <- c(recall, recall_lda)

Note

The model appears to improve over logistic regression, but this result isn’t fully reliable. Some features had to be removed because they showed no variation within one or more classes, which is a requirement for Linear Discriminant Analysis to work properly.

set.seed(32)

qda_model <- qda(status ~., data = parkinsons_data[train, -1])

pred <- predict(qda_model, newdata = parkinsons_data[-train, -1])

conf_mat <- conf_table(parkinsons_data[-train, 18], pred$class)

accuracy_qda <- round(mean(parkinsons_data[-train, 18] == pred$class), 5) * 100

precision_qda <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['0', '1']), 5) * 100

recall_qda <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['1', '0']), 5) * 100

glue("After the evaluation of Quadratic Discriminant analysis regression model the accuracy, precision and recall metrics are the following respectively {accuracy_qda}%, {precision_qda}% and {recall_qda}%")

## After the evaluation of Quadratic Discriminant analysis regression model the accuracy, precision and recall metrics are the following respectively 78.182%, 78.182% and 100%

models <- c(models, "QDA model")
accuracy <- c(accuracy, accuracy_qda)
precision <- c(precision, precision_qda)
recall <- c(recall, recall_qda)

QDA reaches 100% recall, but this is misleading: its accuracy and precision (78.2%) are exactly the always-Parkinson’s baseline, so it labelled every test recording as Parkinson’s and missed all 12 healthy recordings. A model that never says “healthy” is not useful for screening even though its recall looks perfect.

set.seed(32)

naiveBayes_model <- naiveBayes(status ~ ., data = parkinsons_data[train, -1])

pred <- predict(naiveBayes_model, newdata = parkinsons_data[-train, -1])

conf_mat <- conf_table(parkinsons_data[-train, 18], pred)

accuracy_nB <- round(mean(parkinsons_data[-train, 18] == pred), 5) * 100

precision_nB <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['0', '1']), 5) * 100

recall_nB <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['1', '0']), 5) * 100

glue("After the evaluation of Naive Bayes analysis regression model the accuracy, precision and recall metrics are the following respectively {accuracy_nB}%, {precision_nB}% and {recall_nB}%")

## After the evaluation of Naive Bayes analysis regression model the accuracy, precision and recall metrics are the following respectively 87.273%, 90.909% and 93.023%

models <- c(models, "Naive Bayes")
accuracy <- c(accuracy, accuracy_nB)
precision <- c(precision, precision_nB)
recall <- c(recall, recall_nB)

Naive Bayes is the top model on this particular split (87.3% accuracy, 90.9% precision, 93.0% recall), but that should not be over-interpreted: it is a single test set of 9 people, and in the repeated-split check further down Naive Bayes turns out to be one of the two worst models on average. This is a good illustration of why one split cannot be used to pick a winner.

x_train <- model.matrix(status ~ ., data = parkinsons_data[train, -1])[,-1]

y_train <- parkinsons_data[train, ]$status

lasso_model <- cv.glmnet(x_train, y_train, alpha = 1, family = "binomial")

x_test <- model.matrix(status ~ ., data = parkinsons_data[-train, -1])[, -1]

prob <- predict(lasso_model, newx = x_test, s = lasso_model$lambda.min, type = "response")

pred <- ifelse(prob > 0.5, 1, 0)


conf_mat <- conf_table(parkinsons_data[-train, 18], pred)

accuracy_l1 <- round(mean(parkinsons_data[-train, ]$status == pred), 5) *100

precision_l1 <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['0', '1']), 5) * 100

recall_l1 <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['1', '0']), 5) * 100

glue("After the evaluation of lasso regularized model the accuracy, precision and recall metrics are the following respectively {accuracy_l1}%, {precision_l1}% and {recall_l1}%")

## After the evaluation of lasso regularized model the accuracy, precision and recall metrics are the following respectively 69.091%, 76% and 88.372%

models <- c(models, "Lasso regularization")
accuracy <- c(accuracy, accuracy_l1)
precision <- c(precision, precision_l1)
recall <- c(recall, recall_l1)





x_train <- model.matrix(status ~ ., data = parkinsons_data[train, -1])[,-1]

y_train <- parkinsons_data[train, ]$status

ridge_model <- cv.glmnet(x_train, y_train, alpha = 0, family = "binomial")

x_test <- model.matrix(status ~ ., data = parkinsons_data[-train, -1])[, -1]

prob <- predict(ridge_model, newx = x_test, s = ridge_model$lambda.min, type = "response")

pred <- ifelse(prob > 0.5, 1, 0)


conf_mat <- conf_table(parkinsons_data[-train, 18], pred)

accuracy_l2 <- round(mean(parkinsons_data[-train, ]$status == pred), 5) *100

precision_l2 <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['0', '1']), 5) * 100

recall_l2 <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['1', '0']), 5) * 100

glue("After the evaluation of Ridge Regularized model the accuracy, precision and recall metrics are the following respectively {accuracy_l2}%, {precision_l2}% and {recall_l2}%")

## After the evaluation of Ridge Regularized model the accuracy, precision and recall metrics are the following respectively 78.182%, 78.182% and 100%

models <- c(models, "Ridge Regularization")
accuracy <- c(accuracy, accuracy_l2)
precision <- c(precision, precision_l2)
recall <- c(recall, recall_l2)

Regularisation helps a lot compared with plain logistic regression here (Lasso 69.1% and Ridge 78.2% accuracy versus 32.7%), which is what we expect when there are many collinear predictors and few training people. Ridge’s 78.2% accuracy and 100% recall equal the always-Parkinson’s baseline, so it again predicted every recording as Parkinson’s, while Lasso is below that baseline.

Non-parametric models

KNN
Decision Trees
Random Forrest
Boosting

#using caret library.

set.seed(32)

parkinsons_data$status <- as.factor(parkinsons_data$status)

ctrl <- trainControl(method = "cv", number = 10)



knn_model <- train(status ~ ., data = parkinsons_data[train , -1], method = "knn",
                   preProcess = c("center", "scale"), trControl = ctrl)

pred <- predict(knn_model, newdata = parkinsons_data[-train, -1])

conf_mat <- conf_table(parkinsons_data[-train, 18], pred)

accuracy_knn <- round(mean(parkinsons_data[-train, ]$status == pred), 5) *100

precision_knn <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['0', '1']), 5) * 100

recall_knn <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['1', '0']), 5) * 100

glue("After the evaluation of knn model the accuracy, precision and recall metrics are the following respectively {accuracy_knn}%, {precision_knn}% and {recall_knn}%")

## After the evaluation of knn model the accuracy, precision and recall metrics are the following respectively 70.909%, 76.471% and 90.698%

models <- c(models, "KNN")
accuracy <- c(accuracy, accuracy_knn)
precision <- c(precision, precision_knn)
recall <- c(recall, recall_knn)

KNN (70.9% accuracy, 76.5% precision, 90.7% recall, on standardised features with k chosen by 10-fold cross-validation) is below the always-Parkinson’s baseline.

set.seed(32)

tree_model <- tree(status ~ ., data = parkinsons_data[train, -1])

pred <- predict(tree_model, parkinsons_data[-train, -1], type = "class")

conf_mat <- conf_table(parkinsons_data[-train , 18], pred)

accuracy_tree <- round(mean(parkinsons_data[-train, ]$status == pred), 5) *100

precision_tree <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['0', '1']), 5) * 100

recall_tree <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['1', '0']), 5) * 100

glue("After the evaluation of decision tree model the accuracy, precision and recall metrics are the following respectively {accuracy_tree}%, {precision_tree}% and {recall_tree}%")

## After the evaluation of decision tree model the accuracy, precision and recall metrics are the following respectively 74.545%, 78.431% and 93.023%

models <- c(models, "Decision Tree")
accuracy <- c(accuracy, accuracy_tree)
precision <- c(precision, precision_tree)
recall <- c(recall, recall_tree)

The single decision tree (74.5% accuracy, 78.4% precision, 93.0% recall) is below the always-Parkinson’s baseline on this split. Under the earlier recording-level split it had looked like the best model (94.5% accuracy), which suggests it was largely memorising individual voices.

set.seed(32)

rf_model <- randomForest(status ~ ., data = parkinsons_data[train, -1], ntree = 500)


pred <- unname(predict(rf_model, newdata = parkinsons_data[-train, -1], type = "class"))

conf_mat <- conf_table(parkinsons_data[-train, 18], pred)

accuracy_rf <- round(mean(parkinsons_data[-train, ]$status == pred), 5) *100

precision_rf <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['0', '1']), 5) * 100

recall_rf <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['1', '0']), 5) * 100

glue("After the evaluation of random forrest model the accuracy, precision and recall metrics are the following respectively {accuracy_rf}%, {precision_rf}% and {recall_rf}%")

## After the evaluation of random forrest model the accuracy, precision and recall metrics are the following respectively 78.182%, 78.182% and 100%

models <- c(models, "Random Forrest")
accuracy <- c(accuracy, accuracy_rf)
precision <- c(precision, precision_rf)
recall <- c(recall, recall_rf)





set.seed(32)

parkinsons_data$status <- as.numeric(as.character(parkinsons_data$status))

boosting_model <- gbm(status ~ ., data = parkinsons_data[train, -1], distribution = "bernoulli", n.trees = 195, cv.folds = 10)

prob <- predict(boosting_model, newdata = parkinsons_data[-train, -1], type = "response")

pred <- ifelse(prob > 0.5, 1, 0)

conf_mat <- conf_table(parkinsons_data[-train, 18], pred)

accuracy_bst <- round(mean(parkinsons_data[-train, ]$status == pred), 5) *100

precision_bst <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['0', '1']), 5) * 100

recall_bst <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['1', '0']), 5) * 100

glue("After the evaluation of boosting model the accuracy, precision and recall metrics are the following respectively {accuracy_bst}%, {precision_bst}% and {recall_bst}%")

## After the evaluation of boosting model the accuracy, precision and recall metrics are the following respectively 78.182%, 78.182% and 100%

models <- c(models, "Boosting")
accuracy <- c(accuracy, accuracy_bst)
precision <- c(precision, precision_bst)
recall <- c(recall, recall_bst)

Random forest and boosting both end up at exactly the baseline numbers (78.2% / 78.2% / 100%), so on this split neither learned to recognise healthy voices from unseen people. Tree ensembles can memorise the voices of the people in the training set very well, which does not transfer to new people when there are only 23 of them.

Support Vector Machine

set.seed(32)

parkinsons_data$status <- as.factor(parkinsons_data$status)

svm_tune <- tune(svm, status ~ ., data = parkinsons_data[train, -1], kernel = "radial",
                 ranges = list(cost = c(0.01, 0.1, 0.5, 1, 2, 4, 8, 10),
                               gamma = c(0.01, 0.1, 0.5, 1, 2, 4, 8, 10)))

glue("Best SVM parameters from 10-fold cross-validation on the training set: cost = {svm_tune$best.parameters$cost}, gamma = {svm_tune$best.parameters$gamma}")

## Best SVM parameters from 10-fold cross-validation on the training set: cost = 4, gamma = 0.1

svm_best <- svm_tune$best.model

pred <- predict(svm_best, newdata = parkinsons_data[-train, -c(1, 18)])

conf_mat <- conf_table(parkinsons_data[-train, 18], pred)

accuracy_svm <- round(mean(parkinsons_data[-train, ]$status == pred), 5) *100

precision_svm <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['0', '1']), 5) * 100

recall_svm <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['1', '0']), 5) * 100

glue("After the evaluation of SVM model the accuracy, precision and recall metrics are the following respectively {accuracy_svm}%, {precision_svm}% and {recall_svm}%")

## After the evaluation of SVM model the accuracy, precision and recall metrics are the following respectively 76.364%, 77.778% and 97.674%

models <- c(models, "Support Vector Machine")
accuracy <- c(accuracy, accuracy_svm)
precision <- c(precision, precision_svm)
recall <- c(recall, recall_svm)

We tuned a radial-kernel Support Vector Machine over a grid of cost and gamma values using 10-fold cross-validation on the training set. The best combination was a cost of 4 and gamma of 0.1. Its test performance (76.4% accuracy, 77.8% precision, 97.7% recall) is slightly below the baseline.

Although the model’s performance in terms of accuracy and precision was not particularly strong, it stood out in terms of recall (sensitivity) — effectively identifying the positive cases, which is especially important in medical data analysis.

Question 3

Regularization covered in question 2, now performing dimension reduction using PCA

pca_model <- prcomp(parkinsons_data[train, -c(1, 18)], scale. = TRUE)

plot(1:22, summary(pca_model)[[6]][2, ], type = 'l', xlab = "Prin Comp", 
     ylab = "% Var Explained")

print("From the scree plot, we can observe an elbow around the 5th principal component, suggesting it as an optimal cutoff. The first five components together explain about 90% of the variance in the training data, making them a reasonable choice for dimensionality reduction.")

## [1] "From the scree plot, we can observe an elbow around the 5th principal component, suggesting it as an optimal cutoff. The first five components together explain about 90% of the variance in the training data, making them a reasonable choice for dimensionality reduction."

# Best number of principal components would be 5 and it explains about 90% of the variation.

pca_train <- predict(pca_model, newdata = parkinsons_data[train, -c(1, 18)], type = "response")[,1:5]
pca_test <- predict(pca_model, newdata = parkinsons_data[-train, -c(1, 18)], type = "response")[,1:5]

pca_df <- as.data.frame(pca_train[,1:5])

colnames(pca_df) <- paste0("PC", 1:5)

pca_df$status <- parkinsons_data[train, 18]

pca_logit_model <- glm(status ~ ., data = pca_df, family = "binomial")

pca_test_df <- as.data.frame(pca_test)

prob <- predict(pca_logit_model, newdata = pca_test_df, type = "response")

pred <- ifelse(prob > 0.5, 1, 0)

conf_mat <- conf_table(parkinsons_data[-train, 18], pred)

accuracy_pca <- round((conf_mat['1', '1']+conf_mat['0','0'])/(conf_mat['1', '1']+conf_mat['0','0'] + conf_mat['1','0']+conf_mat['0', '1']), 5) *100

precision_pca <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['0', '1']), 5) * 100

recall_pca <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['1', '0']), 5) * 100

glue("After the evaluation of PCA model the accuracy, precision and recall metrics are the following respectively {accuracy_pca}%, {precision_pca}% and {recall_pca}%")

## After the evaluation of PCA model the accuracy, precision and recall metrics are the following respectively 81.818%, 81.132% and 100%

models <- c(models, "Principal Component regression")
accuracy <- c(accuracy, accuracy_pca)
precision <- c(precision, precision_pca)
recall <- c(recall, recall_pca)

data_metrics <- data.frame(models = models, accuracy = accuracy, precision = precision, recall = recall)

data_metrics[models %in% c("Ridge Regularization", "Lasso regularization", "Principal Component regression", "Logistic regression"), ]

##                            models accuracy precision  recall
## 1             Logistic regression   32.727    60.000  41.860
## 5            Lasso regularization   69.091    76.000  88.372
## 6            Ridge Regularization   78.182    78.182 100.000
## 12 Principal Component regression   81.818    81.132 100.000

Lasso, Ridge and PCA all improve substantially on plain logistic regression (32.7% accuracy), because they each tame the multicollinearity that makes the unregularised model unstable. PCA + logistic regression is the best of the four on this split (81.8% accuracy, 81.1% precision, 100% recall), and unlike Ridge it does correctly identify a few of the healthy recordings (2 of 12).

Question 2 (continued): Comparing all models

data_metrics <- data.frame(models = models, accuracy = accuracy, precision = precision, recall = recall)

print(data_metrics, row.names = FALSE)

##                          models accuracy precision  recall
##             Logistic regression   32.727    60.000  41.860
##                       LDA model   70.909    78.723  86.047
##                       QDA model   78.182    78.182 100.000
##                     Naive Bayes   87.273    90.909  93.023
##            Lasso regularization   69.091    76.000  88.372
##            Ridge Regularization   78.182    78.182 100.000
##                             KNN   70.909    76.471  90.698
##                   Decision Tree   74.545    78.431  93.023
##                  Random Forrest   78.182    78.182 100.000
##                        Boosting   78.182    78.182 100.000
##          Support Vector Machine   76.364    77.778  97.674
##  Principal Component regression   81.818    81.132 100.000

On this single split, Naive Bayes (87.3% accuracy) and PCA regression (81.8%) are the only two models above the 78.2% always-Parkinson’s baseline. QDA, Ridge, Random Forest and Boosting all sit exactly on the baseline, and Logistic regression, LDA, Lasso, KNN, the Decision Tree and the SVM are below it. Because the test set holds only 9 people, differences of a few percentage points correspond to one or two recordings and should not be over-read.

Recall alone is not a good yardstick here: four models (QDA, Ridge, Random Forest and Boosting) reach 100% recall simply by never predicting “healthy”. Precision and accuracy have to be read alongside it, and the repeated-split check below gives a more trustworthy ranking than this single split.

results_melted <- melt(data_metrics, id.vars = "models")

ggplot(results_melted, aes(x = models, y = value, fill = variable)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(title = "Model Comparison: Accuracy, Precision, Recall",
       y = "Percentage",
       x = "Model",
       fill = "Metric") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_y_continuous(labels = scales::percent_format(scale = 1))

There is no consistent advantage for parametric or non-parametric models on this dataset: the best and worst models on this split come from both families, and the ordering changes from split to split (see the robustness check).

Question 5

set.seed(32)

knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)


ctrl <- trainControl(method = "LOOCV")

parkinsons_data$status <- as.factor(parkinsons_data$status)

logit_model <- train(status ~ ., data = parkinsons_data[train, -1], method = "glm", family = "binomial", trControl = ctrl)

prob <- predict(logit_model, newdata = parkinsons_data[-train, -c(1, 18)])

pred <- ifelse(prob > 0.5, 1, 0)

conf_mat <- conf_table(parkinsons_data[-train, 18], prob)

accuracy_log_loocv <- round(mean(parkinsons_data[-train, ]$status == prob), 5) *100

precision_log_loocv <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['0', '1']), 5) * 100

recall_log_loocv <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['1', '0']), 5) * 100

glue("After the evaluation of LOOCV logistic regression the accuracy, precision and recall metrics are the following respectively {accuracy_log_loocv}%, {precision_log_loocv}% and {recall_log_loocv}%")

## After the evaluation of LOOCV logistic regression the accuracy, precision and recall metrics are the following respectively 32.727%, 60% and 41.86%

set.seed(32)

ctrl <- trainControl(method = "cv", number = 10)


logit_cv_model <- train(status ~ ., data = parkinsons_data[train, -1], method = "glm",family = "binomial",
                        trControl = ctrl)

prob <- predict(logit_cv_model, newdata = parkinsons_data[-train, -c(1,18)])

pred <- ifelse(prob > 0.5, 1, 0)

conf_mat <- conf_table(parkinsons_data[-train, 18], prob)

accuracy_log_10cv <- round(mean(parkinsons_data[-train, ]$status == prob), 5) *100

precision_log_10cv <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['0', '1']), 5) * 100

recall_log_10cv <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['1', '0']), 5) * 100

glue("After the evaluation of 10-fold CV logistic regression the accuracy, precision and recall metrics are the following respectively {accuracy_log_10cv}%, {precision_log_10cv}% and {recall_log_10cv}%")

## After the evaluation of 10-fold CV logistic regression the accuracy, precision and recall metrics are the following respectively 32.727%, 60% and 41.86%

Cross-validating logistic regression gives exactly the same test metrics as before (32.7% accuracy, 60.0% precision, 41.9% recall), and LOOCV and 10-fold agree. That is expected: logistic regression has no tuning parameters, so `caret` uses cross-validation only to estimate performance and then refits the same model on the full training set.

set.seed(32)

ctrl <- trainControl(method = "loocv")

qda_cv_model <- train(status ~ ., data = parkinsons_data[train, -1], method = "qda", trControl = ctrl)

pred <- predict(qda_cv_model, newdata = parkinsons_data[-train, -c(1, 18)])

conf_mat <- conf_table(parkinsons_data[-train, 18], pred)

accuracy_qda_loocv <- round(mean(parkinsons_data[-train, ]$status == pred), 5) *100

precision_qda_loocv <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['0', '1']), 5) * 100

recall_qda_loocv <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['1', '0']), 5) * 100

glue("After the evaluation LOOCV in QDA model, the accuracy, precision and recall metrics are the following respectively {accuracy_qda_loocv}%, {precision_qda_loocv}% and {recall_qda_loocv}%")

## After the evaluation LOOCV in QDA model, the accuracy, precision and recall metrics are the following respectively 78.182%, 78.182% and 100%

set.seed(32)

ctrl10 <- trainControl(method = "cv", number = 10)

qda_cv_model <- train(status ~ ., data = parkinsons_data[train, -1], method = "qda", trControl = ctrl)

pred <- predict(qda_cv_model, newdata = parkinsons_data[-train, -c(1, 18)])

conf_mat <- conf_table(parkinsons_data[-train, 18], pred)

accuracy_qda_10cv <- round(mean(parkinsons_data[-train, ]$status == pred), 5) *100

precision_qda_10cv <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['0', '1']), 5) * 100

recall_qda_10cv <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['1', '0']), 5) * 100

glue("After the evaluation of 10-fold CV in QDA model, the accuracy, precision and recall metrics are the following respectively {accuracy_qda_10cv}%, {precision_qda_10cv}% and {recall_qda_10cv}%")

## After the evaluation of 10-fold CV in QDA model, the accuracy, precision and recall metrics are the following respectively 78.182%, 78.182% and 100%

Cross-validation leaves the QDA results unchanged for the same reason: QDA has no tuning parameter, so the final model is identical to the one fitted before.

set.seed(32)

ctrl <- trainControl(method = "loocv")

tree_loocv_model <- train(status ~ ., data = parkinsons_data[train, -1], method = "rpart", trControl = ctrl)

pred <- predict(tree_loocv_model, newdata = parkinsons_data[-train, -c(1, 18)])


conf_mat <- conf_table(parkinsons_data[-train, 18], pred)

accuracy_tree_loocv <- round(mean(parkinsons_data[-train, ]$status == pred), 5) *100

precision_tree_loocv <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['0', '1']), 5) * 100

recall_tree_loocv <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['1', '0']), 5) * 100

glue("After the evaluation of LOOCV on Decision tree model, the accuracy, precision and recall metrics are the following respectively {accuracy_tree_loocv}%, {precision_tree_loocv}% and {recall_tree_loocv}%")

## After the evaluation of LOOCV on Decision tree model, the accuracy, precision and recall metrics are the following respectively 80%, 79.63% and 100%

set.seed(32)

ctrl <- trainControl(method = "cv", number = 10)

tree_10cv_model <- train(status ~ ., data = parkinsons_data[train, -1], method = "rpart", trControl = ctrl)

pred <- predict(tree_10cv_model, newdata = parkinsons_data[-train, -c(1, 18)])


conf_mat <- conf_table(parkinsons_data[-train, 18], pred)

accuracy_tree_10cv <- round(mean(parkinsons_data[-train, ]$status == pred), 5) *100

precision_tree_10cv <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['0', '1']), 5) * 100

recall_tree_10cv <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['1', '0']), 5) * 100

glue("After the evaluation of 10-fold CV on Decision tree model, the accuracy, precision and recall metrics are the following respectively {accuracy_tree_10cv}%, {precision_tree_10cv}% and {recall_tree_10cv}%")

## After the evaluation of 10-fold CV on Decision tree model, the accuracy, precision and recall metrics are the following respectively 80%, 79.63% and 100%

For the decision tree, cross-validation (through `caret`’s `rpart`, which tunes the complexity parameter) gives 80.0% accuracy, 79.6% precision and 100% recall with both LOOCV and 10-fold, slightly better than the untuned `tree` model (74.5%). Part of that difference comes from using a different tree implementation with a tuned pruning parameter. The recall of 100% again mostly reflects predicting Parkinson’s for nearly everyone.

ctrl <- trainControl(method = "LOOCV")

parkinsons_data$status <- as.factor(parkinsons_data$status)

pca_logit_model_loocv <- train(status ~ ., data = pca_df, method = "glm", family = "binomial", trControl = ctrl)

prob <- predict(pca_logit_model_loocv, newdata = pca_test_df, type = "raw")

pred <- prob

conf_mat <- conf_table(parkinsons_data[-train, 18], pred)

accuracy_pca_loocv <- round((conf_mat['1', '1']+conf_mat['0','0'])/(conf_mat['1', '1']+conf_mat['0','0'] + conf_mat['1','0']+conf_mat['0', '1']), 5) *100

precision_pca_loocv <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['0', '1']), 5) * 100

recall_pca_loocv <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['1', '0']), 5) * 100

glue("After the evaluation of LOOCV PCA logistic regression the accuracy, precision and recall metrics are the following respectively {accuracy_pca_loocv}%, {precision_pca_loocv}% and {recall_pca_loocv}%")

## After the evaluation of LOOCV PCA logistic regression the accuracy, precision and recall metrics are the following respectively 81.818%, 81.132% and 100%

set.seed(32)

ctrl <- trainControl(method = "cv", number = 10)


pca_logit_10cv <- train(status ~ ., data = pca_df, method = "glm",family = "binomial",
                        trControl = ctrl)

prob <- predict(pca_logit_10cv, newdata = pca_test_df)

pred <- prob

conf_mat <- conf_table(parkinsons_data[-train, 18], pred)

accuracy_pca_10cv <- round((conf_mat['1', '1']+conf_mat['0','0'])/(conf_mat['1', '1']+conf_mat['0','0'] + conf_mat['1','0']+conf_mat['0', '1']), 5) *100

precision_pca_10cv <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['0', '1']), 5) * 100

recall_pca_10cv <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['1', '0']), 5) * 100

glue("After the evaluation of 10-fold CV PCA logistic regression the accuracy, precision and recall metrics are the following respectively {accuracy_pca_10cv}%, {precision_pca_10cv}% and {recall_pca_10cv}%")

## After the evaluation of 10-fold CV PCA logistic regression the accuracy, precision and recall metrics are the following respectively 81.818%, 81.132% and 100%

Cross-validating the PCA logistic regression also leaves the results unchanged under both LOOCV and 10-fold (81.8% accuracy, 81.1% precision, 100% recall), again because the model has no tuning parameter to select.

set.seed(32)

ctrl <- trainControl(method = "LOOCV")

parkinsons_data$status <- as.factor(parkinsons_data$status)

rf_loocv <- train(status ~ ., data = parkinsons_data[train, -1], method = "rf", trControl = ctrl)

pred <- predict(rf_loocv, newdata = parkinsons_data[-train, -c(1, 18)])

conf_mat <- conf_table(parkinsons_data[-train, 18], pred)

accuracy_rf_loocv <- round((conf_mat['1', '1']+conf_mat['0','0'])/(conf_mat['1', '1']+conf_mat['0','0'] + conf_mat['1','0']+conf_mat['0', '1']), 5) *100

precision_rf_loocv <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['0', '1']), 5) * 100

recall_rf_loocv <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['1', '0']), 5) * 100

glue("After the evaluation of LOOCV random forest the accuracy, precision and recall metrics are the following respectively {accuracy_rf_loocv}%, {precision_rf_loocv}% and {recall_rf_loocv}%")

## After the evaluation of LOOCV random forest the accuracy, precision and recall metrics are the following respectively 76.364%, 77.778% and 97.674%

set.seed(32)

ctrl <- trainControl(method = "cv", number = 10)

rf_10cv <- train(status ~ ., data = parkinsons_data[train, -1], method = "rf", trControl = ctrl)

pred <- predict(rf_10cv, newdata = parkinsons_data[-train, -c(1, 18)])

conf_mat <- conf_table(parkinsons_data[-train, 18], pred)

accuracy_rf_10cv <- round((conf_mat['1', '1']+conf_mat['0','0'])/(conf_mat['1', '1']+conf_mat['0','0'] + conf_mat['1','0']+conf_mat['0', '1']), 5) *100

precision_rf_10cv <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['0', '1']), 5) * 100

recall_rf_10cv <- round(conf_mat['1', '1']/(conf_mat['1', '1'] + conf_mat['1', '0']), 5) * 100

glue("After the evaluation of 10-fold CV random forest the accuracy, precision and recall metrics are the following respectively {accuracy_rf_10cv}%, {precision_rf_10cv}% and {recall_rf_10cv}%")

## After the evaluation of 10-fold CV random forest the accuracy, precision and recall metrics are the following respectively 78.182%, 78.182% and 100%

For the random forest, caret tunes mtry by cross-validation. With LOOCV the results move slightly down (76.4% accuracy, 77.8% precision, 97.7% recall), and with 10-fold CV they equal the earlier untuned forest (78.2% / 78.2% / 100%). These differences amount to one test recording.

before <- function(m) c(accuracy[models == m], precision[models == m], recall[models == m])

before_mat <- rbind(
  "Logistic"       = before("Logistic regression"),
  "Decision Tree"  = before("Decision Tree"),
  "Random Forrest" = before("Random Forrest"),
  "QDA"            = before("QDA model"),
  "PCA+Logit"      = before("Principal Component regression")
)

# "After" = the LOOCV-trained versions of the same models
after_mat <- rbind(
  "Logistic"       = c(accuracy_log_loocv, precision_log_loocv, recall_log_loocv),
  "Decision Tree"  = c(accuracy_tree_loocv, precision_tree_loocv, recall_tree_loocv),
  "Random Forrest" = c(accuracy_rf_loocv, precision_rf_loocv, recall_rf_loocv),
  "QDA"            = c(accuracy_qda_loocv, precision_qda_loocv, recall_qda_loocv),
  "PCA+Logit"      = c(accuracy_pca_loocv, precision_pca_loocv, recall_pca_loocv)
)

df <- data.frame(
  Model = rownames(before_mat),
  Accuracy_Before = before_mat[, 1],
  Accuracy_After  = after_mat[, 1],
  Precision_Before = before_mat[, 2],
  Precision_After  = after_mat[, 2],
  Recall_Before = before_mat[, 3],
  Recall_After  = after_mat[, 3]
)

print(df, row.names = FALSE)

##           Model Accuracy_Before Accuracy_After Precision_Before Precision_After
##        Logistic          32.727         32.727           60.000          60.000
##   Decision Tree          74.545         80.000           78.431          79.630
##  Random Forrest          78.182         76.364           78.182          77.778
##             QDA          78.182         78.182           78.182          78.182
##       PCA+Logit          81.818         81.818           81.132          81.132
##  Recall_Before Recall_After
##         41.860       41.860
##         93.023      100.000
##        100.000       97.674
##        100.000      100.000
##        100.000      100.000

delta_df <- df %>%
  mutate(
    Accuracy_Diff = Accuracy_After - Accuracy_Before,
    Precision_Diff = Precision_After - Precision_Before,
    Recall_Diff = Recall_After - Recall_Before
  ) %>%
  dplyr::select(Model, Accuracy_Diff, Precision_Diff, Recall_Diff) %>%
  pivot_longer(
    cols = -Model,
    names_to = "Metric",
    values_to = "Delta"
  )


delta_df$Metric <- gsub("_Diff", "", delta_df$Metric)


ggplot(delta_df, aes(x = Model, y = Delta, fill = Metric)) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(
    title = "Change in Metrics After Cross-Validation",
    y = "Change (%)",
    x = "Model"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Cross-validation changes very little: the models with no tuning parameters (logistic, QDA, PCA + logistic) are unchanged by construction, and the two models with tuned parameters (decision tree and random forest) move by 2 to 5 accuracy points, which on a 55-recording test set is 1 to 3 recordings. These changes are within the noise of a single split and do not support conclusions about which model generalises better. One further caveat: the folds used for tuning are random over recordings, so recordings from the same person can land in both the fitting and validation folds; grouping the folds by person would be more appropriate.

Robustness check: repeated subject-level splits

The test set above contains only 9 people, 2 of whom are healthy, so a single split can say very little about which model is really best. To check how stable the comparison is, the subject-level split is repeated 30 times with different random seeds and the models that are fast to fit are re-evaluated on each split (SVM, KNN and boosting are left out to keep the run time short). A baseline that simply labels every recording as Parkinson’s is included for reference. The table shows the mean and standard deviation of each metric across splits.

score <- function(actual, pred) {
  tp <- sum(actual == 1 & pred == 1)
  fp <- sum(actual == 0 & pred == 1)
  fn <- sum(actual == 1 & pred == 0)
  c(Accuracy = mean(actual == pred) * 100,
    Precision = if (tp + fp > 0) tp / (tp + fp) * 100 else NA,
    Recall = tp / (tp + fn) * 100)
}

as01 <- function(x) as.numeric(as.character(x))

robust_split <- function(seed) {
  set.seed(seed)
  test_sub <- c(sample(names(subject_status)[subject_status == 0], 2),
                sample(names(subject_status)[subject_status == 1], 7))
  tr <- which(!subject %in% test_sub)

  d <- parkinsons_data[, -1]
  d$status <- as01(d$status)
  d_train <- d[tr, ]
  d_test <- d[-tr, ]
  f_train <- d_train
  f_train$status <- factor(f_train$status)
  x_cols <- setdiff(names(d), "status")

  pc <- prcomp(d_train[, x_cols], scale. = TRUE)
  pc_train <- data.frame(pc$x[, 1:5], status = d_train$status)
  pc_test <- as.data.frame(predict(pc, d_test[, x_cols])[, 1:5])

  preds <- list(
    "Baseline: always Parkinson's" = function() rep(1, nrow(d_test)),
    "Logistic regression" = function() as.numeric(predict(glm(status ~ ., data = d_train, family = "binomial"), d_test, type = "response") > 0.5),
    "LDA" = function() as01(predict(lda(status ~ ., data = d_train[, -(1:5)]), d_test)$class),
    "QDA" = function() as01(predict(qda(status ~ ., data = d_train), d_test)$class),
    "Naive Bayes" = function() as01(predict(naiveBayes(status ~ ., data = f_train), d_test)),
    "PCA + Logistic" = function() as.numeric(predict(glm(status ~ ., data = pc_train, family = "binomial"), pc_test, type = "response") > 0.5),
    "Decision Tree" = function() as01(predict(tree(status ~ ., data = f_train), d_test, type = "class")),
    "Random Forest" = function() as01(predict(randomForest(status ~ ., data = f_train, ntree = 200), d_test))
  )

  do.call(rbind, lapply(names(preds), function(m) {
    sc <- tryCatch(score(d_test$status, preds[[m]]()), error = function(e) c(Accuracy = NA, Precision = NA, Recall = NA))
    data.frame(Model = m, seed = seed, t(sc))
  }))
}

robust_results <- do.call(rbind, lapply(1:30, robust_split))

robust_summary <- robust_results %>%
  group_by(Model) %>%
  summarise(
    mean_accuracy = mean(Accuracy, na.rm = TRUE),
    Accuracy = sprintf("%.1f ± %.1f", mean(Accuracy, na.rm = TRUE), sd(Accuracy, na.rm = TRUE)),
    Precision = sprintf("%.1f ± %.1f", mean(Precision, na.rm = TRUE), sd(Precision, na.rm = TRUE)),
    Recall = sprintf("%.1f ± %.1f", mean(Recall, na.rm = TRUE), sd(Recall, na.rm = TRUE)),
    .groups = "drop"
  ) %>%
  arrange(desc(mean_accuracy)) %>%
  dplyr::select(-mean_accuracy)

print(as.data.frame(robust_summary), row.names = FALSE)

##                         Model    Accuracy  Precision      Recall
##                PCA + Logistic  84.0 ± 8.8 89.5 ± 6.7  90.7 ± 8.5
##                 Random Forest  82.6 ± 6.9 86.8 ± 7.2  92.7 ± 6.5
##                           QDA  82.3 ± 5.0 83.0 ± 4.4  97.6 ± 2.8
##                           LDA  81.2 ± 8.4 89.2 ± 6.8  87.1 ± 9.5
##                 Decision Tree  78.7 ± 7.0 85.8 ± 6.5  87.9 ± 7.0
##  Baseline: always Parkinson's  78.2 ± 0.3 78.2 ± 0.3 100.0 ± 0.0
##           Logistic regression 71.5 ± 10.2 84.7 ± 7.5 78.4 ± 13.5
##                   Naive Bayes 71.3 ± 10.5 93.1 ± 7.0 68.6 ± 12.0

ggplot(melt(robust_results[, c("Model", "Accuracy", "Precision", "Recall")], id.vars = "Model"),
       aes(x = Model, y = value, fill = variable)) +
  geom_boxplot() +
  labs(title = "Metrics across 30 random subject-level splits", y = "Percentage", x = "Model", fill = "Metric") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Averaged over 30 random subject-level splits, PCA + logistic regression has the highest mean accuracy (about 84%, with about 90% precision and 91% recall), followed closely by random forest, QDA and LDA (81% to 83%). Every one of these differences is smaller than the split-to-split standard deviation (5 to 9 points), so the four models cannot be told apart reliably, and all of them are only a few points above the always-Parkinson’s baseline (78.2%). QDA has the highest and most stable recall (about 98%) but the lowest precision of the four, in line with its tendency to call most recordings Parkinson’s. Plain logistic regression and Naive Bayes are below the baseline on average, and the single decision tree is about level with it.

Question 4

Interpretability v Predictive performance

Interpretable models like Logistic Regression, LDA, and QDA allow clinicians to understand how individual features contribute to the prediction. This transparency is critical in medical settings, where treatment decisions require justification.

On the other hand, complex models like Random Forest, Boosting, or SVMs often achieve higher raw predictive power, especially in capturing nonlinear relationships and interactions. However, they operate as “black boxes” and offer limited insight into why a prediction was made.

In this study, PCA + Logistic Regression had the highest mean accuracy (about 84%) in the repeated-split check while remaining a simple, transparent model, and its principal components can be traced back to groups of correlated voice features. That is consistent with the idea that a well-regularised interpretable model can match black-box models on a small structured dataset, but the margin over random forest, QDA and LDA is well within the noise, so the evidence here supports “no worse than” rather than “better than”.

Discussion

This project evaluated parametric and non-parametric classifiers for detecting Parkinson’s disease from voice measurements. The most important lesson concerned evaluation. Splitting the data by recording put the same people in both the training and test sets, and under that split several models scored 85% to 95% accuracy. When whole people are held out, which is the realistic setting for a screening tool, performance drops sharply: on the single held-out set of 9 people most models either matched the trivial “everyone has Parkinson’s” baseline (78.2%) or fell below it, and the model that looked best there (Naive Bayes) was among the worst on average across 30 repeated splits.

Across those repeated subject-level splits, PCA + logistic regression had the best mean accuracy (about 84%, with about 90% precision and 91% recall), with random forest, QDA and LDA within one standard deviation of it and only a few points above the baseline. Plain logistic regression and Naive Bayes did worst, and regularisation or PCA (which handle the strong multicollinearity) clearly helped logistic regression. Cross-validation through caret did not change the results of models without tuning parameters.

In short, the voice features do carry some signal that generalises to unseen people, but with only 32 people in the data the top models cannot be reliably ranked, and recall figures near 100% can simply mean that a model calls everyone Parkinson’s. PCA + logistic regression is a reasonable, interpretable choice, not a demonstrated winner.

Limitations

The dataset is small: 195 recordings from only 32 people (8 healthy, 24 with Parkinson’s). Evaluating on held-out people leaves few people in any test set (9 here, only 2 of them healthy), so results vary a lot from split to split (standard deviations of 5 to 10 accuracy points in the robustness check).

Some models, such as LDA, required feature removal due to zero variance, possibly reducing their full potential.

Class imbalance: about 78% of the recordings are from people with Parkinson’s, so 78% accuracy and 100% recall are achievable by labelling everyone as Parkinson’s. Models must be compared with that baseline, and accuracy and recall in particular can be misleading on their own.

The first version of this analysis split the data by recording rather than by person, which leaked information between training and test sets and inflated every score; the results here use a subject-level split instead.

The folds used for hyperparameter tuning and cross-validation inside the training set are still random over recordings, so tuned values (for example the SVM cost and gamma) may be optimistic.

What We Would Do Differently

Use grouped cross-validation (for example leave-one-person-out) for both tuning and evaluation, so that every person is used as a test case and the estimates are far less noisy than a single 9-person test set

Perform best subsets regression to better understand the importance of features and remove redundant features

Use more advanced techniques such as SMOTE (Synthetic Minority Over-sampling) to better balance the classes.

Include hyperparameter tuning via grid search for all models to ensure optimal performance under consistent validation.

Future Directions

Use a larger and more diverse dataset, possibly including multiple modalities (e.g., imaging, wearable sensor data).

Explore ensemble methods combining parametric and non-parametric models for robustness.

Incorporate time-series models or longitudinal patient data to track disease progression.

Build an interactive diagnostic tool or dashboard to demonstrate real-world utility to clinicians.

Project

Dinesh Kumar Gummadavelli

2025-04-22