{"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":87793,"databundleVersionId":11228175,"sourceType":"competition"}],"dockerImageVersionId":30749,"isInternetEnabled":false,"language":"r","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"\n# <div align=\"center\"> Stanford RNA 3D Folding \n    \n\n<div align=\"center\"> <img src=\"https://www.kaggle.com/competitions/87793/images/header\" width=\"300\"></div><br>\n\n**From the competition material**: For each sequence in the test set, you can predict <ins>five structures</ins>. Your notebook should look for a file test_sequences.csv and output submission.csv. This file should contain x, y, z coordinates of the C1' atom in each residue across your predicted structures 1 to 5:\n\n\n**ID**,**resname**,**resid**,**x_1**,**y_1**,**z_1**,... **x_5**,**y_5**,**z_5** <br>\nR1107_1,G,1,-7.561,9.392,9.361,... -7.301,9.023,8.932<br>\nR1107_2,G,1,-8.02,11.014,14.606,... -7.953,10.02,12.127<br>\netc.<br>\n\n\n**Evaluation**: Submissions are scored using <ins>TM-score</ins> (\"template modeling\" score), which goes from 0.0 to 1.0 (higher is better):\n","metadata":{}},{"cell_type":"markdown","source":"# 📚 Libraries","metadata":{}},{"cell_type":"code","source":"library(tidyverse) # metapackage of all tidyverse packages\nlibrary(h2o)\nh2o.init()","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","trusted":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2025-03-06T00:41:09.834995Z","iopub.execute_input":"2025-03-06T00:41:09.837511Z","iopub.status.idle":"2025-03-06T00:41:20.41685Z","shell.execute_reply":"2025-03-06T00:41:20.412683Z"},"_kg_hide-input":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# 🔍 Exploratory Data Analysis (EDA)","metadata":{}},{"cell_type":"code","source":"# Reading in the competition material\ntrain_sequences <- read.csv('/kaggle/input/stanford-rna-3d-folding/train_sequences.csv') \ntrain_labels <- read.csv('/kaggle/input/stanford-rna-3d-folding/train_labels.csv') \nvalidation_sequences <- read.csv('/kaggle/input/stanford-rna-3d-folding/validation_sequences.csv') \nvalidation_labels <- read.csv('/kaggle/input/stanford-rna-3d-folding/validation_labels.csv') ","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-06T00:41:20.423004Z","iopub.execute_input":"2025-03-06T00:41:20.52636Z","iopub.status.idle":"2025-03-06T00:41:21.676056Z","shell.execute_reply":"2025-03-06T00:41:21.674097Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### train_sequences.csv","metadata":{}},{"cell_type":"code","source":"print(paste('train_sequences has', dim(train_sequences)[1], 'rows and', dim(train_sequences)[2], 'columns'))\nprint(paste('train_sequences has', sum(is.na(train_sequences)), 'NAs.'))\n\n# Quick look at the data\ntrain_sequences %>% head(3)\n\n# There are no NAs in this dataset","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-06T00:41:21.679566Z","iopub.execute_input":"2025-03-06T00:41:21.681307Z","iopub.status.idle":"2025-03-06T00:41:21.723894Z","shell.execute_reply":"2025-03-06T00:41:21.721901Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Let's focus on the sequence column...specifically the length of the sequences.","metadata":{}},{"cell_type":"code","source":"options(repr.plot.width=15, repr.plot.height=8)\ntrain_sequences %>% select(sequence) %>% mutate(length = str_count(sequence)) %>% \n    ggplot(aes(x=length)) + geom_boxplot() + \n    xlab(\"Sequence length\") + ylab(\"Sequence\") + ggtitle(\"Boxplot of sequence length values (with outliers)\") + theme_minimal()\n\n# That's kind of messy...let's filter out some of the outliers","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-06T00:41:21.727269Z","iopub.execute_input":"2025-03-06T00:41:21.728799Z","iopub.status.idle":"2025-03-06T00:41:22.579494Z","shell.execute_reply":"2025-03-06T00:41:22.577448Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"options(repr.plot.width=15, repr.plot.height=8)\ntrain_sequences %>% select(sequence) %>% mutate(length = str_count(sequence)) %>% \n    ggplot(aes(x=length)) + geom_boxplot(outliers = FALSE) + \n    xlab(\"Sequence length\") + ylab(\"Sequence\") + ggtitle(\"Boxplot of sequence length values (without outliers)\") + theme_minimal()\n\n# By using the outliers = FALSE argument within geom_boxplot(), we can create something a little easier to make sense of","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-06T00:41:22.582434Z","iopub.execute_input":"2025-03-06T00:41:22.583911Z","iopub.status.idle":"2025-03-06T00:41:22.938078Z","shell.execute_reply":"2025-03-06T00:41:22.936086Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Let's look at the temporal_cutoff column and see how the length of published sequences has evolved over time","metadata":{}},{"cell_type":"code","source":"options(repr.plot.width=15, repr.plot.height=8)\ntrain_sequences %>% mutate(cutoff_year = lubridate::year(ymd(temporal_cutoff)), sequence_length = str_count(sequence)) %>% \n    group_by(cutoff_year) %>% summarize(avg_sequence_length = mean(sequence_length)) %>% \n    ggplot(aes(x=cutoff_year, y=avg_sequence_length)) + geom_smooth(method = \"lm\", formula = y ~ poly(x, 3), se = FALSE) + geom_line() +\n    xlab(\"Temporal Cutoff Year\") + ylab(\"Average Sequence Length discovered\") + ggtitle(\"Average length of published sequence over time\") + theme_minimal()\n\n# The black line is actual sequence length and the blue line is smoothed.  Both show a gradual increase in published sequence length \n# since 2000, with significant published sequence lengths since 2014 timeframe.","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-06T00:41:22.940895Z","iopub.execute_input":"2025-03-06T00:41:22.94228Z","iopub.status.idle":"2025-03-06T00:41:24.683318Z","shell.execute_reply":"2025-03-06T00:41:24.681185Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Focussing on the target_id...from the competition material:\n>\"In train_sequences.csv, this is formatted as pdb_id_chain_id, where pdb_id is the id of the entry in the Protein Data Bank and chain_id is the chain id of the monomer in the pdb file","metadata":{}},{"cell_type":"markdown","source":"#### Let's first look at how many unique pdb_id's and chain_id's there are in the train_sequences","metadata":{"execution":{"iopub.status.busy":"2025-03-02T22:06:18.770555Z","iopub.execute_input":"2025-03-02T22:06:18.772152Z","iopub.status.idle":"2025-03-02T22:06:18.78293Z","shell.execute_reply":"2025-03-02T22:06:18.781167Z"}}},{"cell_type":"code","source":"pdb_id_count <- train_sequences %>% mutate(pdb_id = str_split(target_id, \"_\", simplify = TRUE)[ , 1], chain_id = str_split(target_id, \"_\", simplify = TRUE)[ , 2]) %>% \n    select(pdb_id) %>% unique() %>% count() %>% pull()\n\nchain_id_count <- train_sequences %>% mutate(pdb_id = str_split(target_id, \"_\", simplify = TRUE)[ , 1], chain_id = str_split(target_id, \"_\", simplify = TRUE)[ , 2]) %>% \n    select(chain_id) %>% unique() %>% count() %>% pull()\n\nprint(paste(\"There are \",pdb_id_count,\"unique pdb_id's in the dataset.\"))\nprint(paste(\"There are \",chain_id_count,\"unique pdb_id's in the dataset.\"))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-06T00:41:24.685998Z","iopub.execute_input":"2025-03-06T00:41:24.687368Z","iopub.status.idle":"2025-03-06T00:41:24.731378Z","shell.execute_reply":"2025-03-06T00:41:24.72937Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Now let's look at which pdb_id's are repeated in the train_sequences","metadata":{}},{"cell_type":"code","source":"train_sequences %>% mutate(pdb_id = str_split(target_id, \"_\", simplify = TRUE)[ , 1], chain_id = str_split(target_id, \"_\", simplify = TRUE)[ , 2]) %>% \n    group_by(pdb_id) %>% summarize(total = n()) %>% filter(total > 1) %>% arrange(desc(total)) %>% head(10)\n\n# In total, there are 60 pdb_id's that are repeated in train_sequences, with the top 10 most repeated shown below","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-06T00:41:24.73412Z","iopub.execute_input":"2025-03-06T00:41:24.735556Z","iopub.status.idle":"2025-03-06T00:41:24.784719Z","shell.execute_reply":"2025-03-06T00:41:24.782928Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### train_labels.csv","metadata":{}},{"cell_type":"code","source":"print(paste('train_labels has', dim(train_labels)[1], 'rows and', dim(train_labels)[2], 'columns'))\nprint(paste('train_labels has', sum(is.na(train_labels)), 'NAs.'))\n# Quick look at the data\ntrain_labels %>% head(3)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-06T00:41:24.78741Z","iopub.execute_input":"2025-03-06T00:41:24.788817Z","iopub.status.idle":"2025-03-06T00:41:24.840576Z","shell.execute_reply":"2025-03-06T00:41:24.838681Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# There are a decent number of NAs in train_labels...here are a few:\ntrain_labels[rowSums(is.na(train_labels)) > 0,] %>% head(3)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-06T00:41:24.843389Z","iopub.execute_input":"2025-03-06T00:41:24.844907Z","iopub.status.idle":"2025-03-06T00:41:24.887384Z","shell.execute_reply":"2025-03-06T00:41:24.885335Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Let's create a boxplot of the coordinate values for x, y, & z ","metadata":{}},{"cell_type":"code","source":"options(repr.plot.width=15, repr.plot.height=8)\ntrain_labels %>% select(x_1, y_1, z_1) %>% \n    pivot_longer(everything(), names_to = \"coordinate\", values_to = \"value\") %>% \n    ggplot(aes(x=value, fill = coordinate)) + geom_boxplot() + facet_grid(cols=vars(coordinate)) + \n    coord_flip() + scale_fill_brewer(palette = \"YlGnBu\") + \n    xlab(\"Coordinate Values\") + ylab(\"Coordinates\") + ggtitle(\"Boxplot of coordinate values in train_labels\") + theme_minimal()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-06T00:41:24.890474Z","iopub.execute_input":"2025-03-06T00:41:24.89199Z","iopub.status.idle":"2025-03-06T00:41:27.152511Z","shell.execute_reply":"2025-03-06T00:41:27.150654Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Resname occurrences in train_labels","metadata":{}},{"cell_type":"code","source":"options(repr.plot.width=15, repr.plot.height=8)\ntrain_labels %>% group_by(resname) %>% summarize(total = n()) %>% \n    ggplot(aes(x=reorder(resname, total), y=total)) + geom_col(aes(fill=resname)) + geom_text(aes(label = total), vjust = -0.5) +  scale_fill_brewer(palette = \"Blues\") + \n    xlab(\"Resname\") + ylab(\"Count\") + ggtitle(\"Barplot of resname occurrences in train_labels\") + guides(fill=\"none\") + theme_minimal()\n\n# There are 2 x resnames with 'x' and four that are missing...the rest come from U, A, C, & G.","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-06T00:41:27.15511Z","iopub.execute_input":"2025-03-06T00:41:27.156583Z","iopub.status.idle":"2025-03-06T00:41:27.476453Z","shell.execute_reply":"2025-03-06T00:41:27.474497Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### validation_sequences.csv","metadata":{}},{"cell_type":"code","source":"print(paste('validation_sequences has', dim(validation_sequences)[1], 'rows and', dim(validation_sequences)[2], 'columns'))\nprint(paste('validation_sequences has', sum(is.na(validation_sequences)), 'NAs.'))\n\n# Quick look at the data\nvalidation_sequences %>% head(3)\n\n# This dataset has the same number of columns as train_sequences, and it also has 0 NA's.","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-06T00:41:27.479296Z","iopub.execute_input":"2025-03-06T00:41:27.480785Z","iopub.status.idle":"2025-03-06T00:41:27.512293Z","shell.execute_reply":"2025-03-06T00:41:27.51061Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Let's recreate the boxplot of the sequence lengths we used in train_sequences, using outliers = FALSE","metadata":{}},{"cell_type":"code","source":"options(repr.plot.width=15, repr.plot.height=8)\nvalidation_sequences %>% select(sequence) %>% mutate(length = str_count(sequence)) %>% \n    ggplot(aes(x=length)) + geom_boxplot(outliers = FALSE) + \n    xlab(\"Sequence length\") + ylab(\"Sequence\") + ggtitle(\"Boxplot of sequence length values (without outliers)\") + theme_minimal()\n\n# A quick comparison of the two boxplots shows the medians are drastically different for the sequence lengths between validation_sequences and train_sequences","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-06T00:41:27.514888Z","iopub.execute_input":"2025-03-06T00:41:27.516267Z","iopub.status.idle":"2025-03-06T00:41:27.77864Z","shell.execute_reply":"2025-03-06T00:41:27.776651Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train_sequence_med <- train_sequences %>% select(sequence) %>% mutate(length = str_count(sequence)) %>% select(length) %>% pull() %>% median()\nvalid_sequence_med <- validation_sequences %>% select(sequence) %>% mutate(length = str_count(sequence)) %>% select(length) %>% pull() %>% median()\n\nprint(paste('The train_sequences median sequence length is', train_sequence_med))\nprint(paste('The validation_sequences median sequence length is', valid_sequence_med))\n\n# That's quite a difference!","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-06T00:41:27.781382Z","iopub.execute_input":"2025-03-06T00:41:27.782822Z","iopub.status.idle":"2025-03-06T00:41:27.836506Z","shell.execute_reply":"2025-03-06T00:41:27.83479Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### validation_labels.csv","metadata":{}},{"cell_type":"code","source":"print(paste('validation_labels has', dim(validation_labels)[1], 'rows and', dim(validation_labels)[2], 'columns'))\nprint(paste('validation_labels has', sum(is.na(validation_labels)), 'NAs.'))\n# Quick look at the data\nvalidation_labels %>% head(3)\n\n# Note validation_labels has 123 columns; train_labels has 6","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-06T00:41:27.83917Z","iopub.execute_input":"2025-03-06T00:41:27.840539Z","iopub.status.idle":"2025-03-06T00:41:27.917813Z","shell.execute_reply":"2025-03-06T00:41:27.916083Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### The number of NAs is misleading; from a discussion post [here](https://www.kaggle.com/competitions/stanford-rna-3d-folding/discussion/565746), the -1e+18 should be interpreted as NA/NaN values.\n\n### We will recreate the resname boxplot from train_labels from earlier:","metadata":{}},{"cell_type":"code","source":"options(repr.plot.width=15, repr.plot.height=8)\nvalidation_labels %>% group_by(resname) %>% summarize(total = n()) %>% \n    ggplot(aes(x=reorder(resname, total), y=total)) + geom_col(aes(fill=resname)) + geom_text(aes(label = total), vjust = -0.5) +  scale_fill_brewer(palette = \"Blues\") + \n    xlab(\"Resname\") + ylab(\"Count\") + ggtitle(\"Barplot of resname occurrences in valid_labels\") + guides(fill=\"none\") + theme_minimal()\n\n# Note there are no missing or 'X' resnames.  Also of note is in this datset, 'U' occurs more frequently than 'A'.  'C' and 'G' maintain their ordering.","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-06T00:41:27.920524Z","iopub.execute_input":"2025-03-06T00:41:27.921904Z","iopub.status.idle":"2025-03-06T00:41:28.535109Z","shell.execute_reply":"2025-03-06T00:41:28.532917Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Submission format","metadata":{}},{"cell_type":"code","source":"# This is our submission format\nsample_submission <- read.csv('/kaggle/input/stanford-rna-3d-folding/sample_submission.csv')\nprint(paste('There are ',dim(sample_submission)[1], 'rows and',dim(sample_submission)[2], 'columns.'))\nsample_submission %>% head(6)\n\n# Note: It is required to submit 5 sets of coordinate plots (x,y,z), per ID, resname & resid from the target_id, which is 18 columns.\n# For rows in the submission, we need one row for each resname/resid for each target\n# From the sample submission, we need 2515 rows","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-06T00:41:28.537772Z","iopub.execute_input":"2025-03-06T00:41:28.539208Z","iopub.status.idle":"2025-03-06T00:41:28.607005Z","shell.execute_reply":"2025-03-06T00:41:28.605257Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### Generating the output file from test_sequence","metadata":{}},{"cell_type":"code","source":"test_sequences <- read.csv('/kaggle/input/stanford-rna-3d-folding/test_sequences.csv')\nprint(paste('There are ',dim(test_sequences)[1], 'rows and',dim(test_sequences)[2], 'columns.'))\n\n# Reading in the test sequences, we can see we have 12 x rows and the same 5 x columns found in train & validation sequences.","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-06T00:41:28.609595Z","iopub.execute_input":"2025-03-06T00:41:28.610978Z","iopub.status.idle":"2025-03-06T00:41:28.636276Z","shell.execute_reply":"2025-03-06T00:41:28.634447Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# By summing each sequence length (equivalent to one resname/resid), for the entire test_sequences, we can arrive at our number of 2515.\ntest_sequences %>% mutate(length = nchar(sequence)) %>% select(length) %>% sum(.)\n\n# Now, we need to create the df to house our submission","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-06T00:41:28.6391Z","iopub.execute_input":"2025-03-06T00:41:28.640542Z","iopub.status.idle":"2025-03-06T00:41:28.66147Z","shell.execute_reply":"2025-03-06T00:41:28.659623Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# This function will accept a test_sequence ID & sequence and return df in format of submission (one resid/resname per row in sequence)\nparse_target <- function(tmp_ID, tmp_sequence){\n    seq_length <- nchar(tmp_sequence)\n    tmp_df=data.frame(matrix(ncol = 3, nrow = seq_length)) \n    colnames(tmp_df) <- c('ID','resname','resid')\n    tmp_df$resname <- unlist(strsplit(tmp_sequence,split=''))\n    tmp_df$ID <- tmp_ID\n    tmp_df$resid <- seq(1:seq_length)\n\n    return(tmp_df)\n}\n# Create the test_clean df in the format of the submission\ntest_id_seq <- test_sequences %>% select(target_id, sequence)\ntest_clean <- data.frame(matrix(ncol = 3, nrow = 0)) \ncolnames(test_clean) <- c('ID','resname','resid')\n\n# For each target_id / sequence, apply the function and rbind to previous results\nfor(i in 1:nrow(test_id_seq)){\n    tmp_df <- as.data.frame(mapply(parse_target,test_id_seq$target_id[i],test_id_seq$sequence[i], SIMPLIFY=FALSE))\n    colnames(tmp_df) <- c('ID','resname','resid')\n    test_clean <- rbind(test_clean,tmp_df)\n}\n\nprint(paste('There are',nrow(test_clean),'rows in the test_clean df.'))\nprint(head(test_clean))\n\n# That works...now let's get back to the training data, add some features, & train a model.","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-06T00:41:28.664095Z","iopub.execute_input":"2025-03-06T00:41:28.665463Z","iopub.status.idle":"2025-03-06T00:41:28.722038Z","shell.execute_reply":"2025-03-06T00:41:28.720261Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# 🎯 Model Training","metadata":{}},{"cell_type":"code","source":"# Function to create some features based off resname combinations & sequence length\nfeat_eng <- function(df){\n    df <- df %>% mutate(seq_length = nchar(sequence), \n               A_cnt = str_count(sequence,\"A\"), \n               C_cnt = str_count(sequence,\"C\"), \n               U_cnt = str_count(sequence,\"U\"), \n               G_cnt = str_count(sequence,\"G\"),                  \n               AC_cnt = str_count(sequence,\"AC\"), \n               AU_cnt = str_count(sequence,\"AU\"), \n               AG_cnt = str_count(sequence,\"AG\"),\n               CA_cnt = str_count(sequence,\"CA\"), \n               CU_cnt = str_count(sequence,\"CU\"), \n               CG_cnt = str_count(sequence,\"CG\"),\n               UA_cnt = str_count(sequence,\"UA\"), \n               UC_cnt = str_count(sequence,\"UC\"), \n               UG_cnt = str_count(sequence,\"UG\"),\n               GA_cnt = str_count(sequence,\"GA\"), \n               GC_cnt = str_count(sequence,\"GC\"), \n               GU_cnt = str_count(sequence,\"GU\"), \n               AA_cnt = str_count(sequence,\"AA\"), \n               CC_cnt = str_count(sequence,\"CC\"), \n               UU_cnt = str_count(sequence,\"UU\"), \n               GG_cnt = str_count(sequence,\"GG\"),\n               begin_seq = substr(sequence,1,1),\n               end_seq = substr(sequence,nchar(sequence),nchar(sequence))) %>% select(-c(sequence, temporal_cutoff, description, all_sequences))\n    return(df)\n}\ntrain_sequences <- feat_eng(train_sequences)\n\n# Create the target_id from ID to facilitate joining with the meta_train data we just created\ntrain_labels <- train_labels %>% group_by(ID) %>% mutate(target_id = paste(unlist(strsplit(ID,'_'))[1],unlist(strsplit(ID,'_'))[2],sep='_'))\n\n# Create the final df to train\ntrain_data_clean <- merge(train_labels, train_sequences, by=\"target_id\", all.x = TRUE)\n\n# For now, we will impute the missing x_1,y_1,z_1 values with the group average for that target_id, resname\ntrain_data_clean <- train_data_clean %>% group_by(target_id, resname) %>% mutate(x_1 = ifelse(is.na(x_1),mean(x_1, na.rm=TRUE),x_1), y_1 = ifelse(is.na(y_1),mean(y_1, na.rm=TRUE),y_1), z_1 = ifelse(is.na(z_1),mean(z_1, na.rm=TRUE),z_1)) %>% ungroup() %>% select(-target_id)\n\n# The above imputation for train_data_clean doesn't take care of all NAs...so for now we will remove them(~ 1979 rows)\ntrain_data_clean <- train_data_clean %>% na.exclude() \nprint(paste('There are ',sum(is.na(train_data_clean)),\"NA's in the dataset!\"))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-06T00:41:28.724751Z","iopub.execute_input":"2025-03-06T00:41:28.726138Z","iopub.status.idle":"2025-03-06T00:41:34.461925Z","shell.execute_reply":"2025-03-06T00:41:34.460054Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"features <- train_data_clean %>% select(-c(x_1,y_1,z_1, ID)) %>% colnames()  ","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-06T00:41:34.464673Z","iopub.execute_input":"2025-03-06T00:41:34.466186Z","iopub.status.idle":"2025-03-06T00:41:34.480787Z","shell.execute_reply":"2025-03-06T00:41:34.47902Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### For now, our model approach will consist of 3 separate models for x,y,z values","metadata":{}},{"cell_type":"code","source":"# Train model for X\n\nSEED=42\nMAX_DEPTH=15\nNTREES=1100\nLRATE=0.15\n\ntrain_hframe = train_data_clean %>%  select(-c(ID,y_1,z_1)) %>% as.h2o() # will need to separate out x_1,y_1,z_1 & ID for each iteration\ntrain_hframe$resname <- as.factor(train_hframe$resname)\ntrain_hframe$begin_seq <- as.factor(train_hframe$begin_seq)\ntrain_hframe$end_seq <- as.factor(train_hframe$end_seq)\n# split the data into training and validation sets:\nsplits <- h2o.splitFrame(data = train_hframe, ratio = 0.8, seed = SEED)\ntrain <- splits[[1]]\nvalid <- splits[[2]]\nxgb_x <- h2o.xgboost(x = features,\n        y = 'x_1',\n        nfolds = 5,\n        seed = SEED,\n        booster = 'gbtree',\n        ntrees=NTREES,\n        max_depth=MAX_DEPTH,\n        learn_rate=LRATE,\n        keep_cross_validation_predictions = TRUE,\n        training_frame = train_hframe,\n        validation_frame = valid)\n\nprint(h2o.mse(xgb_x,train=TRUE,valid=TRUE))\n\n# Train model for Y\ntrain_hframe = train_data_clean %>% select(-c(ID,x_1,z_1)) %>% as.h2o() # will need to separate out x_1,y_1,z_1 & ID for each iteration\ntrain_hframe$resname <- as.factor(train_hframe$resname)\ntrain_hframe$begin_seq <- as.factor(train_hframe$begin_seq)\ntrain_hframe$end_seq <- as.factor(train_hframe$end_seq)\n# split the data into training and validation sets:\nsplits <- h2o.splitFrame(data = train_hframe, ratio = 0.8, seed = SEED)\ntrain <- splits[[1]]\nvalid <- splits[[2]]\nxgb_y <- h2o.xgboost(x = features,\n        y = 'y_1',\n        nfolds = 5,\n        seed = SEED,\n        booster = 'gbtree',\n        ntrees=NTREES,\n        max_depth=MAX_DEPTH,\n        learn_rate=LRATE,\n        keep_cross_validation_predictions = TRUE,\n        training_frame = train_hframe,\n        validation_frame = valid)\n\nprint(h2o.mse(xgb_y,train=TRUE,valid=TRUE))\n\n# Train model for Z\ntrain_hframe = train_data_clean %>% select(-c(ID,y_1,x_1)) %>% as.h2o() # will need to separate out x_1,y_1,z_1 & ID for each iteration\ntrain_hframe$resname <- as.factor(train_hframe$resname)\ntrain_hframe$begin_seq <- as.factor(train_hframe$begin_seq)\ntrain_hframe$end_seq <- as.factor(train_hframe$end_seq)\n# split the data into training and validation sets:\nsplits <- h2o.splitFrame(data = train_hframe, ratio = 0.8, seed = SEED)\ntrain <- splits[[1]]\nvalid <- splits[[2]]\nxgb_z <- h2o.xgboost(x = features,\n        y = 'z_1',\n        nfolds = 5,\n        seed = SEED,\n        booster = 'gbtree',\n        ntrees=NTREES,\n        max_depth=MAX_DEPTH,\n        learn_rate=LRATE,\n        keep_cross_validation_predictions = TRUE,\n        training_frame = train_hframe,\n        validation_frame = valid)\n\nprint(h2o.mse(xgb_z,train=TRUE,valid=TRUE))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-06T00:41:34.483445Z","iopub.execute_input":"2025-03-06T00:41:34.484924Z","iopub.status.idle":"2025-03-06T00:48:19.336088Z","shell.execute_reply":"2025-03-06T00:48:19.334311Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Create features on test sequence\ntest_sequences <- feat_eng(test_sequences)\n\n# Merge test_clean we created earlier with test_sequences\ntest_data_clean <- merge(test_clean, test_sequences, by.x=\"ID\",by.y=\"target_id\", all.x = TRUE)\ntest_data_clean$resname <- as.factor(test_data_clean$resname)\ntest_data_clean$begin_seq <- as.factor(test_data_clean$begin_seq)\ntest_data_clean$end_seq <- as.factor(test_data_clean$end_seq)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-06T00:48:19.338643Z","iopub.execute_input":"2025-03-06T00:48:19.340032Z","iopub.status.idle":"2025-03-06T00:48:19.375976Z","shell.execute_reply":"2025-03-06T00:48:19.374197Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"preds_x = h2o.predict(xgb_x,newdata = as.h2o(test_data_clean %>% select(-ID)))\npreds_y = h2o.predict(xgb_y,newdata = as.h2o(test_data_clean %>% select(-ID)))\npreds_z = h2o.predict(xgb_z,newdata = as.h2o(test_data_clean %>% select(-ID)))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-06T00:48:19.378547Z","iopub.execute_input":"2025-03-06T00:48:19.379907Z","iopub.status.idle":"2025-03-06T00:48:21.189773Z","shell.execute_reply":"2025-03-06T00:48:21.186923Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# 🤞 Submission","metadata":{}},{"cell_type":"markdown","source":"### For now, we will use the same predictions for all x,y,z combinations.  Future version will incorporate different methods for these predictions.","metadata":{}},{"cell_type":"code","source":"submission <- data.frame(test_data_clean$ID, \n                         test_data_clean$resname, \n                         test_data_clean$resid) \n\nsubmission <- cbind(submission, \n                         as.data.frame(preds_x),\n                         as.data.frame(preds_y),\n                         as.data.frame(preds_z),\n                    \n                         as.data.frame(preds_x),\n                         as.data.frame(preds_y),\n                         as.data.frame(preds_z),\n                    \n                         as.data.frame(preds_x),\n                         as.data.frame(preds_y),\n                         as.data.frame(preds_z),\n                    \n                         as.data.frame(preds_x),\n                         as.data.frame(preds_y),\n                         as.data.frame(preds_z),\n                    \n                         as.data.frame(preds_x),\n                         as.data.frame(preds_y),\n                         as.data.frame(preds_z))\ncolnames(submission) <- sample_submission %>% colnames()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-06T00:48:21.193985Z","iopub.execute_input":"2025-03-06T00:48:21.196388Z","iopub.status.idle":"2025-03-06T00:48:21.582721Z","shell.execute_reply":"2025-03-06T00:48:21.579926Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### My first two submissions in R errored at the scoring phase.  This was due to the ID columns not being in the same order as the sample_submission, so we will create a simple sort column based off the sample_submission, and apply that to our submission df.","metadata":{}},{"cell_type":"code","source":"# Update ID to match submission ID column\nsubmission <- submission %>% mutate(ID = paste(ID,resid,sep=\"_\"))\n\n# Create a sort_order column within the sample_submission so we can apply it to our submission\nsample_submission$sort_order <- seq(1:nrow(sample_submission))\nsubmission <- merge(submission,sample_submission %>% select(ID,sort_order), by='ID',all.x=TRUE) %>% arrange(sort_order) %>% select(-sort_order)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-06T00:48:21.586654Z","iopub.execute_input":"2025-03-06T00:48:21.58914Z","iopub.status.idle":"2025-03-06T00:48:21.645747Z","shell.execute_reply":"2025-03-06T00:48:21.642775Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"submission %>% head(3)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-06T00:48:21.649972Z","iopub.execute_input":"2025-03-06T00:48:21.652358Z","iopub.status.idle":"2025-03-06T00:48:21.700928Z","shell.execute_reply":"2025-03-06T00:48:21.69817Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"write_csv(submission, 'submission.csv')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-03-06T00:48:21.70514Z","iopub.execute_input":"2025-03-06T00:48:21.707362Z","iopub.status.idle":"2025-03-06T00:48:21.91548Z","shell.execute_reply":"2025-03-06T00:48:21.913353Z"}},"outputs":[],"execution_count":null}]}