S-shaped data: Smoothing with quasibinomial distribution

Figure 1: Synthetic data and fitted curves.

S-shaped distributed data can be found in many applications. Such data can be approximated with logistic distribution function [1].  Cumulative distribution function of logistic distribution function is a logistic function, i.e., logit.

To demonstrate this, in this short example, after generating a synthetic data, we will fit quasibinomial regression model to different observations.

ggplot [2], an implementation of grammar of graphics [3], provides capability to apply regression or customised smoothing onto a raw data during plotting.

Generating Synthetic Data

Let generate set of $n$ observation over time $t$, denoted, $X_{1}, X_{2}, ..., X_{n}$ for $k$ observation $X=(x_{1}, x_{2}, ..., x_{k})$. We will use cumulative function for logistic distribution [4],
$$F(x;\mu,s) = \frac{1}{2} + \frac{1}{2} \tanh((x-\mu)/2s)$$, adding some random noise to make
it realistic.

Let's say there are $k=6$ observations with the following parameter sets, $\mu = \{9,2,3,5,7,5\}$  and $s=\{2,2,4,3,4,2\}$, we will utilise `mapply` [5] in generating a syntetic data frame.

generate_logit_cdf <- function(mu, s, 
                               sigma_y=0.1, 
                               x=seq(-5,20,0.1)) {
  x_ms <- (x-mu)/s
  y    <- 0.5 + 0.5 * tanh(x_ms)  
  y    <- abs(y + rnorm(length(x), 0, sigma_y))
  ix   <- which(y>=1.0)
  if(length(ix)>=1) { 
    y[ix] <- 1.0
  }
  return(y)
}
set.seed(424242)
x      <- seq(-5,20,0.025) # 1001 observation
mu_vec <- c(1,2,3,5,7,8)   # 6 variables
s_vec  <- c(2,2,4,3,4,2)
# Syntetic variables
observations_df<- mapply(generate_logit_cdf, 
                              mu_vec, 
                              s_vec, 
                              MoreArgs = list(x=x))
# Give them names
colnames(observations_df) <- c("Var1", "Var2", "Var3", "Var4", "Var5", "Var6")
head(observations_df)  

Smoothing of observations

Using the syntetic data we have generated, `observations_df`,
we can noq use `ggplot` and `quasibinomial` `glm` to visualise
and smooth the variables.

library(ggplot2)
 library(reshape2)
 df_all <- reshape2:::melt(observations_df)
 colnames(df_all) <- c("x", "observation", "y")
 df_all$observation <- as.factor(df_all$observation)
 p1<-ggplot(df_all, aes(x=x, y=y, colour=observation)) + geom_point() + 
        scale_color_brewer(palette = "Reds") +
        theme(
              panel.background = element_blank(), 
              axis.text.x      = element_text(face="bold", color="#000000", size=11),
              axis.text.y      = element_text(face="bold", color="#000000", size=11),
              axis.title.x     = element_text(face="bold", color="#000000", size=11),
              axis.title.y     = element_text(face="bold", color="#000000", size=11)
#              legend.position = "none"
              )
 l1<-ggplot(df_all, aes(x=x, y=y, colour=observation)) +
        geom_point(size=3) + scale_color_brewer(palette = "Reds") + 
        scale_color_brewer(palette = "Reds") + 
        #geom_smooth(method="loess", se = FALSE, size=1.5) +
        geom_smooth(aes(group=observation),method="glm", family=quasibinomial(), formula="y~x",
                       se = FALSE, size=1.5) +
        xlab("x") + 
        ylab("y") +
        #scale_y_continuous(breaks=seq(0.0,1,0.1)) +
        #scale_x_continuous(breaks=seq(0.0,230,20)) +
        #ggtitle("")  + 
        theme(
              panel.background = element_blank(), 
              axis.text.x      = element_text(face="bold", color="#000000", size=11),
              axis.text.y      = element_text(face="bold", color="#000000", size=11),
              axis.title.x     = element_text(face="bold", color="#000000", size=11),
              axis.title.y     = element_text(face="bold", color="#000000", size=11)
              )
 library(gridExtra)
 gridExtra:::grid.arrange(p1,l1)

References
[1] https://en.wikipedia.org/wiki/Logistic_distribution#Applications
[2] http://www.ggplot.org
[3] The Grammar of Graphics, L. Wilkinson, http://www.amzn.com/038724544
[4] http://en.wikipedia.org/wiki/Logistic_distribution#Cumulative_distribution_function.
[5] https://stat.ethz.ch/R-manual/R-devel/library/base/html/mapply.html

Notes:
* Quasibinomial distribution is as a term used in R's GLM implementation context, see here.

Scale back or transform back multiple linear regression coefficients: Arbitrary case with ridge regression

Summary

The common case in data science or machine learning applications, different features or predictors manifest them in different scales. This could bring difficulty in interpreting the resulting coefficients of linear regression, such as one feature having very large or small values compare to other predictors and being in different units first of all. The common approach to overcome this to use z-score's for each features, centring and scaling.

This approach would allow us to interpret the effects.  A possible question however is how could we map regression coefficients obtained with the scaled data back to original coefficients. In the context of ridge regression, this question is posed by Mark Seeto in the R mailing list and provided a solution for two predictor case with an R code. In this post we formalize his approach. Note that, in the case of how to scale, Professor Gelman suggests dividing them by two standard deviations. In this post we won't cover that approach and use usual approach.

Algebraic Solution: No error term

An arbitrary linear regression for $n$ variable reads as follows
$$y=(\Sigma_{i=1}^n   \beta_{i} x_{i}) + \beta_{0}$$
here, $y$ is being response variable, $x_{i}$ are the predictors, $n=1,..,n$. Let's use primes for the scaled regression equation  for $n$ variable.
$$y'=(\Sigma_{i=1}^n   \beta_{i}' x_{i}') + \beta_{0}'$$
We would like to express $\beta_{i}$ by only using $\beta_{i}'$ and two statistic from the data, namely mean and standard deviations, $\mu_{x_{i}}$, $\mu_{y}$, $\sigma_{x_{i}}$ and $\sigma_{y}$.

The following transformation can be shown by using the z-scores and some algebra,

$$\beta_{0}=\beta_{0}' \sigma_{y} + \mu_{y} - \Sigma_{i=1}^{n} \frac{\sigma_{y}}{\sigma_{x_{i}}}\beta_{i}' \mu_{x_{i}}$$
$$\beta_{i} = \beta_{i}' \frac{\sigma_{y}}{\sigma_{x_{i}}}$$ 

Ridge regression in R

There are many packages and tools in R to perform ridge regression. One of the prominent one is glmnet. Following Mark Seeto's example, here we extent that in to many variate case with a helper function scaleBack.lm from R1magic package.  Function provides a transform utility for $n$-variate case. Here we demo this using 6 predictors, also available as gist,

rm(list=ls())
library(glmnet)
library(R1magic) # https://github.com/msuzen/R1magic
set.seed(4242)
n <- 100 # observations
X <- model.matrix(~., data.frame(x1 = rnorm(n, 1, 1),
                                 x2 = rnorm(n, 2, 2),
                                 x3 = rnorm(n, 3,2),
                                 x4 = rnorm(n, 4,2),
                                 x5 = rnorm(n, 5,1),
                                 x6 = rnorm(n, 6,1)
                                ))[,-1] # glmnet adds the intercept
Y          <- matrix(rnorm(n, 1, 2),n,1)
# Now apply scaling 
X.s        <- scale(X)
Y.s        <- scale(Y)
# Ridge regression & coefficients with scaled data
glm.fit.s    <- glmnet(X.s, Y.s, alpha=0)
betas.scaled <- as.matrix(as.vector(coef(glm.fit.s)[,80]), 1, 7)
# trasform the coefficients 
betas.transformed <- scaleBack.lm(X, Y, betas.scaled)
# Now verify the correctness of scaled coefficients: 
# ridge regression & coefficients
glm.fit    <- glmnet(X, Y, alpha=0)
betas.fit  <- as.matrix(as.vector(coef(glm.fit)[,80]), 1, 7)
# Verify correctness: Difference is smaller than 1e-12
sum(betas.fit-betas.transformed) < 1e-12 # TRUE

Conclusion

Multiple regression is used by many practitioners. In this post we have shown how to scale continuous predictors and transform back the regression coefficients to original scale. Scaled coefficients would help us to better interpret the results. The question of when to standardize the data is a different issue.

Distributed processing with BASH via LSF' blaunch

LSF is currently a de-facto standard in managing computing queues. Many inexpensive CPU or GPU clusters uses LSF. Majority of current academic HPC codes are using MPI for their core computations. While the parallelisation/distribution of the computational tasks are achived by the core MPI codes, LSF batch job script, usually can invoke via shell, is still a serial code. See here. This may bring a bottleneck in case that the batch script needs to run large set of commands in sequence, that can be done in distributed way.  

Idle CPUs

An example task could be that your parallel code generates large set of output files and you require to further process them, such as compressing all data files produced by processes. Once you launch $n$ processes, after completion of the MPI code, running compression over the $n$ data files will be processed in a single process. This will not only waste $n-1$ cpu time waiting idle, it will waste our wall clock time $n-1$ fold while compression is performed sequentially. To overcome this limitation LSF has introduced a tool called blaunch, for distributing the tasks. The documentation provided by IBM is really fantastic. However it lacks a simple example for novice users. In this post I provide a simple example that compresses $n$ data files using $n$ distributed processes. 

Job distribution via blaunch

The first step in doing this task, we need to figure our what hostnames we are currently running with the LSF job script. This can be obtained using $LSF_HOSTS environment variable. The second step is to write this information into hosts file that can be read by blaunch. A compression script should be provided in a seperate file. This is an important point. Because whole concept of blaunch is that it distributes the given command to all available processes. Note that, LSF process is running your script not the distributed host processes, at least what it is transparently. Knowing this fact, we can only differentiate the process ID outside of our job submit script using $LSF_PM_TASKID.

Here is a possible portion of the job script, using BASH,

1 2 3 4 5

cd /mydata rm hosts.txt export myhosts=($LSB_HOSTS) for hostName in ${myhosts[*]};do echo $hostName >> hosts.txt ; done blaunch -u hosts.txt compress.sh

As you can see we don't use  $LSF_PM_TASKID in the job script, rather in the compression script that is going to be launch in one of the distributed process:

1 2 3 4 5

export trajF=(`find . -name '*.dat'`) # ${trajF[ii]}; ii=`expr $ii + 1`;done # This line should be removed (see comments) export ii=`expr $LSF_PM_TASKID - 1` echo ${trajF[$ii]} xz -4e --format=lzma ${trajF[$ii]}

We assume that commands are available in the BASH environment. Note that we use bash array to identify file names and your output convention should match with the identification.

 Conclusion

We have provided a simple working example of launching parallel shell script tasks using blaunch. A similar approach is presented here. An alternating approach to this could be using GNU Parallel , which is discussed here.

Euclid Algorithm for Set of Integers: 'Reduce' vs. trees in R

The Euclid Algorithm provides a solution to the greatest common divisor (GCD) of two natural numbers $x_{1}$ and $x_{-2}$, denoted by $GCD(x_{1}, x_{2})$. This will produce the largest integer that divides $x_{1}$ and $x_{2}$. Solution is proposed by Euclid of Ancient Greece. This can be formulated as a recurrence relation:
$$x_{k} = x_{k-1}  modulo  x_{k-2}.$$
'modulo' binary operation returns the remainder from a given division. Stopping criterion for the recurrence relation is reached when $x_{k-2}=0$ and the result of GCD will be the current value of $x_{k-1}$. This process can be visualised as successive divisions. Let's implement this in R in a naive way.

# Naive Euclid algorithm by msuzen
gcd <- function(a, b)
{
  rk_1 <- a;
  rk_2 <- b;
  # Recurrence Formula:  r_k =  r_k-1 modulo r_k-2
  # Increment k until r_k-2 == 0 
  while(rk_2 != 0) {
    rk      <- rk_1%%rk_2; # remainder
    rk_1    <- rk_2;       # proceed in recurrence
    rk_2    <- rk;
  }
  return(rk_1)
}

This is a straight forward task. Let's make the problem little more generic. What happends if we would like to know GCD of $n$ natural numbers, $x_{1},..., x_{n}$? Than, a solution is to apply GCD operation pairwise, for example if $n=3$:
$$GCD(x_{1}, GCD(x_{2}, x_{3})) = GCD(GCD(x_{1}, x_{2}), x_{3})$$
How can we implement this for a vector of non-negative integers?

Tree Approach

The simplest way to reach GCD of $n$ numbers is probably thinking of this process as a binary tree, formed by pairing elements of set of integers as we obtain GCDs. It is relatively easy to implement this because ordering of pairs is not important. We can start from the beginning and pair up as we obtain the results. Here is the naive implementation.

1 2 3 4 5 6 7 8

gcdN <- function(X) { n <- length(X) gcdEnd <- X[1] for(i in 2:n) { gcdEnd <- gcd(gcdEnd, X[i]) } return(gcdEnd) }

'Reduce' operation 

So called tree approach we have given above is actually noting but a Reduce operation in the context of MapReduce. The function gcdN can be replaced with a single line.

Reduce("gcd", X)

This example looks trivial but having Reduce and friends in our programming toolbox makes our life a little easier. Noting so novel here! But it could saves us time not to implement tree like algorithm.