For my internship I have to perform certain analysis to determine residuals & outliers. The table I'm currently using has over 12million records with 130+ columns.
My first tests take approx. 5398 seconds or 1.5hours to apply the module and establish outliers and put them into a plot. The goal of my project is performing analysis on 14-20 similar models within an hour(on a better server tho).
Right now I've scripted the following code and I was wondering. Where could the code below be improved the most?
<br />(I know I'm not the best coder, but always trying to learn and improve :))
I am currently reading the following pdf's about performance but I'd still like some feedback:)
http://www.burns-stat.com/pages/Tutor/R_inferno.pdf <br />
http://www.bioconductor.org/help/course-materials/2010/BioC2010/EfficientRProgramming.pdf
Thank you in advance,
**Edit1**
I just changed the glm() to lm() and the process time went down from 1.5hr to 15 minutes.
I also got the following details from Rprof()
$by.self
self.time self.pct total.time total.pct
"lm.fit" 920.64 99.16 921.56 99.26
".Call" 1.48 0.16 1.48 0.16
"plot.xy" 1.22 0.13 1.22 0.13
".External2" 1.10 0.12 1.52 0.16
"colnames<-" 0.92 0.10 0.92 0.10
".External" 0.88 0.09 0.88 0.09
"as.character" 0.66 0.07 0.66 0.07
"unique.default" 0.18 0.02 0.24 0.03
"strsplit" 0.18 0.02 0.18 0.02
"[.data.frame" 0.16 0.02 0.18 0.02
"match" 0.12 0.01 0.12 0.01
"sort.int" 0.12 0.01 0.12 0.01
"Determine_bounds" 0.08 0.01 927.86 99.94
"na.omit" 0.08 0.01 0.42 0.05
"FUN" 0.08 0.01 0.34 0.04
"na.omit.data.frame" 0.08 0.01 0.34 0.04
"is.na" 0.08 0.01 0.08 0.01
"is.factor" 0.06 0.01 0.06 0.01
"lapply" 0.04 0.00 0.38 0.04
"anyNA" 0.04 0.00 0.04 0.00
"readLines" 0.04 0.00 0.04 0.00
"model.matrix.default" 0.02 0.00 1.24 0.13
"model.response" 0.02 0.00 0.68 0.07
"anyDuplicated.default" 0.02 0.00 0.02 0.00
"as.double" 0.02 0.00 0.02 0.00
"lazyLoadDBfetch" 0.02 0.00 0.02 0.00
"max" 0.02 0.00 0.02 0.00
"pmatch" 0.02 0.00 0.02 0.00
"table" 0.02 0.00 0.02 0.00
$by.total
total.time total.pct self.time self.pct
"Determine_bounds" 927.86 99.94 0.08 0.01
"main" 927.86 99.94 0.00 0.00
"lm" 924.00 99.53 0.00 0.00
"lm.fit" 921.56 99.26 920.64 99.16
".External2" 1.52 0.16 1.10 0.12
".Call" 1.48 0.16 1.48 0.16
"dbGetQuery" 1.48 0.16 0.00 0.00
"DbQuery" 1.48 0.16 0.00 0.00
"doTryCatch" 1.48 0.16 0.00 0.00
".local" 1.48 0.16 0.00 0.00
".oci.GetQuery" 1.48 0.16 0.00 0.00
"standardGeneric" 1.48 0.16 0.00 0.00
"try" 1.48 0.16 0.00 0.00
"tryCatch" 1.48 0.16 0.00 0.00
"tryCatchList" 1.48 0.16 0.00 0.00
"tryCatchOne" 1.48 0.16 0.00 0.00
"plot" 1.26 0.14 0.00 0.00
"plot.default" 1.26 0.14 0.00 0.00
"model.matrix.default" 1.24 0.13 0.02 0.00
"model.matrix" 1.24 0.13 0.00 0.00
"plot.xy" 1.22 0.13 1.22 0.13
"colnames<-" 0.92 0.10 0.92 0.10
".External" 0.88 0.09 0.88 0.09
"dev.off" 0.84 0.09 0.00 0.00
"model.response" 0.68 0.07 0.02 0.00
"as.character" 0.66 0.07 0.66 0.07
"summaryRprof" 0.52 0.06 0.00 0.00
"<Anonymous>" 0.44 0.05 0.00 0.00
"eval" 0.44 0.05 0.00 0.00
"model.frame.default" 0.44 0.05 0.00 0.00
"na.omit" 0.42 0.05 0.08 0.01
"lapply" 0.38 0.04 0.04 0.00
"FUN" 0.34 0.04 0.08 0.01
"na.omit.data.frame" 0.34 0.04 0.08 0.01
"unique.default" 0.24 0.03 0.18 0.02
"factor" 0.20 0.02 0.00 0.00
"strsplit" 0.18 0.02 0.18 0.02
"[.data.frame" 0.18 0.02 0.16 0.02
"[" 0.18 0.02 0.00 0.00
"quantile" 0.16 0.02 0.00 0.00
"quantile.default" 0.16 0.02 0.00 0.00
"match" 0.12 0.01 0.12 0.01
"sort.int" 0.12 0.01 0.12 0.01
"sort" 0.12 0.01 0.00 0.00
"sort.default" 0.12 0.01 0.00 0.00
"is.na" 0.08 0.01 0.08 0.01
"as.factor" 0.08 0.01 0.00 0.00
"diff" 0.08 0.01 0.00 0.00
".getXlevels" 0.08 0.01 0.00 0.00
"IQR" 0.08 0.01 0.00 0.00
"levels" 0.08 0.01 0.00 0.00
"unique" 0.08 0.01 0.00 0.00
"is.factor" 0.06 0.01 0.06 0.01
"anyNA" 0.04 0.00 0.04 0.00
"readLines" 0.04 0.00 0.04 0.00
"pdf" 0.04 0.00 0.00 0.00
"sapply" 0.04 0.00 0.00 0.00
"anyDuplicated.default" 0.02 0.00 0.02 0.00
"as.double" 0.02 0.00 0.02 0.00
"lazyLoadDBfetch" 0.02 0.00 0.02 0.00
"max" 0.02 0.00 0.02 0.00
"pmatch" 0.02 0.00 0.02 0.00
"table" 0.02 0.00 0.02 0.00
"anyDuplicated" 0.02 0.00 0.00 0.00
"deparse" 0.02 0.00 0.00 0.00
".deparseOpts" 0.02 0.00 0.00 0.00
"paste" 0.02 0.00 0.00 0.00
"range" 0.02 0.00 0.00 0.00
"xy.coords" 0.02 0.00 0.00 0.00
$sample.interval
[1] 0.02
$sampling.time
[1] 928.4
The R code:
#include libraries
library(DBI)
library(ROracle)
library(outliers)
library(Rmisc)
library(data.table)
#set params
options(width=10000, error=traceback)
args <- commandArgs()
#print(args)
debug <-F
### Connect to database
drv <- dbDriver("Oracle")
con <- dbConnect(drv, username="DB USER", password="DB PASS", dbname="DB NAME")
main <- function(){
#Data_table = "BPM_FINISHING_DATA"
Data_table = "TBL_TRAINING_SET"
print(paste("Performing queries against", Data_table))
AllTraits = "BACKFAT_CARCASS"
if(debug){cat("\n\n")}
ptm <- proc.time()
a <- Determine_bounds(Data_table)
tijd <- proc.time() - ptm
cat("\n")
cat(paste("\t",tijd['elapsed'][1]))
}
Determine_bounds <- function(Data_table){
### Create query
sqlstr <- paste("
select PIG_ID,
coalesce(to_char(trial_tstart), to_char(trial_tmed), to_char(trial_tend)) \"trial\",
bln_tend_sex_lgp_itn \"farm_tend\",
gain_tstart_tend \"gain\",
BREEDING_LINE,
SEX_CODE,
bln_cmt_tstart_itn \"farmcomp\",
--coalesce(bln_cmt_tstart_itn, hys_birth_tend, HYS_CGG) \"farmcomp\",
birth_weight_scaled
from
",Data_table,"
where
EBT_TSTART = 'Y'
AND EBT_TEND = 'Y'
")
if(Data_table == "BPM_FINISHING_DATA"){
sqlstr <- paste(sqlstr, " AND FINISHER_YN = 'N'")
}
cat(format(sqlstr))
### Run query against database
res <- DbQuery(sqlstr)
res$trial[is.na(res$trial)] <-"XXX"
res$farm_trend[is.null(res$farm_trend)] <- -9
res$birth_weight_scaled[is.null(res$birth_weight_scaled)] <- -9
lmfit = lm(res$gain ~ res$trial + res$farm_trend * res$BREEDING_LINE* res$SEX_CODE + res$birth_weigh t_scaled)#+res$farmcomp)
kwant <- quantile(resid(lmfit), probs= c(0.25, 0.75))
Q1 <- kwant[1]
Q3 <- kwant[2]
sigma <- IQR(resid(lmfit))
upp_multi <- 3.5 ##Amount of times sigma for outlier calculation
low_multi <- 1.5 ##Amount of times sigma for outlier calculation
upp_fence <- Q3+(upp_multi * sigma)
low_fence <- Q1-(low_multi * sigma)
print(paste("Upper fence: ", upp_fence, " \t Lower Fence: ", low_fence ,"\n"))
pdf('~/scripts/PlotResiduals.pdf')
plot(resid(lmfit),ylim=c(low_fence*2,upp_fence *2), axes=F )
title(main="Outliers", sub="X axis title")
abline(h= upp_fence, col="blue")
abline(h= low_fence, col="red")
dev.off()
}
DbQuery <- function(sqlstr){
res <- dbGetQuery(con, sqlstr)
return(res)
}
main()
dbDisconnect(con)
q()