{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Linear Regression Tutorial in R\nThis is a chance for me to practice with R. It probably won't give a very good prediction score but it should help me understand the relationships in the data better."},{"metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","trusted":true},"cell_type":"code","source":"library(tidyverse) # metapackage with lots of helpful functions\n\nlist.files(path = \"../input\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Load the data\nUnfortunately the training data exists in two different files, train.csv and structures.csv, so they have to be loaded and pre-processed a bit before the dataset is ready."},{"metadata":{"trusted":true},"cell_type":"code","source":"data <- read.csv(\"../input/train.csv\")\nstructures <- read.csv(\"../input/structures.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Code from https://www.kaggle.com/stoney71/eda-in-r\n\nnames(structures) <- c(\"molecule_name\",\n                          \"atom_index_0\",\n                          \"atom_0\",\n                          \"ai0_str_x\",\n                          \"ai0_str_y\",\n                          \"ai0_str_z\")\n# merge structures data with first atom in the row (\"atom_index_0\")\ndata <- left_join(data,\n                      structures,\n                      by = c(\"molecule_name\", \"atom_index_0\"))\n\nnames(structures) <- c(\"molecule_name\",\n                          \"atom_index_1\",\n                          \"atom_1\",\n                          \"ai1_str_x\",\n                          \"ai1_str_y\",\n                          \"ai1_str_z\")\n# merge structures data with first atom in the row (\"atom_index_1\")\ndata <- left_join(data,\n                      structures,\n                      by = c(\"molecule_name\", \"atom_index_1\"))\n\n# calculate distance between each pair of atoms\ndata <- mutate(data,\n                 coupling_distance = sqrt((ai1_str_x - ai0_str_x)^2 +\n                                          (ai1_str_y - ai0_str_y)^2 +\n                                          (ai1_str_z - ai0_str_z)^2))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"head(data)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Visualize\nFirst, I'll look at the correlations among the quantitative variables. I have this handy code which I found on the internet somewhere that plots the pair plot with R values for the correlations. If any variables are correlated, we will see it here."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Correlation panel\npanel.cor <- function(x, y){\n  usr <- par(\"usr\"); on.exit(par(usr))\n  par(usr = c(0, 1, 0, 1))\n  r <- round(cor(x, y), digits=2)\n  txt <- paste0(\"R = \", r)\n  cex.cor <- 0.8/strwidth(txt)\n  text(0.5, 0.5, txt, cex = cex.cor * r)\n}\n# Customize upper panel\nupper.panel<-function(x, y){\n  points(x,y, pch = 19)\n}\n\nsubset = c(1:1000)\npairs(data[subset,c(3,4,8,9,10,12,13,14,15)],\n      lower.panel=panel.cor,\n      upper.panel=upper.panel,\n      cex.labels=0.7)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Doesn't look like there are any significant correlations so we can include all of the variables in our analysis if needed. "},{"metadata":{},"cell_type":"markdown","source":"# Fit model\nWe will fit a regular linear regression model with no interaction and see what the coefficients look like. I'll start by using only the coupling_distance and type."},{"metadata":{"trusted":true},"cell_type":"code","source":"TestSamples = sample(1:nrow(data),size = 0.20*nrow(data))\ndata.val = data[TestSamples,]\ndata.train = data[-TestSamples,]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"lm.fit=lm(scalar_coupling_constant~coupling_distance+type,data=data.train)\nsummary(lm.fit)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"You can see from the summary that this basic analysis gives us a pretty good fit (R^2 = 0.95). All of the coefficients are significant."},{"metadata":{},"cell_type":"markdown","source":"Let's get a baseline performance metric."},{"metadata":{"trusted":true},"cell_type":"code","source":"mae <- function(error)\n{\n    mean(abs(error))\n}\n\nlm.pred = predict.lm(lm.fit,newdata=data.val)\nactual = data.val$scalar_coupling_constant\n\nprint(\"baseline\")\nlog(mae(lm.pred-actual))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This would put you near the bottom of the leaderboard but at least you understand the data! Let's see if we can improve using more variables."},{"metadata":{"trusted":true},"cell_type":"code","source":"head(data)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"lm.fit=lm(scalar_coupling_constant~atom_index_0+atom_index_1+type+ai0_str_x+ai0_str_y+ai0_str_z+atom_1+ai1_str_x+ai1_str_y+ai1_str_x+coupling_distance,data=data.train)\nsummary(lm.fit)\nlm.pred = predict.lm(lm.fit,newdata=data.val)\nactual = data.val$scalar_coupling_constant\n\nlog(mae(lm.pred-actual))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Nope, in fact our prediction gets slightly worse! Let's try some interaction."},{"metadata":{"trusted":true},"cell_type":"code","source":"lm.fit=lm(scalar_coupling_constant~\n          coupling_distance\n          +type\n          +atom_index_0\n          +atom_index_1\n          +coupling_distance*atom_index_0\n          +coupling_distance*atom_index_1\n          +type*coupling_distance # This one matters big time\n          +type*atom_index_0 # good, as long as type*atom_index_1 also included\n          +type*atom_index_1 # good\n          +atom_index_0*atom_index_1 # good\n          ,data=data.train)\n\nsummary(lm.fit)\nlm.pred = predict.lm(lm.fit,newdata=data.val)\nactual = data.val$scalar_coupling_constant\n\nlog(mae(lm.pred-actual))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There's definitely some interaction! But guessing every combination is arduous, that's why God invented neural networks. Let me know if I made any mistakes, I'm still learning R!"},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"R","language":"R","name":"ir"},"language_info":{"mimetype":"text/x-r-source","name":"R","pygments_lexer":"r","version":"3.4.2","file_extension":".r","codemirror_mode":"r"}},"nbformat":4,"nbformat_minor":1}