liureg : A Comprehensive R Package for the Liu Estimation of Linear Regression Model with Collinear Regressors

The Liu regression estimator is now a commonly used alternative to the conventional ordinary least squares estimator that avoids the adverse effects in the situations when there exists a considerable degree of multicollinearity among the regressors. There are only a few software packages available for estimation of the Liu regression coefficients, though with limited methods to estimate the Liu biasing parameter without addressing testing procedures. Our liureg package can be used to estimate the Liu regression coefficients utilizing a range of different existing biasing parameters, to test these coefficients with more than 15 Liu related statistics, and to present different graphical displays of these statistics.


Introduction
For data collected either from a designed experiment or from an observational study, the ordinary least square (OLS) method does not provide precise estimates of the effect of any explanatory variable (regressor) when regressors are interdependent (collinear with each other).Consider a multiple linear regression (MLR) model, where y is an n × 1 vector of observations on dependent variable, X is known design matrix of order n × p, β is a p × 1 vector of unknown parameters, and ε is an n × 1 vector of random errors with mean zero and variance σ 2 I n , where I n is an identity matrix of order n.
The OLS estimator (OLSE) of β is given by β = (X X) −1 X y, which depends on the characteristics of the matrix X X.If X X is ill-conditioned (near dependencies among various regressors of X X exist) or det(X X) ≈ 0, then the OLS estimates are sensitive to a number of errors, such as non-significant or imprecise regression coefficients (Kmenta, 1980) with wrong sign and non-uniform eigenvalues spectrum.Moreover, the OLS method, for example, can yield a high variance of estimates, large standard errors, and wide confidence intervals.
Researchers may be tempted to eliminate regressor(s) causing problems by consciously removing regressor from the model or by using some screening method such as stepwise and best subset regression etc.However, these methods may destroy the usefulness of the model by removing relevant regressor(s) from the model.To control variance and instability of the OLS estimates, one may regularize the coefficients, with some regularization methods such as the ridge regression (RR), Lasso regression and Liu regression (LR) methods etc., as alternative to the OLS.Computationally, the RR ( βr = (X X + kI) −1 X y) suppresses the effects of collinearity and reduces the apparent magnitude of the correlation among regressors in order to obtain more stable estimates of the coefficients than the OLS estimates and it also improves the accuracy of prediction (see Hoerl and Kennard, 1970;Montgomery and Peck, 1982;Myers, 1986;Rawlings et al., 1998;Seber and Lee, 2003;Tripp, 1983, etc.).However, the ridge coefficient is a complicated function of k when some popular methods (such as given in Golub et al. (1979), Mallows (1973) and McLeod and Xu (2017) etc.) are used for (optimal) selection of k.Different applications can yield values for k which are too small to correct the problem of the ill-conditioned product, X X.In such cases, the RR may still be unstable.Similarly, the choice of k belongs to the researcher, there being no consensus regarding how to select optimal k.As such, other innovative methods were needed to deal with collinear data.Liu (1993) proposed another biased estimator to mitigate the collinearity effect on regressors.They also discussed some of the properties and methods for suitable selection of biasing parameter used in LR.For further detail, see Section "Liu regression estimator." We have developed the liureg (Imdadullah and Aslam, 2017) package to provide the functionality of Liu related computations.The package provides the most complete suite of tools for the LR available in R. Table 1 provides a comparison with other alternatives.For package development and R documentation, we followed Wickham (2015); Leisch (2008); R Core Team (2015).The ridge package by Cule and De Iorio (2012), lmridge by Imdadullah and Aslam (2016a) and lm.ridge from the MASS by Venables and Ripley (2002)   In the available literature, there are only two R packages capable of estimating and/or testing of the Liu coefficients.The R packages mentioned in Table 1 are compared with our liureg package.The lrmest package (Dissanayake and Wijekoon, 2016) computes different estimates such as the OLS, ordinary ridge regression (ORR), Liu estimator (LE), LE type-1, 2, 3, adjusted Liu estimator (ALTE), and their type-1, 2, 3 etc.Moreover, lrmest provides scalar mean square error (MSE), prediction residual error sum of squares (PRESS) values of some of the estimators.The testing of ridge coefficient is performed only on scalar k, however, for a vector of d, the function liu() of lrmest package returns only MSE along with value of the biasing parameter used.The ltsbase package (Kan et al., 2013) computes ridge and Liu estimates based on the least trimmed squares (LTS) method.The MSE value from four regression models can be compared graphically if the argument plot=TRUE is passed to the ltsbase() function.There are three main functions, (i) ltsbase() computes the minimum MSE values for six methods: OLS, ridge, ridge based on LTS, LTS, Liu, and Liu based on LTS method for sequences of biasing parameters ranging from 0 to 1, (ii) the ltsbaseDefault() function returns the fitted values and residuals of the model having minimum MSE, and (iii) the ltsbaseSummary() function returns the regression coefficients and the biasing parameter for the best MSE among the four regression models.
It is important to note that the ltsbase package displays these statistics for models having minimum MSE (bias and variance are not displayed in their output), while our package, liureg, computes these and all other statistics not only for scalar but also for vector biasing parameter.This paper outlines the collinearity detection methods available in the existing literature and uses the mctest (Imdadullah and Aslam, 2016b) package through an illustrative example.To overcome the issues of the collinearity effect on regressors a thorough introduction to Liu regression, properties of the Liu estimator, different methods for the selecting values of d, and testing of the Liu coefficients is presented.Finally, estimation of the Liu coefficients, methods of selecting a biasing parameter, testing of the Liu coefficients, and different Liu related statistics are implemented in R within the liureg package.

Collinearity detection
Diagnosing collinearity is important to many researchers.It consists of two related but separate elements: (1) detecting the existence of collinear relationship among regressors and (2) assessing the extent to which this relationship has degraded the parameter estimates.There are many diagnostic measures used for detection of collinearity in the existing literature provided by various authors (Belsley et al., 1980;Curto and Pinto, 2011;Farrar and Glauber, 1967;Fox and Weisberg, 2011;Gunst and Mason, 1977;Klein, 1962;Koutsoyiannis, 1977;Kovács et al., 2005;Marquardt, 1970;Theil, 1971).These diagnostic methods assist in determining whether and where some corrective action is necessary (Belsley et al., 1980).Widely used, and the most suggested diagnostics, are the value of pair-wise correlations, the variance inflation factor (VIF)/ tolerance (TOL) (Marquardt, 1970), the eigenvalues and eigenvectors (Kendall, 1957), the CN & CI (Belsley et al., 1980;Chatterjee and Hadi, 2006;Maddala, 1988), Leamer's method (Greene, 2002), Klein's rule (Klein, 1962), the tests proposed by Farrar and Glauber (Farrar and Glauber, 1967), the Red indicator (Kovács et al., 2005), the corrected VIF (Curto and Pinto, 2011), and Theil's measures (Theil, 1971), (see also Imdadullah et al. (2016)).All of these diagnostic measures are implemented in a the R package mctest (Imdadullah and Aslam, 2016b).Below, we use the Hald dataset (Hald, 1952), for testing collinearity among regressors.We then use the liureg package to compute the Liu regression coefficients for different Liu related statistics and methods of selection of Liu biasing parameter is performed.For optimal choice of biasing parameter,a graphical representation of the Liu coefficients is considered, along with a bias variance trade-off plot.
In additino, model selection criteria is also performed.The Hald data are about heat generated during setting of 13 cement mixtures of 4 basic ingredients and used by Hoerl et al. (1975).Each ingredient percentage appears to be rounded down to a full integer.The data set is included in both the mctest and liureg packages.The results from all overall collinearity diagnostic measures indicate the existence of collinearity among regressor(s).These results do not tell which regressor(s) are reasons of collinearity.The individual collinearity diagnostic measures can be obtained though:

Collinearity detection: An example
> mctest(x = x, y, all = TRUE, type = "i") Call: imcdiag(x = x, y = y, method = method, corr = FALSE, vif = vif, tol = tol, conf = conf, cvif = cvif, leamer = leamer, all = all) The R Journal Vol.9/2, December 2017 ISSN 2073-4859 All Individual Multicollinearity Diagnostics in 0 or 1 in not detected by the test X1 , X2 , X3 , X4 , coefficient(s) are non-significant may be due to multicollinearity R-square of y on all x: 0.9824 * use method argument to check which regressors may be the reason of collinearity The results from most of the individual collinearity diagnostics suggest that all of the regressors are the reason for collinearity among regressors.The last line of the imcdiag() function's output suggests that method argument should be used to check which regressors may be the reason of collinearity among different regressors.This finding suggest that one should use regularization method such as LR.

Liu regression estimator
To deal with multicollinear data, Liu (1993) formulated a new class of biased estimators that has combined benefits of ORR by Hoerl and Kennard (1970) and the Stein type estimator Stein (1956), βS = c β, where c is parameter 0 < c < 1 to avoid their disadvantages.The Liu estimator (LE) can be defined as, where d is the Liu parameter also known as the biasing (tuning or shrinkage) parameter and lies between 0 and 1 (i.e., 0 ≤ d ≤ 1), I p is the identity matrix of order p × p, and β is OLSE.and other statistical areas, the LE has produced a number of new techniques and ideas, see for example Akdeniz and Kac¸iranlar (2001); Hubert and Wijekoon (2006); Jahufer and Chen (2009, 2011, 2012); Kac¸iranlar et al. (1999); Kac¸iranlar and Sakalho glu (2001); Torigoe and Ujiie (2006).
However, Liu (2011) and Druilhet and Mom (2008) have made statements that the biasing parameter d may lie outside the range given by Liu (1993), that is, it may be less than 0 or greater than 1.The LE is a linear transformation of the OLSE, βd = βols .
the main interest of LE lies in the suitable selection of d for which MSE is minimum and that the efficiency of estimators improves, as compared to other values of d.The βd is named as the LE by Akdeniz and Kac¸iranlar (1995) and Gruber (1998).Liu (1993), in applications to econometrics and engineering, provided some important methods for the selection of d and also provided numerical examples using an iterative minimum MSE method to get the smallest possible value to overcome the problem of collinearity in an effective manner.

Reparameterization
The design matrix X n×p and response variable y n×1 should be standardized, scaled or centered first such that information matrix X X is in the correlation form and vector X y is in the form of the correlation among regressors and the response variable.Consider the regression model, y = β 0 1 + Xβ 1 + ε, where X is centered and 1 = c(1, 1, • • • , 1) .The value for β 0 can be estimated by using be the ordered eigenvalues of the matrix X X and q 1 , q 2 , • • • , q p be the eigenvectors corresponding to their eigenvalues, such that Q = (q 1 , q 2 , • • • , q p ) is an orthogonal matrix of X X and The R Journal Vol.9/2, December 2017 ISSN 2073-4859 therefore, the model can be rewritten in canonical form as y = β 0 1 + Zα + ε, where Z = XQ and Similarly, Eq. 1 can be written in canonical form as, The corresponding estimates of β1 and βd can be obtained by following the relations β1 = Q α and βd = Q αd , respectively.For simplification of notations, X and α will be represented as X and β, respectively.
The fitted values of the LE can be found using Eq. 1, where, H d is LE the matrix (Liu, 1993;Walker and Birch, 1988).It is worthwhile to note that H d is not idempotent because it is not a projection matrix, therefore it is called quasi-projection matrix.
As βd is computed on centered variables, they need to be converted back to the original scale: where S xj is the scaling method of regressors.
The intercept term for the LE ( β0d ) can be estimated using the following relation: (2)

Properties of the Liu estimator
Like the linear RR, the Liu regression is also the most popular method among biased methods, because of its relation to OLS.Its statistical properties have been studied by Akdeniz andKac¸iranlar (1995, 2001), Arslan and Billor (2000), Kac¸iranlar and Sakalho glu (2001), Kac¸iranlar et al. (1999) and Sakalho glu et al. ( 2001) among many others.Due to comprehensive properties of the LE, researchers have been attracted towards this area of research.
Let X j denote the jth column of As already discussed, the regressors are centered, thus, the intercept will be zero and can thereby be removed from the model.However, it can be estimated from relation given in Eq. 2. Table 2, lists the Liu properties that are implemented in our liureg package.
Theoretically and practically, LR is used to propose new methods for the choice of the biasing parameter d to investigate the properties of LR, since the biasing parameter plays a key role while the optimal choice of d is the main issue in this context.In the literature, many methods for the selection of an appropriate biasing parameter d have been studied by Akdeniz and Özkale (2005), Arslan and Billor (2000), Akdeniz et al. (2006), Özkale and Kac¸iranlar (2007), and Liu (1993).

Methods of selecting values of d
The existing methods to select biasing parameter in the LR may not fully address the problem of ill-conditioning when there exists severe multicollinearity, while the appropriate selection of biasing parameter d also remains a problem of interest.The parameter d should be selected when there are improvements in the estimates (have stable estimates) or prediction is improved.
The optimal value of d is one which gives minimum MSE.There is one optimal d for any problem The R Journal Vol.9/2, December 2017 ISSN 2073-4859

1) Linear transformation
The LE is a linear transformation of the OLSE ( βd = F d β) 2) Wide range d Wide range of d have smaller MSE than the OLS 3) Optimal d An optimal d always exists that gives minimum MSE 4) Mean Table 2: Properties of Liu estimator.
by the analogy with the estimate of k in RR, a wide range of d (−∞ < d < 1) can give smaller MSE as compared to that of the OLS.For collinear data, a small change in d varies the LR coefficients rapidly.Therefore, a disciplined way of selecting the shrinkage parameter is required that minimizes the MSE.The biasing parameter d depends on the true regression coefficients (β) and the variance of the residuals σ 2 , unfortunately these are unknown, but they can be estimated from the sample data.
We classified estimation methods as (i) Subjective or (ii) Objective

Subjective methods
In these methods, the selection of d is subjective or of judgmental nature and provides graphical evidence of the effect of collinearity on the regression coefficient estimates and also accounts for variation by the LE as compared to the OLSE.In these methods, the reasonable choice of d is done using the Liu trace and the plotting of bias, variance, and MSE.Like ridge trace, the Liu trace is also a graphical representation of the regression coefficients, βd , as a function of d over the interval (−∞, ∞).
Similarly, the plotting of bias, variance, and MSE from the LE may also be helpful in selecting an appropriate value of d.At the cost of bias, optimal d can be selected at which MSE is minimum.All these graphs can be used for selection of optimal (but judgmental) value of d from the horizontal axis to assess the effect of collinearity on each of the coefficients.These graphical representations do not provide a unique solution, rather they render a vaguely defined class of acceptable solutions.However, these traces are still useful graphical representations to check for some optimal d.

Objective methods
Objective methods, to some extent, are similar to judgmental methods for selection of biasing parameter d, but they require some calculations to obtain these biasing parameters.Table 3 lists widely used methods to estimate the biasing parameter d already available in the existing literature.Table 3 also lists other statistics that can be used for the selection of the biasing parameter d.

Testing of the Liu coefficients
Testing of the Liu coefficients is performed by following Aslam (2014) and Halawa and El-Bassiouni (2000).For testing H 0 : β dj = 0 against β dj = 0, the non-exact t-statistics defined by Halawa and El-Bassiouni (2000) are, where βdj is the jth Liu coefficient estimate and SE( βdj ) is an estimate of standard error, which is the square root of the jth diagonal element of the covariance matrix of LE (see property # 6 in Table 2).
The statistics T dj are assumed to follow Student's t distribution with (n − p) df (Halawa and El-Bassiouni, 2000).Hastie and Tibshirani (1990) and Cule and De Iorio (2012) suggest using the df from (n − trace(H d )).For large sample size, the asymptotic distribution of this statistic is normal (Halawa and El-Bassiouni, 2000).

Sr.# Formula Reference
1) where êd and Table 3: Different available methods to estimate d.
For testing overall significance of vector of LE ( βd ) with E( βd ) = F d β and Cov( βd ), the F-statistic is, The standard error of βd is computed by considering the variance of the estimator, given in Eq. 2, and then taking the square root of this variance, that is: (3) The R package liureg Our R package liureg contains functions related to fitting of the LR model and provides a simple way of obtaining the estimates of LR coefficients, testing the Liu coefficients, and the computation of different Liu related statistics, which prove helpful for selection of optimal biasing parameter d.The package computes different Liu related measures available for the selection of biasing parameter d, and computes value of different biasing parameters proposed by some researchers in the literature.
The liureg objects contain a set of standard methods such as print(), summary(), plot(), and predict().Therefore, inferences can be made easily using the summary method for assessing the estimates of regression coefficients, their standard errors, t-values and their respective p-values.The default function liu which calls liuest() to perform required computations and estimation for given The R Journal Vol.9/2, December 2017 ISSN 2073-4859 values of non-stochastic biasing parameter d.The syntax of default function is, liu(formula,data,scaling=("centered","sc","scaled"),d,...) The four arguments of liu() function are described in Table 4.

Argument Description formula
Symbolic representation for LR model of the form, response ∼ predictors.data Contains the variables that have to be used in LR model.d The biasing parameter, may be a scalar or vector.If a d value is not provided, d = 1 will be used as the default value, i.e., the OLS results will be produced.scaling The methods for scaling of predictors.The centered option, centers the predictors, suggested by Liu (1993), and uses the default scaling option; the sc option scales the predictors in correlation form as described in Belsley (1991); Draper and Smith (1998); and the scaled option standardizes the predictors having mean and unit variance.The liu() function returns an object of class "liu".The functions summary(), dest(), and lstats() etc., are used to compute and print a summary of the LR results, list of biasing parameter by Liu (1993Liu ( , 2011) ) and Liu related statistics such as estimated squared bias, R 2 and variance etc., after bias is introduced in regression model.An object of class "liu" is a list, the components of which are described in Table 5.The scales used to standardize the predictors.xs

Object
The scaled matrix of the predictors.scaling The method of scaling used to standardize the predictors.d The LR biasing parameter(s).

Inter
Whether an intercept is included in the model or not.call The matched call.terms The terms object used.

The Liu package implementation in R
The use of liureg is explained through examples using the Hald dataset.
> library(liureg) > mod <-liu(y ~X1 + X2 + X3 + X4, data = as.data.frame(Hald), + scaling = "centered", d = seq(0, 1, 0.01) ) The output of linear LR from liu() function is assigned to an object mod.The first argument of the function is formula, which is used to specify the required LR model for the data provided as second argument.The print method for mod, an object of class "liu", will display the de-scaled coefficients.The output (de-scaled coefficients) from the above command is only for a few selected biasing parameter values.

> summary(mod)
The R Journal Vol.9/2, December 2017 ISSN 2073-4859 Call: liu.default(formula = y ~., data = as.data.frame(Hald) The dest() function, which works with Liu fitted models, computes different biasing parameters developed by researchers, see Table 3.The list of different d values (5 in number) may help in deciding the amount of bias needs to be introduced in LR.The biasing parameters by Liu (1993Liu ( , 2011) )  The lstats() function can be used to compute different statistics for a given Liu biasing parameter specified in a call to liu.The Liu statistics are MSE, squared bias, F-statistics, Liu variance, degrees of freedom (df) by Hastie and Tibshirani (1990), and R 2 etc.Following are results using lstats() for some d = −1.47218,−0.06, 0, 0.1, 0.5, 1.
> lstats(mod) The lstats() also displays the value of d which produces minimum MSE among all provided values of d as argument in liu() function.
The residuals, fitted values from the LR, and predicted values of the response variable y can be computed using the fuctions residuals(), fitted(), and predict(), respectively.To obtain the Var-Cov and Hat matrices, the functions vcov() and hatl() can be used.The df are computed by following Hastie and Tibshirani (1990).The results for Var-Cov and diagonal elements of the hat matrix from vcov() and hatl() functions are given below for d = −1.47218.
> vcov(liu(y ~., as.data.frame(Hald) For given values of X, such as for first five rows of X matrix, the predicted values for some d = −1.47218,−0.06, 0, 0.1, 0.5, 1 will be computed by predict(): > predict(mod, newdata = as.data.frame(Hald[1 : 5, The model selection criteria's of AIC and BIC can be computed using infoliu() function for each value of d used in argument of liu().For some d = −1.47218,−0.06, 0.5, 1, the AIC and BIC values are: > infoliu(liu(y ~., as.data.frame(Hald),0.5The effect of multicollinearity on the coefficient estimates can be identified by using different graphical displays such as the Liu trace (see Figure 1); the plotting of bias, variance, and MSE against d (see Figure 2); and plotting the information criteria against d f (Figure 3).These graphical displays are (subjective) methods for selection of the optimal biasing parameter d.

Summary
The liureg package provides the most complete suite of tools for LR available in R, comparable to those available as listed in Table 1.We have implemented functions to compute the Liu coefficients, the testing of these coefficients, the computation of different Liu related statistics and the computation of the biasing parameter for different existing methods by various authors (see Table 3).We have The R Journal Vol.9/2, December 2017 ISSN 2073-4859  The R Journal Vol.9/2, December 2017 ISSN 2073-4859
also provided guidance in coding.

Table 5 :
Components of the "liu" class.Table 6 lists the functions and methods available in liureg package.
The main model fitting function for implementation of LR models in R. coef() Display de-scaled Liu coefficients.liu() Generic function and default method that calls liuest() and returns an object of S3 class "liu" with different set of methods to standard generics.It has a print method for display of Liu de-scaled coefficients.Bias, variance, and MSE plot as a function of d. plot.infoliu()Plot of AIC and BIC against d.

Table 6 :
Functions and methods in liureg package.To obtain Liu scaled coefficients mod$coef can be used: The summary() function also displays Liu related R 2 , adjusted-R 2 , F-test, AIC, BIC, and minimum MSE at certain d given in liu().
include d CL , d mm , d opt , d ILE , and GCV for the appropriate selection of d.
Following are possible uses of some functions to compute different Liu related statistics.For a detailed description of these functions/commands, see the liureg package documentation.