For my internship I have to perform certain analysis to determine residuals and outliers. The table I'm currently using has over 12 million records with 130+ columns.
My first tests take approx. 5398 seconds or 1.5 hours 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, though).
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()
Right now I've scripted the following code and I was wondering. Where could the code below be improved the most? (I know I'm not the best coder, but always trying to learn and improve.)
I am currently reading [these][1] [two][2] PDFs about performance but I'd still like some feedback.
[1]: http://www.burns-stat.com/pages/Tutor/R_inferno.pdf
[2]: http://www.bioconductor.org/help/course-materials/2010/BioC2010/EfficientRProgramming.pdf