# Split data for validation
n_train <- round(0.75 * n)
train_idx <- sample(n, n_train)
test_idx <- setdiff(1:n, train_idx)
# Refit on training data
cvfit_train <- cv.missoNet(
X = sim$X[train_idx, ],
Y = sim$Z[train_idx, ],
kfold = 5,
lambda.beta = lambda.beta.refined,
lambda.theta = lambda.theta.refined,
verbose = 0
)
# Predictions
Y_pred <- predict(cvfit_train, newx = sim$X[test_idx, ])
Y_test <- sim$Y[test_idx, ] # True complete values
# Calculate performance metrics
mse_per_cpg <- colMeans((Y_pred - Y_test)^2)
cor_per_cpg <- sapply(1:q, function(j) cor(Y_pred[,j], Y_test[,j]))
# Visualization
par(mfrow = c(1, 2))
# MSE vs missing rate
plot(miss_by_cpg, mse_per_cpg, pch = 19, col = "darkblue",
xlab = "Missing Rate", ylab = "Prediction MSE",
main = "Prediction Error vs Missing Rate")
lines(lowess(miss_by_cpg, mse_per_cpg), col = "red", lwd = 2)
# Correlation vs missing rate
plot(miss_by_cpg, cor_per_cpg, pch = 19, col = "darkgreen",
xlab = "Missing Rate", ylab = "Prediction Correlation",
main = "Prediction Accuracy vs Missing Rate")
lines(lowess(miss_by_cpg, cor_per_cpg), col = "red", lwd = 2)
{
cat("\nPrediction Performance:\n")
cat(" Mean MSE:", round(mean(mse_per_cpg), 4), "\n")
cat(" Mean correlation:", round(mean(cor_per_cpg), 3), "\n")
cat(" Worst CpG MSE:", round(max(mse_per_cpg), 4), "\n")
cat(" Best CpG correlation:", round(max(cor_per_cpg), 3), "\n")
}
#>
#> Prediction Performance:
#> Mean MSE: 0.9618
#> Mean correlation: 0.842
#> Worst CpG MSE: 1.4674
#> Best CpG correlation: 0.965
