The Comparison of Some Methods in Analysis of Linear Regression Using R Software

Introduction

The processes of model fitting are: model selection, parameters estimation, checking the adequacy of the model, calculating confidence intervals and testing hypotheses about the parameters in the model and interpreting the results. In the chosen model we have to consider: the assumption of observation independence or at least unrelated, a single error term in the model, the choice of scale for the analysis, the choice of variables that will be included in the model.

If the variability comes only from variability of errors and variance of the errors term does not depend from the values of the explanatory variable X, thus have equal variance. This property is called homoscedasticity.

A test that controls the dispersion of errors in a regression is called "Breusch-Pagan test" by Breusch, and Pagan, 1979, which control the hypothesis if the variance of errors is constant versus the alternative the error variance changes with the level of response. With command “ncvTest”, which is into the R program we compute a score test of the hypothesis of constant error variance.

In issue 1.1 is treated least squares method to estimate the parameters of regression, hat matrix, leverages, studentized residuals, variance of parameters and the standard error for a particular component.

In issue 1.2 is treated weight least squares method. In this case the errors have unequal variance where the form is known. A test that controls the dispersion of errors in a regression is called "Breusch-Pagan test" by Breusch, and Pagan, 1979, which control the hypothesis if the variance of errors is constant versus the alternative the error variance changes with the level of response.

In issue 1.3 is treated bootstrap methods to estimate parameters of linear regression. Finally in issue 1.4 is simulated a linear regression with heteroscedasticity errors. Using R software we have found confidence intervals, estimated coefficients, plots and results for any case. In this issue are shown the commands that are used in R to get the results.

1.1 Least Squares Linear Regression

A vector of observations having n component is assumed to be a realization of a random variable Y whose components are independently. If we repeat the observation for a value given by the explanatory variable, we will not necessarily have the same value from the explained variable, so we may have another value because they explained variable will be considered a random variable (Ronald E. Walpole, Raymond H. Myers, Sharon L. Myers, Keying Ye, “Probability and Statistics for engineers and Scientists”).

Explanatory variable X will be considered as fixed values, matrix with size.

The linear regression model is:

(1)

where: , . is covariance matrix, I_n is matrix unit ,

.

To estimate the parameters , we use the least squares method. The sum of the squared errors (SSE) is:

(2)

The least squares estimators, which minimizes are:

(3)

if exists. The equation is called the theoretical regression equation, while the fitted least squares regression equation.

(4)

is called hat matrix. Residuals: .

Studentized residuals are:

(5)

which if the model is assumed correctly they have variance equal to 1 and tend to be small. The quantities are known as leverages, the main diagonal elements of hat matrix H. Studentization of residuals can only correct the natural non-constant variance in residuals when the errors have constant variance. If there is some underlying heteroscedascity in the errors, studentization can not correct for it.

, and

Residual sum of squares:

.

If and , then

, (6)

(7)

The standard error for a particular component is: .

If errors are independent and identically normally distributed with mean 0 and variance, which means, then using the fact that linear combinations of normally distributed values are also normal:

(8)

We can use these results to test hypothesis about parameters and to construct confidence intervals for. The estimators have normal distribution, if the errors have normal distribution. If the errors do not have normal distribution, then this method cannot be safe for accurate conclusions.

1.2 Weighted Least Squares Method

The basic tool for examining the fit is the residuals, and we have already looked for patterns in residuals and assessed the normality of their distribution. Where the errors are uncorrelated, but have unequal variance where the form is known we can use weighted least squares (WLS). A simple case:

(9)

(10)

where for i = 1,…, n are known. In this situation the model is:

, (11)

The sum of the squared errors in this case is:

(12)

Weighted least squares estimations of parameters are:

(13)

Weighted least squares is appropriate when the form of the non-constant variance is either known exactly. A test that controls the dispersion of errors in a regression is called "Breusch-Pagan test" by Breusch, and Pagan, 1979, which control the hypothesis if the variance of errors is constant versus the alternative the error variance changes with the level of response (fitted values), or with a linear combination of predictors.

Suppose we have a regression model:

(14)

(15)

If the function h (.) is linear in relation to independent variables, then the last equation can be written in the form:

We build auxiliary regression: . Estimate the regression coefficients with the least square method and check the hypothesis that all coefficients of the least regression equation, in addition to the first coefficient, get the zero value.

versus

1.3 Re-Sampling Errors and Re-Sampling Cases in Regression

Re-sampling errors are a bootstrap method to provide more accurate analysis. The application of bootstrap in linear regressions was originally made by Bradley Efron and Robert J Tibshirani (1993). The algorithm for model based re-sampling in simple linear regression (A.C.Davison and D.V.Hinkley 1997, p. 262). Another bootstrap method is re-sampling cases (A.C.Davison and D.V.Hinkley 1997, p. 264). In this method we make no assumption if variance is constant.

1.4 Simulation: True Regression, OLS Regression, WLS Regression and Bootstrap Methods. Using R Software

We start with regression where we know the true linear regression and simulation the errors of regression:

where Y- define the response variable. - predictor variable (explanatory variable).

We simulate the errors in regression such that their variances are linearly dependent on variable. The aim here is to see the difference between results derived from different methods. At first download packages in R program: library(car), library(MASS), library(nlme), library(boot).

The blue and bold words written below are the commands we used in the program R.

x<-seq(0,29, by=1) # explanatory variable

y = 20+10*x + rnorm(30,0, sapply(x, function(x){1+5*x})) # response variable

Sim.Data<-Data.Frame(Y, X)

OLS Regression

fit.ols = lm(y~x, data=sim.data) # OLS regression

Summary(Fit.Ols)

Residual standard error: 74.56 on 28 degrees of freedom. Multiple R-squared: 0.6028, Adjusted R-squared: 0.5886. F-statistic: 42.48 on 1 and 28 DF, p-value: 0.0001305.

We find the confidence intervals for the coefficients of the OLS regression using: confint(fit.ols).

WLS Regression

With command ncvTest computes a score test of the hypothesis of constant error variance against the alternative that the error variance changes with the level of the response (fitted values), or with a linear combination of predictors. Breusch, T. S. and Pagan, A. R.1979.

ncvTest(fit.ols) # test for heteroscedasticity

Non-constant Variance Score Test. Chisquare = 5.75913, Df = 1, p = 0.01640319. The p-values suggestion to reject the hypothesis of homoscedasticity and we accept heteroscedasticity random errors. We expected this result because such as we designed simulated errors. Now we use weight least square method. Try fitting the regression with weights and see what the difference is.

Fit.Wls = Lm(Y~X, Weights=1/(1+5*X), Data=Sim.Data)

Summary(Fit.Wls)

Residual standard error: 7.674 on 28 degrees of freedom. Multiple R-squared: 0.7805. Adjusted R-squared: 0.7726. F-statistic: 99.55 on 1 and 28 DF, p-value: 1.012e-10. We find the confidence intervals for the coefficients of the WLS regression using: confint (fit.wls).

We use the following commands to make the graphics of the regressions:

Name <-C("True Regression","OLS Regression", "WLS Regression")

Colors.Reg <- C("Red", "Blue", "Black")

Type.Line <- C(1, 2, 4)

Line.Width <- C(3, 2, 4)

plot(x, y, main = "True regression, OLS regression and WLS regression", lwd=1, bty="o")

Abline(a=20, B=10, Col=Colors.Reg[1], Lty = Type.Line[1], Lwd =Line.Width [1])

abline(fit.ols$coefficients, col=colors.reg[2], lty = type.line[2], lwd =line.width[2])

abline(fit.wls$coefficients, col=colors.reg[3], lty = type.line[3], lwd =line.width[3])

legend ("topleft", legend = name, col = colors.reg, lty = type.line, lwd = line.width, bty ="n")

Figure 1 Truth regression, OLS regression and WLS regression for the simulate data.

Bootstrap method: re-sampling errors.

Another method for to estimating regression coefficients and standard deviation is the bootstrap method of re-sampling errors.

Sim.Lm<-Glm(Y~X, Data=Sim.Data) # the Regression

Sim.Fit<-Function(Data) Coef(Glm(Data$Y~Data$X))

sim.diag<-glm.diag.plots(sim.lm, ret=T) # Diagnostics plots

Sim.Res<-Sim.Diag$Res*Sim.Diag$Sd

Sim.Res<-Sim.Res-Mean(Sim.Res)

sim.df<-data.frame(sim.data, res=sim.res, fit=fitted(sim.lm))

Sim.Model<-Function(Data, I){

D<-Data

D$Y<-D$Fit+D$Res[I]

Sim.Fit(D)

}

Fit.Boot<-Boot(Sim.Df, Sim.Model, R=9999)

Figure 2 Diagnostics plots for sim.lm regression.

The plot on the top left is a plot of the jackknife deviance residuals against the fitted values. The plot on the top right is a normal QQ plot of the standardized deviance residuals. The dotted line is the expected line if the standardized residuals are normally distributed. The bottom two panels are plots of the Cook statistics. On the left is a plot of the Cook statistics against the standardized leverages. The final plot again shows the Cook statistic this time plotted against case number enabling us to find which observations are influential (Davison, A.C. and Snell, E.J., 1991).

Summary (Fit.Boot)

Bootstrap percentile interval and adjusted bootstrap percentile (BCa) interval.

Calculations based on 9999 bootstrap replicates.

boot.ci(fit.boot, index = 1:min(2, length(fit.boot$t0)), type = c('perc', 'bca'))

Level Percentile BCa

95% (-41.597, 62.66) (-39.150, 65.582)

boot.ci(fit.boot, index = 2:min(2, length(fit.boot$t0)), type = c('perc','bca'))

Level Percentile BCa

95% (6.98, 13.08) (7.00, 13.11)

Bootstrap method: re-sampling case.

Sim.Fit<-Function(Data)Coef(Glm(Data$Y~Data$X))

Sim.Case<-Function(Data, I)Sim.Fit(Data[I,])

Fit.Boot.Case<-Boot(Sim.Data, Sim.Case, R=9999)

Fit.Boot.Case

Bootstrap method: weight error resamling.

Plot(Fit.Boot.Case, Index=2, Jack=T)

Sim.Lm.W<-Glm(Y~X,Weight=1/(1+5*X),Data=Sim.Data)

Sim.Fit.W<-Function(Data) Coef(Glm(Data$Y~Data$X))

Sim.Diag.W<-Glm.Diag.Plots(Sim.Lm.W, Ret=T)

The last command makes plot of jackknife deviance residuals against linear predictor, normal scores plots of standardized deviance residuals, plot of approximate Cook statistics against where are the leverages.

sim.res.w<-sim.res-mean(sim.res.w)

sim.df.w<-data.frame(sim.data, res=sim.res.w, fit=fitted(sim.lm.w))

Sim.Model<-Function(Data, I){ D<-Data

D$Y<-D$Fit+D$Res[I]

Sim.Fit(D)}

fit.boot.w<-boot(sim.df.w, sim.model, R=9999); fit.boot.w

We find bootstrap percentile interval and adjusted bootstrap percentile (BCa) interval only for the second coefficient of the regression fit.boot.w.

Level Percentile BCa

95% (6.318, 12.342) (6.471, 12.518)

Bootstrap Methods: Weight Errors and Re-Sampling Case

sim.fit.w<-function(data)coef(glm(data$y~data$x,weight=1/(1+5*x)))

Sim.Case.W<-Function(Data, I)Sim.Fit(Data[I,])

fit.boot.case.w<-boot(sim.data, sim.case, R=9999); summary(fit.boot.case.w)

Summary table: Estimate errors and standard deviations.

Method	Coefficients	Std. Error
True regression	=20	-
	=10	-
OLS regression	9.348	25.969
	10.024	1.538
WLS regression	19.5242	6.5710
	9.3218	0.9343
Bootstrap method resampling errors without weight	9.038062	26.3335
	10.045030	1.5606
Bootstrap method resampling cases without weight	9.3481021	18.40179
	10.023602	1.64067
Bootstrap method resampling errors with weights	19.289300	26.0353
	8.957616	1.5364
Bootstrap method resampling case with weights	9.348102	18.5540
	10.02364	1.6552

Conclusions

For linear regression with normal random errors having constant variance, the least squares method of the coefficients estimators and standard deviation are accurate, but when the random errors having non constant variance, the estimators with this method are not accurate. The bootstrap methods improved quantile estimation. The WLS method is the most accurate where we know weights errors. Where we don’t know completely the form of the variance of error, we can use bootstrap method resampling errors. We can use packages inside R software to get results.

References