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 two PDFs about performance but I'd still like some feedback.
help("Rprof")and several packages for the R part of your code. \$\endgroup\$lminstead ofglm? \$\endgroup\$speedglmpackage (cran.r-project.org/web/packages/speedglm/index.html). Worth a try. Note how it says "High performances can be obtained especially if R is linked against an optimized BLAS, such as ATLAS". You will find many articles showing you how to do that if you googleR blas atlas. \$\endgroup\$