3
\$\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\$
15
  • 2
    \$\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, 2015 at 13:45
  • \$\begingroup\$ Thank you @Roland ! Will give it a try, I was actually just reading about it :) \$\endgroup\$
    – Bas
    Sep 21, 2015 at 13:47
  • \$\begingroup\$ try lm instead of glm ? \$\endgroup\$
    – Ben Bolker
    Sep 21, 2015 at 13:49
  • 3
    \$\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, 2015 at 23:08
  • 2
    \$\begingroup\$ @Flodel, if you create your comment about the HPC as answer I will accept it, you can't believe how much it helped me :) \$\endgroup\$
    – Bas
    Oct 20, 2015 at 9:03

1 Answer 1

2
\$\begingroup\$

Comment by @flodel:

The CRAN High-Performance Tasks Views mention the speedglm package. It is 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.

I'll point you to these results showing how switching from the default blas shipped with R to OpenBLAS improved this person's qr decomposition (what lm uses) computation times by a factor of ~4 (from 417 to 113 ms). So regardless of whether you choose to try speedglm, it is definitely worth looking into what blas you are currently using and possibly switching to a better one.

\$\endgroup\$
1
  • 2
    \$\begingroup\$ @flodel If you care to repost this answer under your own name, we can delete this placeholder. \$\endgroup\$ Jan 3, 2023 at 18:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.