Tutorial on univariant time series analysis using R for monthly rainfall predictions.

Introduction

When it is monsoon time, news papers are flooded with predictions about average rain fall during upcoming monsoon. I was wondering how it is done. Through a little bit of net surfing I learned about Time Series Analysis (TSA) . Prime reason we use TSA is to perform time series forecasting. It is primarily based on assumption that, we could predict/forecast future values based on past observations. The two approaches used in TSA are 1) Time domain approach and 2) Frequency domain approach. In this post we will limit yourself to the concept of time domain approach especially with uni variant data ( means, we only have one set of observations. This approach, focuses on predicting future values of a time series as a function of previous observations.Let us explore how this can be applied to predict rainfall in a given month based on average monthly rain fall recorded between 1901-2013. Data required for predicting rain fall is obtained from Indian Data repository.

Acquire data and transform

Setup your working environment and read average monthly rainfall data which is stored as a CVS file. Note the row.names option. I set this as first column in my data frame. This is important because we need to convert “monrain” object into a time series object for all other steps.

library(forecast)
monrain <- read.csv(file = "monRain.csv", header = T, row.names = 1)
head(monrain, n = 5) # identify

##       JAN  FEB  MAR  APR  MAY   JUN   JUL   AUG   SEP   OCT  NOV  DEC
## 1901 34.7 38.6 17.8 38.9 50.6 113.2 241.4 271.6 124.7  52.4 38.7  8.2
## 1902  7.4  4.2 19.0 44.1 48.8 111.7 284.9 201.0 200.2  62.5 29.4 25.2
## 1903 16.7  8.0 31.1 17.1 59.5 120.3 293.2 274.0 198.1 119.5 40.3 18.0
## 1904 14.9  9.7 31.4 33.7 73.8 165.5 260.3 207.7 130.8  69.8 11.2 16.4
## 1905 24.7 20.3 41.8 33.8 55.8  93.7 253.0 201.7 178.1  54.9  9.6 10.1

tail(monrain, n = 5)

##       JAN  FEB  MAR  APR  MAY   JUN   JUL   AUG   SEP   OCT  NOV  DEC
## 2009 12.0 12.0 14.2 25.1 56.0  85.7 280.7 192.5 139.4  71.4 53.7 11.1
## 2010  7.5 17.0 14.0 39.0 73.8 138.1 300.7 274.7 197.7  69.0 61.4 22.7
## 2011  6.8 25.8 22.4 41.1 53.1 183.5 246.0 284.9 186.9  38.1 20.1  7.6
## 2012 26.5 12.7 11.3 47.5 31.7 117.8 250.2 262.4 193.5  58.7 30.7 11.7
## 2013 11.3 40.1 15.7 30.4 57.8 219.8 310.0 254.7 152.7 129.4 14.0  6.7

While there are several ways in which you could convert a regular data frame into a time series object. I am using the transpose option. To do so, you need to know when you started first observation and recorded last observation along with the frequency of observations. In my data the first observation is taken on 1901 Jan and last observation is taken on 2013 Dec.at monthly frequency. With this information I will convert the data frame into a time series Object using ts() function from R base package.Examine the ts object with summary and other functions.

tsmonrain <- ts(as.vector(t(as.matrix(monrain))), start = c(1901, 1), end = c(2013, 12), frequency = 12)
str(tsmonrain) # to understand the structure of the tsmonrain data object

##  Time-Series [1:1356] from 1901 to 2014: 34.7 38.6 17.8 38.9 50.6 ...

summary(tsmonrain) # to obtain the summary of observations

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.60   22.80   49.95   98.23  169.80  375.50

start(tsmonrain) # to explore the start of the observation time

## [1] 1901    1

end(tsmonrain) # to explore the end of observations

## [1] 2013   12

frequency(tsmonrain) # to explore the frequency of observations

## [1] 12

Explore data for trend and seasonality

Like any other data analysis, first and foremost thing to perform in any TSA analysis is to explore the data. However, in TSA you will plot your data with specific objective. Primary objective of this is to identify components such as trend, seasonality, and irregularities if any in our data. Why is it so important to find them? Well, if your data has any of those components it is said to be non-stationary and data may not be suitable for forecasting.IN statistical terms you need to fix the variance and mean in our data series.

Let us plot our data with time on x axis and observations(average monthly rainfall) on Y axis and also add a trend line.

plot(tsmonrain, ylab= "Average rain fall (mm/Month)", main = "Monthly average rain fall from 1901-2013", type = "l" )
lines(lowess(tsmonrain), col = "blue")

From the above plot, it appears that our data “does not” have any trend (neither upward nor downward) and stationary. That is good.
But from historic knowledge we know that monsoon rain in India is a seasonal phenomenon. There has to be fluctuations in periodic manner within any given year. Let us plot a small window of the time series data and see if we could pick it. We will use the window() function to subset the data for chosen period.

plot(window(tsmonrain, start= 2000, end = 2002, frequency = 12), type = "b", col = "brown",ylab= "Average rain fall (mm/Month)", main = "Monthly average rain fall from 2000-2002")

Yes indeed we could see the fluctuation in the monthly ran fall. Let us explore this phenomenon a bit more. How does the Average rainfall in a given month across years is distributed. We will use the box plots.

boxplot(tsmonrain~cycle(tsmonrain), notch = T,col = "purple", xlab = "Calender month", ylab = "average rain fall (mm)", main = " Monthly distribution of rain between 1901-2013 ")

From the above plot it is evident that Over past 100 years or so we have received more rain during the second and third quarters or 6, 7, 8 and 9th Months with few outliers here and there. This tells us that there is a seasonal component to our data. So our data is not a true stationary data.

Decompose data to make it stationary

To move ahead with TSA, we need to remove the seasonality component make our data stationary. This process is often referred to as “Decomposition”. There are several ways to decompose your time series data depending on cause of non-stationarity.

Let us decompose data to eliminate the trend, if any, using diff() function and averaging the monthly data. I am only differentiating the data by 4 and 8 months because the plot will not change even after decomposition with Single month or 12 months WRT seasonality.

par(mfrow = c(2,1))
plot(diff(tsmonrain, 4), ylab = "differenced on 4 month mean", main = "Decomposition on 4month Mean") # four month mean
plot(diff(tsmonrain,8), ylab = "differenced on 8 month mean", main = "Decomposition on 8 month Mean") # 8 month mean

From the above plots it is evident that irrespective of differentiating it is hard to find a trend. so WRT to trend the data is stationary.

But we are pretty sure that, our data has the seasonality (From Box plots). So we must decompose our data for seasonality. Let us decompose data to eliminate the variance (seasonality) using the log transformation

plot(log10(tsmonrain))

We can combine diff and log to eliminate both variance and mean in our time series data.

plot(log10(diff(tsmonrain)))

## Warning in plot(log10(diff(tsmonrain))): NaNs produced

Alternatively we could use the Seasonal Trend Decomposition using Loess (stl()) available in base package or decompose().

plot(stl(tsmonrain, s.window = "periodic"), main = "Decomposition of additive time series")

Have a look at this plot. On top you have the original data plotted. plot in panel labelled as seasonal we have the data series plotted after removing the seasonal component. Below that is the data after removal of trend component. The final plot is the one which is called as white noise or residuals.

Model building

Now we could go ahead to model building. The type of models could be AR (Auto regression), MA (Moving average), ARMA, ARIMA seasonal ARIMA etc.. To build models you need to provide, p,d,q for the AR and MA separately for a given time point.To obtain this information we will use the Auto Correlation (ACF()) and Partial Auto-Correlation (pacf()) functions. Partial auto correlation function will provide the lags (p) in AR(p) models and p, d,q values for the ARIMA models. Please refer to the above links for details explanation of the same.

par(mfrow = c(1,2))
acf(log10(tsmonrain), main = "auto correlaiton plot")
pacf(log10(tsmonrain), main = "partial auto correlation plot")

IN above plots you could figure out the lag values to be used in model building if you are using AR(p)/MA or ARMA/ARIMA. Thanks to R community. Here i would like to use the “auto.arima()” function available in the forecast package. One good reason why I picked the package is i don't need to figure these options. The function call it self will try several options and provide me a best model based on the Akaike information criterion (AIC) . so i will pass the log10 transformed data to the auto.arima() function.

# ARIMA models _ Box-Jenkins approach. 1) Model identification, 2) Parameter estimation 3) Diagnostic checking
arimaFit <- auto.arima(log10(tsmonrain), approximation = F, trace = F)
summary(arimaFit)

## Series: log10(tsmonrain) 
## ARIMA(2,0,0)(1,0,0)[12] with non-zero mean 
## 
## Coefficients:
##          ar1      ar2    sar1  intercept
##       0.1301  -0.0601  0.8590     1.7436
## s.e.  0.0304   0.0280  0.0157     0.0464
## 
## sigma^2 estimated as 0.05575:  log likelihood=27.13
## AIC=-44.25   AICc=-44.21   BIC=-18.19
## 
## Training set error measures:
##                         ME      RMSE       MAE       MPE     MAPE
## Training set -0.0001024415 0.2357753 0.1715335 -3.755715 13.42712
##                   MASE        ACF1
## Training set 0.9892733 0.001706884

The seasonal arima (1) [sar1] model looking good, we could proceed to predictions. Explore the training set error messages.

Predict for future

pred <- predict(arimaFit, n.ahead = 48)
pred$pred

##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 2014 1.1292263 1.6404653 1.2766853 1.5190472 1.7590512 2.2576379 2.3859335
## 2015 1.2158589 1.6550068 1.3425244 1.5507102 1.7568707 2.1851505 2.2953547
## 2016 1.2902753 1.6674979 1.3990794 1.5779084 1.7549976 2.1228847 2.2175488
## 2017 1.3541980 1.6782275 1.4476593 1.6012712 1.7533887 2.0693991 2.1507144
##            Aug       Sep       Oct       Nov       Dec
## 2014 2.3126189 2.1217574 2.0599926 1.2303773 0.9554556
## 2015 2.2323784 2.0684308 2.0153756 1.3027464 1.0665920
## 2016 2.1634529 2.0226239 1.9770502 1.3649105 1.1620569
## 2017 2.1042467 1.9832764 1.9441291 1.4183087 1.2440600

plot(tsmonrain, type = "l", xlim=c(2012,2017), xlab = "year", ylab = "rain")
lines(10^(pred$pred), col= "blue")
lines(10^(pred$pred + 2*pred$se),col= "red", type = "c", lty=2)
lines(10^(pred$pred-2*pred$se),col= "green", type = "c", lty=2)

Now let us take a closer look at the 2016 monthly forecasts

#Try some beutification
current <- 10^( window(pred$pred, start = 2016.000, end = 2016.917, frequency = 12))
current

##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 2016  19.51081  46.50481  25.06567  37.83627  56.88498 132.70420 165.02463
##            Aug       Sep       Oct       Nov       Dec
## 2016 145.69777 105.34743  94.85280  23.16917  14.52302

plot(10^(window(pred$pred, start = 2016.000, end = 2016.917, frequency = 12)), type = "b", col = "blue", xlab = "Year.Calender Month", ylab = "Rainfall in mm", main = "Monthly predicted rainfall for 2016")
lines(10^(window(pred$pred + pred$se, start = 2016.000, end = 2016.917, frequency = 12)),col= "brown", type = "c", lty=2)
lines(10^(window(pred$pred - pred$se, start = 2016.000, end = 2016.917, frequency = 12)),col= "red", type = "c", lty=2)

Metabolite Biomarkers for Wheat - Conventional Vs Organic farming systems

While I am sharpening my machine learning skills and looking for non gene expression data from plant sciences, one recent article caught my attention. Authors were attempting to use metabolite data from multi-year, multi-varietal trail conducted with both conventional and Organic farming systems of wheat to predict farming system. It is an excellent research article with nice hypothesis, systematic data collection and analysis. Please refer to the article

Authors concluded article with nice suggestions for further analysis. I thought I could try my machine learning skills and apply Random Forest (rf) modeling with some feature engineering (Authors have done rf modeling already, but said the performance of model was only slightly better than the SVM and not many model tuning parameters were tested.)

Summary points from Supervised Learning/Classification perspective: 1. SVMs trained on entire data set (all years, all culitvars, both treatments) gave an accuracy of 0.9032 at p-value of 1.486e-11) 2. Better accuracies could be obtained while using the subsets of data. Ex. 2007 year data gave an accuracy of 0.9677 p-value of 3.746e-08. 3. Reducing the feature set to verified biological compounds minimizes the risk of systematic errors through compounds minimizes risk of systematic errors through background noise.

One thing I feel would make a big impact in biomarker discovery: Location would play siginificant component as trait is a function of G X E [Genotype by Environment]. So I would love to use the location data to make the biomarkers identified more robust is available. (However for now i am assuming that, all the trials were done at the same location across years.)

Set R environment and read files

setwd("C:/Users/rduvv/Desktop/Ind_Data/wheat")
library(caret)
library(ggplot2)
library(reshape2)
owm <- read.csv(file = "owm_m.csv", header = T)

Let us quickly explore the data:

dim(owm)

## [1] 313  40

colnames(owm)

##  [1] "Growth"                    "Variety"                  
##  [3] "Year"                      "X2.Methylcitrate"         
##  [5] "X2.Phospho.D.glycerate"    "X4.Aminobutanoate"        
##  [7] "Adenosine"                 "Citrate"                  
##  [9] "D.Gluconic.acid"           "D.Glucono.1.5.lactone"    
## [11] "Ethanolamine"              "Fructose"                 
## [13] "Fumarate"                  "Glucose"                  
## [15] "Glycine"                   "Inositol.1.phosphate"     
## [17] "L.Alanine"                 "L.Arginine"               
## [19] "L.Asparagine"              "L.Aspartate"              
## [21] "L.Citrulline"              "L.Isoleucine"             
## [23] "L.Leucine"                 "L.Methionine"             
## [25] "L.Phenylalanine"           "L.Proline"                
## [27] "L.Serine"                  "L.Threonine"              
## [29] "L.Tryptophan"              "L.Valine"                 
## [31] "Lysine"                    "Malate"                   
## [33] "Pantothenate"              "Ribitol"                  
## [35] "Succinate"                 "Sucrose"                  
## [37] "alpha.alpha.Trehalose"     "beta.Alanine"             
## [39] "myo.Inositol"              "trans.4.Hydroxy.L.proline"

table(owm$Year)

## 
## y2007 y2009 y2010 
##   160    16   137

table(owm$Variety)

## 
## Antonius  Caphorn      CCP       DJ      MC2   Probus      RdB    Runal 
##       31       31       29       14       15       31       31       45 
## Sandomir    Scaro   Titlis 
##       29       27       30

table(owm$Growth)

## 
## conventional      organic 
##          156          157

I am going to use the caret package as authors did for modeling

One of the key dimentionality reduction method available in Caret package is near zero variance and zero variance. Let us use this to find if there are any features with zero (or) near zero variance.

nzv <- nearZeroVar(owm [,4:40], saveMetrics = T)
table(nzv$nzv)

## 
## FALSE 
##    37

table(nzv$zeroVar)

## 
## FALSE 
##    37

As you could see there are no features with zero Variance or Ner zero variance (as all the 37 features are false for zero and nearZero Variance)

Now let us look for highly correlated features and eliminate them.

corOWM <- cor(owm[,4:40]) # find correlations between the features 4 through 40
summary(corOWM[upper.tri(corOWM)]) # Look for the correlations.

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -0.3207  0.3298  0.5216  0.5076  0.7034  0.9925

It is clear that there are features that are highly correlated as indicated by Max correlation i.e 0.9925. Now I want to use the findCorrelations funciton to find the colIDs of the features and eliminate them from the data set.

highcor <- findCorrelation(corOWM, cutoff = 0.75)
sort(highcor)

##  [1]  1  3  8 11 12 14 15 19 20 22 23 24 25 27 30 32 35 36 37

By carefully looking through the above file I will remove the columns from my “owm_m.csv” and save as “owm_cor_rem.csv”. Note: This is because I could not figure out how to merge the features and response. I will work on it later. For now this will work. Now read the file into R and store in a variable called owm.data

owm.data <- read.csv(file = "owm_cor_rem.csv", header = T)
dim(owm.data)

## [1] 313  24

Now you could see that we have eliminated Sixteen features that are highly correlated. Let us do the modeling with this enhanced set of features. Removal of those highly correlated features will reduce noise in model and might improve performance.

Data Partition

set.seed(2)
owm.data.inTrain <- createDataPartition(y = owm.data$Growth, p = 0.8, list = F)
owm.data.train <- owm.data[owm.data.inTrain,]
owm.data.test <- owm.data[-owm.data.inTrain,]
dim(owm.data.train)

## [1] 251  24

dim(owm.data.test)

## [1] 62 24

Just check how the data partitioning worked.Training and testing data sets has 251 and 62 observations respectively.

table(owm.data.train$Growth)

## 
## conventional      organic 
##          125          126

table(owm.data.test$Growth)

## 
## conventional      organic 
##           31           31

Nice. Now let us begin modeling. In the research article, LGOCV is used and I would like to use the repeatedcv. Besides this I want to use preProc function.

# define model controls
ctrl.rf <- trainControl(method = "repeatedcv", repeats = 3, classProbs = TRUE, summaryFunction = twoClassSummary, allowParallel = TRUE)
set.seed(351)
#Build the model
owm.rf <- train(owm.data.train$Growth~., data = owm.data.train, method = "rf",tuneLength = 15, trControl = ctrl.rf, preProc = c("center", "scale"))

## Warning in train.default(x, y, weights = w, ...): The metric "Accuracy" was
## not in the result set. ROC will be used instead.

owm.rf # check  the model

## Random Forest 
## 
## 251 samples
##  23 predictor
##   2 classes: 'conventional', 'organic' 
## 
## Pre-processing: centered (33), scaled (33) 
## Resampling: Cross-Validated (10 fold, repeated 3 times) 
## Summary of sample sizes: 225, 225, 226, 227, 226, 225, ... 
## Resampling results across tuning parameters:
## 
##   mtry  ROC        Sens       Spec       ROC SD      Sens SD   
##    2    0.9414544  0.8940171  0.8863248  0.04076668  0.07222488
##    4    0.9487673  0.9123932  0.8916667  0.03720240  0.06984316
##    6    0.9493740  0.9153846  0.8914530  0.03684804  0.06518537
##    8    0.9522929  0.9230769  0.8863248  0.03547748  0.06718917
##   10    0.9491515  0.9205128  0.8779915  0.03721082  0.06241365
##   13    0.9482406  0.9228632  0.8858974  0.03904680  0.06353603
##   15    0.9480290  0.9200855  0.8777778  0.04103770  0.06560428
##   17    0.9452683  0.9175214  0.8696581  0.04272970  0.06697231
##   19    0.9460970  0.9096154  0.8666667  0.03911553  0.07363573
##   21    0.9434706  0.9147436  0.8777778  0.04341363  0.06871451
##   24    0.9422317  0.9094017  0.8696581  0.04633356  0.07456094
##   26    0.9413578  0.9040598  0.8613248  0.04749999  0.07334910
##   28    0.9392101  0.9036325  0.8566239  0.04914945  0.08040815
##   30    0.9384937  0.9091880  0.8670940  0.05028611  0.07173966
##   33    0.9405462  0.9038462  0.8617521  0.04888764  0.08537425
##   Spec SD   
##   0.09733191
##   0.08233865
##   0.09953435
##   0.10011631
##   0.09659468
##   0.10370989
##   0.10852397
##   0.11432777
##   0.10666004
##   0.09759360
##   0.09509934
##   0.11179812
##   0.10145962
##   0.10779033
##   0.10077579
## 
## ROC was used to select the optimal model using  the largest value.
## The final value used for the model was mtry = 8.

plot(owm.rf) # plot the model

x<-varImp(owm.rf) # Ccheck Variable importance
plot(x, top = 20) # plot Variable importance

The model #8 is the best model at ROC of 0.95229 with Sens of 0.9230 and Spec of 0.8863.

And the metabolites ** myo.Inositol, L.Tryptophan, Phospho-D-Glycerate, Aminobutanoate, LArginine,** seems to contribute to the more than 60 % of model performance.

Now validate the model with test data set.

owm.rf.predicted <- predict (owm.rf, owm.data.test)
confusionMatrix(owm.data.test$Growth, owm.rf.predicted)

## Confusion Matrix and Statistics
## 
##               Reference
## Prediction     conventional organic
##   conventional           29       2
##   organic                 1      30
##                                          
##                Accuracy : 0.9516         
##                  95% CI : (0.865, 0.9899)
##     No Information Rate : 0.5161         
##     P-Value [Acc > NIR] : 5.105e-14      
##                                          
##                   Kappa : 0.9032         
##  Mcnemar's Test P-Value : 1              
##                                          
##             Sensitivity : 0.9667         
##             Specificity : 0.9375         
##          Pos Pred Value : 0.9355         
##          Neg Pred Value : 0.9677         
##              Prevalence : 0.4839         
##          Detection Rate : 0.4677         
##    Detection Prevalence : 0.5000         
##       Balanced Accuracy : 0.9521         
##                                          
##        'Positive' Class : conventional   
## 

Cool. Note the Model accuracy 0.9516 and p- value of 5.105e-14. Check the ROC curve.

#caret's ROC curve procedure
library(pROC)
roc0 <- roc(owm.data.test$Growth, predict (owm.rf, owm.data.test, type = "prob")[,1], levels = rev(levels(owm.data.test$Growth)))
plot(roc0, print.thres = c(0.5), type = "S", print.thres.pattern = "%.3f(Spec = %.2f, Sens = %.2f)", print.thres.cex = .8, legacy.axes =T, col = '#0080ff', grid = TRUE)

## 
## Call:
## roc.default(response = owm.data.test$Growth, predictor = predict(owm.rf,     owm.data.test, type = "prob")[, 1], levels = rev(levels(owm.data.test$Growth)))
## 
## Data: predict(owm.rf, owm.data.test, type = "prob")[, 1] in 31 controls (owm.data.test$Growth organic) < 31 cases (owm.data.test$Growth conventional).
## Area under the curve: 0.9698

Nice to see the model perfomance with 0.9516 Accuracy over 0.9667 Sensitivity and 0.9375 Sensitivity.This looks randomForest could be used to refine model and our effrots of feature selection have paid off with slight increase in model accuracy to 0.9516 (as compared to 0.9032 from paper-across years).

So the top three Meabolites (80% VarImp) that could distinguish Farming system in the Varieties tested in this study are :

Myo.Inositol L.Tryptophan Phospho.D.Glycerate

melted <- melt(owm.data)
head(melted)

##    Growth  Variety  Year         variable      value
## 1 organic Antonius y2007 X2.Methylcitrate 0.01343727
## 2 organic Antonius y2007 X2.Methylcitrate 0.01374540
## 3 organic Antonius y2007 X2.Methylcitrate 0.01023164
## 4 organic Antonius y2007 X2.Methylcitrate 0.01012037
## 5 organic Antonius y2007 X2.Methylcitrate 0.01206357
## 6 organic Antonius y2007 X2.Methylcitrate 0.01116994

small.melted <- subset(melted, melted$variable == 'myo.Inositol' | melted$variable == 'L.Tryptophan')
qplot(small.melted$variable, y = small.melted$value, data = small.melted, color = Growth, facets = ~Variety, geom = "boxplot")

Look at the box plots for individual metabolite summary statistic except for few cultivars these two metabolites can different between Organic and conventional farming systems. Hence I would recommend use of Myo-Inositol and L.Tryptophan as bio markers. Well these are also reported in research article. Only improvement is model accuracy.

There are still some more parameters that could be fine tuned. I will post them at a later date if time permits.

But one information I feel that could significantly enhance the Biomarker discovery could be the location information.

Disclaimer: Observations made here are purely based on my analysis and understanding of data. Any applications derived based on this should be carefully validated experimentally.

Thanks for reading. Your suggestions and feedback are more than welcome.