{"metadata":{"kernelspec":{"name":"ir","display_name":"R","language":"R"},"language_info":{"name":"R","codemirror_mode":"r","pygments_lexer":"r","mimetype":"text/x-r-source","file_extension":".r","version":"4.4.0"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":84896,"databundleVersionId":10305135,"sourceType":"competition"}],"dockerImageVersionId":30749,"isInternetEnabled":true,"language":"r","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"## Poor RMSE with linear regression\nMuch below suggests this is a poor model such as RMSE = 863.853 and the R-squared 0.0026. \n\nStill, a fun experiment.\n\nWe may not be able to have R notebooks scored.","metadata":{}},{"cell_type":"code","source":"suppressPackageStartupMessages({\nlibrary(tidyverse, quietly = TRUE)\nlibrary(janitor)\nlibrary(MASS)      # For stepAIC function\nlibrary(car)       # For VIF calculation\nlibrary(lmtest)    # For Breusch-Pagan test\nlibrary(caret)     # For preprocessing\nlibrary(nortest)\n})\n\nselect <- dplyr::select","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","trusted":true,"execution":{"iopub.status.busy":"2024-12-10T05:24:54.863134Z","iopub.execute_input":"2024-12-10T05:24:54.865756Z","iopub.status.idle":"2024-12-10T05:24:58.839663Z","shell.execute_reply":"2024-12-10T05:24:58.837766Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"df_train=read.csv('/kaggle/input/playground-series-s4e12/train.csv') %>% clean_names()\ndf_test=read.csv('/kaggle/input/playground-series-s4e12/test.csv') %>% clean_names()\nsample_sub=read.csv('/kaggle/input/playground-series-s4e12/sample_submission.csv')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-10T05:24:58.843728Z","iopub.execute_input":"2024-12-10T05:24:58.878854Z","iopub.status.idle":"2024-12-10T05:25:31.658157Z","shell.execute_reply":"2024-12-10T05:25:31.656335Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"colSums(is.na(df_train))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-10T05:25:31.661401Z","iopub.execute_input":"2024-12-10T05:25:31.662919Z","iopub.status.idle":"2024-12-10T05:25:32.586545Z","shell.execute_reply":"2024-12-10T05:25:32.584836Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## NA annual_income, make the median (about $25000)","metadata":{}},{"cell_type":"code","source":"df_train$annual_income[is.na(df_train$annual_income)] <- 25000","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-10T05:25:32.589137Z","iopub.execute_input":"2024-12-10T05:25:32.590578Z","iopub.status.idle":"2024-12-10T05:25:32.608519Z","shell.execute_reply":"2024-12-10T05:25:32.606641Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Median per income bucket for missing values ","metadata":{}},{"cell_type":"code","source":"df_train <- df_train %>%\n  mutate(income_bucket = cut(annual_income, breaks = seq(0, max(annual_income, na.rm = TRUE) + 25000, by = 25000), include.lowest = TRUE))\n\n# Impute NA values with median for each numeric column per income bucket\ndf_train <- df_train %>%\n  group_by(income_bucket) %>%\n  mutate(across(where(is.numeric), ~ifelse(is.na(.), median(., na.rm = TRUE), .))) %>%\n  ungroup() %>%\n  select(-income_bucket) \n\ndf_test <- df_test %>%\n  mutate(income_bucket = cut(annual_income, breaks = seq(0, max(annual_income, na.rm = TRUE) + 25000, by = 25000), include.lowest = TRUE))\n\n# Impute NA values with median for each numeric column per income bucket\ndf_test <- df_test %>%\n  group_by(income_bucket) %>%\n  mutate(across(where(is.numeric), ~ifelse(is.na(.), median(., na.rm = TRUE), .))) %>%\n  ungroup() %>%\n  select(-income_bucket) ","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-10T05:25:32.611106Z","iopub.execute_input":"2024-12-10T05:25:32.612558Z","iopub.status.idle":"2024-12-10T05:25:35.103118Z","shell.execute_reply":"2024-12-10T05:25:35.101396Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Preprocess the data: Normalize predictors only\npreprocessor1 <- preProcess(df_train %>% select_if(is.numeric) %>% select(-premium_amount), method = c(\"center\", \"scale\"))\npreprocessor2 <- preProcess(df_test %>% select_if(is.numeric), method = c(\"center\", \"scale\"))\n\ndf_train_normalized <- df_train\ndf_train_normalized[df_train_normalized %>% select_if(is.numeric) %>% select(-premium_amount) %>% colnames()] <- predict(preprocessor1, df_train %>% select_if(is.numeric) %>% select(-premium_amount))\n\ndf_test_normalized <- df_test\ndf_test_normalized[df_test_normalized %>% select_if(is.numeric) %>% colnames()] <- predict(preprocessor2, df_test %>% select_if(is.numeric))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-10T05:25:35.105896Z","iopub.execute_input":"2024-12-10T05:25:35.107351Z","iopub.status.idle":"2024-12-10T05:25:46.273755Z","shell.execute_reply":"2024-12-10T05:25:46.271836Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"numeric_columns <- df_train %>% select_if(is.numeric) %>% names()\ncorrelations <- cor(df_train %>% select_if(is.numeric), df_train$premium_amount)\n\ncorrelation_df <- data.frame(Variable = numeric_columns, Correlation = correlations)\n\ncorrelation_df <- correlation_df[order(correlation_df$Correlation), ]\nprint(correlation_df)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-10T05:25:46.276697Z","iopub.execute_input":"2024-12-10T05:25:46.27821Z","iopub.status.idle":"2024-12-10T05:25:46.929296Z","shell.execute_reply":"2024-12-10T05:25:46.927345Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"summary(df_train)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-10T05:25:46.932171Z","iopub.execute_input":"2024-12-10T05:25:46.93372Z","iopub.status.idle":"2024-12-10T05:25:47.616136Z","shell.execute_reply":"2024-12-10T05:25:47.614369Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Forward Selection Linear Regression - Model","metadata":{}},{"cell_type":"code","source":"# Define the initial model (only intercept)\ninitial_model <- lm(premium_amount ~ 1, data = df_train_normalized)\n\n# Define the full model (all normalized predictors)\nfull_model <- lm(premium_amount ~ ., data = df_train_normalized %>% select(-id) %>% select_if(is.numeric))\n\n# Perform forward selection using stepAIC\nforward_model <- stepAIC(initial_model, \n                         scope = list(lower = initial_model, upper = full_model), \n                         direction = \"forward\", \n                         trace = FALSE)\n\n# Summary of the final selected model\ncat(\"Summary of the final selected model:\\n\")\nsummary(forward_model)\n\n# Check for multicollinearity using Variance Inflation Factor (VIF)\nvif_values <- vif(forward_model)\ncat(\"Variance Inflation Factor (VIF):\\n\")\nprint(vif_values)\n\n# Identify predictors with high VIF\nif (any(vif_values > 5)) {\n  cat(\"Warning: High multicollinearity detected for the following predictors:\\n\")\n  print(names(vif_values[vif_values > 5]))\n}\n\n# Performance metrics for the final model\n# Get predicted values\npredictions <- predict(forward_model, newdata = df_train_normalized)\n\n# Calculate RMSE (Root Mean Squared Error)\nrmse <- sqrt(mean((df_train$premium_amount - predictions)^2))  # Use original target for residuals\ncat(\"Root Mean Squared Error (RMSE):\", rmse, \"\\n\")\n\n# Calculate R-squared\nr_squared <- summary(forward_model)$r.squared\ncat(\"R-squared:\", r_squared, \"\\n\")\n\n# Diagnostic plots\npar(mfrow = c(2, 2)) # Set plotting area to show multiple diagnostic plots\nplot(forward_model)\npar(mfrow = c(1, 1)) # Reset plotting area\n\n# Residual Analysis\nresiduals <- resid(forward_model)\ncat(\"Residuals Summary:\\n\")\nprint(summary(residuals))\n\n# Check for heteroscedasticity using the Breusch-Pagan test\nbp_test <- bptest(forward_model)\ncat(\"Breusch-Pagan Test for Heteroscedasticity:\\n\")\nprint(bp_test)\n\n# Check for normality of residuals \nad_test <- ad.test(residuals)\nprint(ad_test)\n\n# Cross-validation to validate the model\nset.seed(123)\ntrain_control <- trainControl(method = \"cv\", number = 10)\ncross_val_model <- train(premium_amount ~ ., \n                         data = df_train_normalized %>% select_if(is.numeric), \n                         method = \"lm\", \n                         trControl = train_control)\n\ncat(\"Cross-Validation Results:\\n\")\nprint(cross_val_model)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-10T05:25:47.618891Z","iopub.execute_input":"2024-12-10T05:25:47.620424Z","iopub.status.idle":"2024-12-10T05:29:57.236085Z","shell.execute_reply":"2024-12-10T05:29:57.217729Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"summary(forward_model)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Check for multicollinearity using Variance Inflation Factor (VIF)\nvif_values <- vif(forward_model)\ncat(\"Variance Inflation Factor (VIF):\\n\")\nprint(vif_values)\n\n# Identify predictors with high VIF\nif (any(vif_values > 5)) {\n  cat(\"Warning: High multicollinearity detected for the following predictors:\\n\")\n  print(names(vif_values[vif_values > 5]))\n}\n\n# Performance metrics for the final model\n# Get predicted values\npredictions <- predict(forward_model, newdata = df_train)\n\n# Calculate RMSE (Root Mean Squared Error)\nrmse <- sqrt(mean((df_train$premium_amount - predictions)^2))\ncat(\"Root Mean Squared Error (RMSE):\", rmse, \"\\n\")\n\n# Calculate R-squared\nr_squared <- summary(forward_model)$r.squared\ncat(\"R-squared:\", r_squared, \"\\n\")\n\n# Diagnostic plots\npar(mfrow = c(2, 2)) # Set plotting area to show multiple diagnostic plots\nplot(forward_model)\npar(mfrow = c(1, 1)) # Reset plotting area\n\n# Residual Analysis\nresiduals <- resid(forward_model)\ncat(\"Residuals Summary:\\n\")\nprint(summary(residuals))\n\n# Check for heteroscedasticity using the Breusch-Pagan test\nbp_test <- bptest(forward_model)\ncat(\"Breusch-Pagan Test for Heteroscedasticity:\\n\")\nprint(bp_test)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"predictions <- predict(forward_model, newdata = df_test)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Make predictions using the forward_model\npredictions_test <- predict(forward_model, newdata = df_test_normalized)\n\n# Update the 'Premium.Amount' column in the sample submission file\nsample_sub$`Premium Amount` <- predictions_test\n\n# Save the updated submission file\nwrite.csv(sample_sub, \"submission.csv\", row.names = FALSE)\n\ncat(\"Submission file 'submission.csv' created successfully.\\n\")","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}