1
\$\begingroup\$

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.

\$\endgroup\$

migrated from stackoverflow.com Sep 21 '15 at 22:47

This question came from our site for professional and enthusiast programmers.

  • \$\begingroup\$ There is a technique called "profiling", which you can use to find out which parts of your code are actually consuming most of the time. Learn that technique. See help("Rprof") and several packages for the R part of your code. \$\endgroup\$ – Roland Sep 21 '15 at 13:45
  • \$\begingroup\$ Thank you @Roland ! Will give it a try, I was actually just reading about it :) \$\endgroup\$ – Bas Sep 21 '15 at 13:47
  • \$\begingroup\$ try lm instead of glm ? \$\endgroup\$ – Ben Bolker Sep 21 '15 at 13:49
  • \$\begingroup\$ or even lm.fit if you only need the residuals \$\endgroup\$ – Roland Sep 21 '15 at 14:03
  • 2
    \$\begingroup\$ The CRAN High-Performance Tasks Views (cran.r-project.org/web/views/HighPerformanceComputing.html) mention the speedglm package (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 google R blas atlas. \$\endgroup\$ – flodel Sep 21 '15 at 23:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.