3.11 Regression with large number of regressors
In Practice, we often encounter data with a large number of regressors. However, the problem with that kind of data is we do not know what variables should include and what should not. That is because some of the non-significant variables might be related to other significant ones. Another critical issue is related to the difference between inferential models predictive models. In the former case, we often prefer unbiasedness over variability so that the resulted models will be more complex. Whereas, with predictive models, we prefer less complex ones to get less variability, hence more accurate prediction.
That is why researchers have proposed some strategies for the trade-off between these two goals. We will demonstrate those strategies through an example with simulated data. But, for pedagogical purposes, we will use a small number of regressors.
The regressors can also be the result of some transformations or combinations of the original ones such as: \(x^2\), \(2x_1-3x_2\), \(log(x)\), \(\sqrt x\), etc.
Since those strategies are based on some well-known criteria to choose between models, we will introduce the most popular ones.
The coefficient of determination \(R^2\):
As shown in (3.59) or in (3.60).
The model with higher value for \(R^2\) will be preferred.
The adjusted coefficient of determination \(\overline R^2\):
Since the previous \(R^2\) always prefers larger models, This one, shown in (3.61), tends to penalize larger models against smaller ones.
The mean squared error \(MSE\):
It is the sum squared errors \(SSR\), defined earlier, divided by the number of observations.
The Akaike information criterion \(AIC\):
It is given by:
\[\begin{equation} AIC=-2l_{(\beta,\sigma^2)}^{max}+2(k+1) \tag{3.91} \end{equation}\]
Where \(l_{(\beta,\sigma^2)}^{max}\) is the maximum likelihood discussed in the section 3.3, and \(K\) is the numder of parameters (+1 for the variance). Since this criterion (as well as the following ones) does not provide a measure for goodness of fit, it is used only to rank different models. That is why, that formula will be reduced to the following (after dsicarding constant terms:
\[\begin{equation*} \begin{split} AIC&=-2\Bigg(-\frac{n}{2}log(2\pi)-\frac{n}{2}log(s^2)-\frac{1}{2s^2}\overbrace{(y-X\widehat\beta )^t(y-X\widehat\beta )}^{=(n-k)s^2}\Bigg)+2(k+1)\\ &=-2\Bigg(\overbrace{-\frac{n}{2}log(2\pi)}^{constant}-\frac{n}{2}log(s^2)\overbrace{-\frac{(n-k)}{2}}^{constant}\Bigg)+2(k+1)\\ &=nlog\Bigg(\frac{SSR}{n-k}\Bigg)+2(k+1) \end{split} \end{equation*}\]
Notice that we have used the unbiased standard error \(s^2\) (divided by \((n-k)\)) and not \(\sigma_{Ml}\) (divided by \(n\)). However, that matters only in small samples. Otherwise, they converge to each other.
The model with the smallest \(AIC\) value will be preferred, such that each time we decrease \(k\) (the number of parameters), the \(AIC\) value decreases, and hence the model will get a higher rank of preference.
There exist a corrected version for Akaike used in small samples given by:
\(AIC_{corrected}=AIC+\frac{2(k+1)(k+2)}{n-k-2}\)The Bayes information criterion \(BIC\):
The same as the previous one except for the penalty term. It is given by:
\[\begin{equation} BIC=nlog(\frac{SSR}{n})+(k+1)log(k) \tag{3.92} \end{equation}\]
We can compare these two criteria by only comparing their penalty terms as follows:
\[\begin{equation*} \begin{split} BIC>AIC&\implies(k+1)log(k)>2(k+1)\\ &\implies log(k)>2\\ &\implies k>e^2\\ &\implies k>7 \end{split} \end{equation*}\]
That means that \(BIC\) is parsimonious than \(AIC\) if the number of the model parameters larger than 7. That is why in predictive models, practitioners prefer \(BIC\) since models with less complexity tend to give accurate predictions (because they have less variability, as discussed earlier).
The Hannon and Quin \(HQ\):
This criterion is also base don the maximum likelihood with a different penalty term than the previous ones. it is given by:
\[\begin{equation} HQ=nlog\Bigg(\frac{SSR}{n}\Bigg)+2c\Big(k+1\Big)log\Big(log(n)\Big) \tag{3.93} \end{equation}\]
Where \(c\) is any constant greater than one.
This test has the smallest penalty term so that it tends to prefer more complex models.
The Bridge criterion \(BC\):
This criterion tries to obtain the advantages of both, \(AIC\) and \(BIC\), and used only over models, where the number of parameters of each is less than that of the selected model by \(AIC\). Because \(AIC\) always selects models with larger number of parameter, this number, thus, will be used as an upper bound for the candidate models in \(BC\). It is given by:
\[\begin{equation} BC=nlog\Bigg(\frac{SSR}{n}\Bigg)+c\Big(1+2^{-1}+..+(k+1)^{-1}\Big) \tag{3.94} \end{equation}\]
Where the suggested value for \(c\) is \(n^{2/3}\).
The root mean squared error \(RMSE\)
In machine learning, The way to evaluate the performance of a model is by holding out a small set of observations from the current sample to compare the predicted values to the actual ones using some appropriate metrics. In econometrics, if the sample size is large enough so that we can hold out some observation, we can use the \(RMSE\) test, which works precisely as \(MSE\), but applies only to the hold-out samples.
say we have used only \((n-10)\) for fitting our model, then \(RMSE\) criterion can be used to assess the model performance by making use the 10 hold-out samples as follows:
\[\begin{equation} RMSE=\sqrt{\frac{1}{10}\sum\limits^{n}_{t=n-10}(y_t-\widehat y_t)^2} \tag{3.95} \end{equation}\]
For more depth into those criteria, read this paper written by(Ding et al. 2018) in the following link.
3.11.1 All possible regression
This option is the most expensive choice, so that if we have \(k\) potential regressors, then we should fit \(2^k\) models. For instance, with 3 regressors \(x_2\), \(x_3\), and \(x_4\), we should fit 8 models: (\(x_2\)), (\(x_3\)), (\(x_4\)), (\(x_2,x_3\)), (\(x_2,x_4\)), (\(x_3,x_4\)), (\(x_2,x_3,x_4\)), and the null model with only the constant term.
To implement this method in R, we make use of the R package olsrr, but first, let us generate some data.
# generate the relevant regressors
df_true <- tibble(matrix(nrow = 1000, ncol=5))
for (i in seq(2,4)) {
set.seed(i)
df_true[,i] <- rnorm(1000, mean=i*5, sd=1)
names(df_true)[i] <- paste0('x',i)
}
# generate an irrelevant regressor
set.seed(5)
df_true$x5 <- 20*rbeta(1000, 6, 8)
# generate the error term
set.seed(6)
err <-rnorm(1000, 0, 5)
# specify the true model
y <- 1.5+2.9*df_true$x2+2.3*df_true$x3+2.1*df_true$x4+err
df_true[,1] <- y
names(df_true)[1] <- 'y'
# take out a random sample of 70 observations
set.seed(11)
df <- df_true %>% sample_n(70)
head(df)| y | x2 | x3 | x4 | x5 |
|---|---|---|---|---|
| 100.57 | 10.34 | 13.97 | 19.55 | 12.59 |
| 99.46 | 9.72 | 15.74 | 18.52 | 8.35 |
| 106.64 | 10.05 | 16.41 | 20.68 | 9.21 |
| 116.71 | 10.65 | 15.64 | 20.37 | 9.96 |
| 114.82 | 10.09 | 14.85 | 20.62 | 7.07 |
| 101.16 | 10.47 | 14.59 | 20.87 | 3.07 |
We expect then the variable \(x_5\) to be non significant.
Now we can fit all the possible regressions by calling the function ols_step_all_possible that displays all the models along with a bunch of criteria.
library(olsrr)
mod_R_all <- lm(y~., data=df)
output_R_all <- ols_step_all_possible(mod_R_all)
output_R_all| n | predictors | rsquare | adjr | aic | sbc | |
|---|---|---|---|---|---|---|
| 3 | 1 | x4 | 0.1290 | 0.1162 | 452.6337 | 459.3791 |
| 2 | 1 | x3 | 0.1212 | 0.1083 | 453.2538 | 459.9993 |
| 1 | 1 | x2 | 0.0563 | 0.0424 | 458.2427 | 464.9882 |
| 4 | 1 | x5 | 0.0030 | -0.0116 | 462.0872 | 468.8327 |
| 8 | 2 | x3 x4 | 0.2361 | 0.2133 | 445.4498 | 454.4438 |
| 6 | 2 | x2 x4 | 0.2157 | 0.1923 | 447.2892 | 456.2832 |
| 5 | 2 | x2 x3 | 0.1728 | 0.1481 | 451.0195 | 460.0135 |
| 10 | 2 | x4 x5 | 0.1300 | 0.1041 | 454.5490 | 463.5430 |
| 9 | 2 | x3 x5 | 0.1216 | 0.0954 | 455.2243 | 464.2183 |
| 7 | 2 | x2 x5 | 0.0620 | 0.0340 | 459.8183 | 468.8122 |
| 11 | 3 | x2 x3 x4 | 0.3156 | 0.2845 | 439.7570 | 450.9995 |
| 14 | 3 | x3 x4 x5 | 0.2361 | 0.2014 | 447.4498 | 458.6923 |
| 13 | 3 | x2 x4 x5 | 0.2188 | 0.1833 | 449.0181 | 460.2605 |
| 12 | 3 | x2 x3 x5 | 0.1744 | 0.1368 | 452.8877 | 464.1302 |
| 15 | 4 | x2 x3 x4 x5 | 0.3161 | 0.2740 | 441.7051 | 455.1961 |
Now if we want to retrieve the best subset of regressors using one of the above criteria, say \(AIC\), then we should call the ols_step_best_subset function as follows:
| n | predictors | rsquare | adjr | aic | sbc | |
|---|---|---|---|---|---|---|
| 3 | 3 | x2 x3 x4 | 0.3156 | 0.2845 | 439.7570 | 450.9995 |
| 4 | 4 | x2 x3 x4 x5 | 0.3161 | 0.2740 | 441.7051 | 455.1961 |
| 2 | 2 | x3 x4 | 0.2361 | 0.2133 | 445.4498 | 454.4438 |
| 1 | 1 | x4 | 0.1290 | 0.1162 | 452.6337 | 459.3791 |
The true model (number 3) has been recognized by all The criteria. The smallest value for \(AIC\) is \(439.7570\), and for \(BIC\) is \(450.9995\).
Unfortunately, statsmodels Python package does not have a selection method. However, we can create a function to fit all possible regression as follows. But first, let us retrieve the data into python.
Then we use the itertools package to generate all the possible subsets from the set of the available regressors.
from itertools import chain, combinations
# create a function to display all the subsets
def all_subset(arr):
s = list(arr)
return chain.from_iterable(combinations(s,r) for r in range(1,len(s)+1))
subsets = list(all_subset(df_py.columns[1:,].values))Once we have got all the possible subsets, we can loop over each subset to fit the associated model.
import statsmodels.formula.api as smf
import statsmodels.api as sm
predictors = []
rsquare = []
adjr = []
aic = []
sbc = []
for i in range(len(subsets)):
data = pd.concat([df_py[['y']],df_py[list(subsets[i])]], axis=1)
mod_py_all = smf.ols('y~'+'+'.join(list(subsets[i])), data=data).fit()
predictors.append(subsets[i])
rsquare.append(mod_py_all.rsquared)
adjr.append(mod_py_all.rsquared_adj)
aic.append(mod_py_all.aic)
sbc.append(mod_py_all.bic)
# put the results in a data frame
out_py = pd.DataFrame({'predictors':predictors, 'rsquare': rsquare, 'adjr':adjr, 'aic':aic,
'sbc':sbc})
# order the results by AIC criterion and display the 5 best models
out = out_py.sort_values(by='aic').head(5)
out| predictors | rsquare | adjr | aic | sbc | |
|---|---|---|---|---|---|
| 10 | x2, x3, x4 | 0.3156 | 0.2845 | 437.7570 | 446.7510 |
| 14 | x2, x3, x4, x5 | 0.3161 | 0.2740 | 439.7051 | 450.9476 |
| 7 | x3, x4 | 0.2361 | 0.2133 | 443.4498 | 450.1953 |
| 5 | x2, x4 | 0.2157 | 0.1923 | 445.2892 | 452.0347 |
| 13 | x3, x4, x5 | 0.2361 | 0.2014 | 445.4498 | 454.4438 |
As expected, python has recognized the best model.
3.11.2 compare two specific models
If we want to compare two specific nested models we can use ANOVA table.
suppose we want to compare the full model with the right model (without \(x_5\))
| term | df.residual | rss | df | sumsq | statistic | p.value |
|---|---|---|---|---|---|---|
| y ~ x2 + x3 + x4 | 66 | 1900.750 | NA | NA | NA | NA |
| y ~ x2 + x3 + x4 + x5 | 65 | 1899.342 | 1 | 1.4077 | 0.0482 | 0.827 |
Since the p-value (\(0.8269568\)) is far larger than the most used threshold (\(0.05\)), we can not reject the null hypothesis that says the two models are similar, and hence we keep the simpler model (without \(x_5\)), which is, in our case, a correct decision.
mod_py1 = smf.ols('y~x2+x3+x4', data=df_py).fit()
mod_py2 = smf.ols('y~x2+x3+x4+x5', data=df_py).fit()
out_py1 = sm.stats.anova_lm(mod_py1,mod_py2)
out_py1| df_resid | ssr | df_diff | ss_diff | F | Pr(>F) |
|---|---|---|---|---|---|
| 66 | 1900.750 | 0 | NaN | NaN | NaN |
| 65 | 1899.342 | 1 | 1.4077 | 0.0482 | 0.827 |
3.11.3 Forward selection
In this method, we start with the null model (intercept only) and add one regressor at each time. The regressor added should be the one that has the lowest p-value resulted from the regression of the response on that regressor. Then, we continue that process until all the remaining regressors are non-significant (depending on the level used, say 0.05 ). However, we can use another criterion in place of the p-value, such as \(AIC\) or \(BIC\), ..etc.
To get the best model using the above method, We will call the olsrr package function ols_step_forward_aic() that uses \(AIC\), or if we want to instead use p-value then we should call the function ols_step_forward_p().
mod_R_forward <- lm(y~., data=df)
# set details to TRUE to display all intermediate steps
output_R_forward <- ols_step_forward_aic(mod_R_forward, details = TRUE)[out] Forward Selection Method
[out] ------------------------
[out]
[out] Candidate Terms:
[out]
[out] 1. x2
[out] 2. x3
[out] 3. x4
[out] 4. x5
[out]
[out]
[out] Step => 0
[out] Model => y ~ 1
[out] AIC => 460.3004
[out]
[out] Initiating stepwise selection...
[out]
[out] Table: Adding New Variables
[out] -----------------------------------------------------------------------
[out] Predictor DF AIC SBC SBIC R2 Adj. R2
[out] -----------------------------------------------------------------------
[out] x4 1 452.634 459.379 253.317 0.12899 0.11618
[out] x3 1 453.254 459.999 253.902 0.12124 0.10831
[out] x2 1 458.243 464.988 258.617 0.05632 0.04244
[out] x5 1 462.087 468.833 262.254 0.00304 -0.01162
[out] -----------------------------------------------------------------------
[out]
[out] Step => 1
[out] Added => x4
[out] Model => y ~ x4
[out] AIC => 452.6337
[out]
[out] Table: Adding New Variables
[out] ----------------------------------------------------------------------
[out] Predictor DF AIC SBC SBIC R2 Adj. R2
[out] ----------------------------------------------------------------------
[out] x3 1 445.450 454.444 246.581 0.23608 0.21328
[out] x2 1 447.289 456.283 248.265 0.21574 0.19233
[out] x5 1 454.549 463.543 254.930 0.13004 0.10407
[out] ----------------------------------------------------------------------
[out]
[out] Step => 2
[out] Added => x3
[out] Model => y ~ x4 + x3
[out] AIC => 445.4498
[out]
[out] Table: Adding New Variables
[out] ----------------------------------------------------------------------
[out] Predictor DF AIC SBC SBIC R2 Adj. R2
[out] ----------------------------------------------------------------------
[out] x2 1 439.757 450.999 241.703 0.31559 0.28448
[out] x5 1 447.450 458.692 248.509 0.23608 0.20136
[out] ----------------------------------------------------------------------
[out]
[out] Step => 3
[out] Added => x2
[out] Model => y ~ x4 + x3 + x2
[out] AIC => 439.757
[out]
[out] Table: Adding New Variables
[out] ----------------------------------------------------------------------
[out] Predictor DF AIC SBC SBIC R2 Adj. R2
[out] ----------------------------------------------------------------------
[out] x5 1 441.705 455.196 243.811 0.31609 0.27401
[out] ----------------------------------------------------------------------
[out]
[out]
[out] No more variables to be added.
[out]
[out] Variables Selected:
[out]
[out] => x4
[out] => x3
[out] => x2[out]
[out]
[out] Stepwise Summary
[out] -------------------------------------------------------------------------
[out] Step Variable AIC SBC SBIC R2 Adj. R2
[out] -------------------------------------------------------------------------
[out] 0 Base Model 460.300 464.797 260.983 0.00000 0.00000
[out] 1 x4 452.634 459.379 253.317 0.12899 0.11618
[out] 2 x3 445.450 454.444 246.581 0.23608 0.21328
[out] 3 x2 439.757 450.999 241.703 0.31559 0.28448
[out] -------------------------------------------------------------------------
[out]
[out] Final Model Output
[out] ------------------
[out]
[out] Model Summary
[out] ---------------------------------------------------------------
[out] R 0.562 RMSE 5.211
[out] R-Squared 0.316 MSE 27.154
[out] Adj. R-Squared 0.284 Coef. Var 5.015
[out] Pred R-Squared 0.227 AIC 439.757
[out] MAE 4.078 SBC 450.999
[out] ---------------------------------------------------------------
[out] RMSE: Root Mean Square Error
[out] MSE: Mean Square Error
[out] MAE: Mean Absolute Error
[out] AIC: Akaike Information Criteria
[out] SBC: Schwarz Bayesian Criteria
[out]
[out] ANOVA
[out] --------------------------------------------------------------------
[out] Sum of
[out] Squares DF Mean Square F Sig.
[out] --------------------------------------------------------------------
[out] Regression 876.442 3 292.147 10.144 0.0000
[out] Residual 1900.750 66 28.799
[out] Total 2777.193 69
[out] --------------------------------------------------------------------
[out]
[out] Parameter Estimates
[out] ---------------------------------------------------------------------------------------
[out] model Beta Std. Error Std. Beta t Sig lower upper
[out] ---------------------------------------------------------------------------------------
[out] (Intercept) 5.418 18.632 0.291 0.772 -31.782 42.618
[out] x4 2.612 0.704 0.383 3.710 0.000 1.206 4.017
[out] x3 2.045 0.659 0.317 3.103 0.003 0.729 3.361
[out] x2 1.852 0.669 0.285 2.769 0.007 0.516 3.187
[out] ---------------------------------------------------------------------------------------As we can see, in the first step, we regressed the response on every single regressor separately then selected the one with the smallest \(AIC\) value, which was \(x_4\) (\(AIC= 452.634\)). Next, we regressed \(y\) on \(x_4\) and every regressor, then selected the best one in terms of \(AIC\), which was \(x_3\) (\(AIC=445.450\)). In the following step \(x_2\) was selected (\(AIC=439.757\)). Finally, when we add the last regressor \(x_5\), the \(AIC\) of this model was \(AIC=441.705\), which is larger than the \(AIC\) of the previous model (\(AIC=439.757\)), that is why we have not included this last regressor. Again, we have selected the right model.
If we want to display only the final model, we should set the argument details=FALSE (which is the default)
To conduct the forward selection method, we will use the following function found here, but it uses the p-value.
import statsmodels.formula.api as smf
def forward(data, y):
init_reg = data.columns[1:].tolist()
best_reg = []
while (len(init_reg)>0):
remaining_reg = list(set(init_reg)-set(best_reg))
new_p = pd.Series(index=remaining_reg, dtype='float64')
for new_column in remaining_reg:
model = smf.ols('y~'+ '+'.join(init_reg+[new_column]), data=data).fit()
new_p[new_column] = model.pvalues[new_column]
min_p = new_p.min()
if(min_p<0.1):
best_reg.append(new_p.idxmin())
else:
break
return best_reg
This function uses p-value to select variables.
[out] ['x4', 'x3', 'x2']3.11.4 Backawrd elimination
This method is the inverse of the previous one, so that we start with the full model, then remove the regressor that has the highest p-value (if we use the p-value) at each step. Then the process continues until we left with only significant regressors.
[out]
[out]
[out] Stepwise Summary
[out] -------------------------------------------------------------------------
[out] Step Variable AIC SBC SBIC R2 Adj. R2
[out] -------------------------------------------------------------------------
[out] 0 Full Model 441.705 455.196 243.811 0.31609 0.27401
[out] 1 x5 439.757 450.999 241.703 0.31559 0.28448
[out] -------------------------------------------------------------------------
[out]
[out] Final Model Output
[out] ------------------
[out]
[out] Model Summary
[out] ---------------------------------------------------------------
[out] R 0.562 RMSE 5.211
[out] R-Squared 0.316 MSE 27.154
[out] Adj. R-Squared 0.284 Coef. Var 5.015
[out] Pred R-Squared 0.227 AIC 439.757
[out] MAE 4.078 SBC 450.999
[out] ---------------------------------------------------------------
[out] RMSE: Root Mean Square Error
[out] MSE: Mean Square Error
[out] MAE: Mean Absolute Error
[out] AIC: Akaike Information Criteria
[out] SBC: Schwarz Bayesian Criteria
[out]
[out] ANOVA
[out] --------------------------------------------------------------------
[out] Sum of
[out] Squares DF Mean Square F Sig.
[out] --------------------------------------------------------------------
[out] Regression 876.442 3 292.147 10.144 0.0000
[out] Residual 1900.750 66 28.799
[out] Total 2777.193 69
[out] --------------------------------------------------------------------
[out]
[out] Parameter Estimates
[out] ---------------------------------------------------------------------------------------
[out] model Beta Std. Error Std. Beta t Sig lower upper
[out] ---------------------------------------------------------------------------------------
[out] (Intercept) 5.418 18.632 0.291 0.772 -31.782 42.618
[out] x2 1.852 0.669 0.285 2.769 0.007 0.516 3.187
[out] x3 2.045 0.659 0.317 3.103 0.003 0.729 3.361
[out] x4 2.612 0.704 0.383 3.710 0.000 1.206 4.017
[out] ---------------------------------------------------------------------------------------As shown in the table, the variable \(x5\) has been removed, so the final model is the right model.
As we did earlier, we will create a function to enable this method.
import statsmodels.formula.api as smf
def backward(data, y):
reg = data.columns[1:].tolist()
while (len(reg)>0):
model = smf.ols('y~'+ '+'.join(reg), data=data).fit()
p = model.pvalues[1:]
max_p = p.max()
if(max_p>=0.05):
rm_reg = p.idxmax()
reg.remove(rm_reg)
else:
break
return reg
[out] ['x2', 'x3', 'x4']3.11.5 Stepwise selection
This method combines both the forward and backward, such that at each step some variables could be included in the model and some others could be removed until all the regressors are significant.
Unlike the previous R examples, This time we will use the \(p-value\).
[out]
[out]
[out] Stepwise Summary
[out] -------------------------------------------------------------------------
[out] Step Variable AIC SBC SBIC R2 Adj. R2
[out] -------------------------------------------------------------------------
[out] 0 Base Model 460.300 464.797 260.983 0.00000 0.00000
[out] 1 x4 (+) 452.634 459.379 253.317 0.12899 0.11618
[out] 2 x3 (+) 445.450 454.444 246.581 0.23608 0.21328
[out] 3 x2 (+) 439.757 450.999 241.703 0.31559 0.28448
[out] -------------------------------------------------------------------------
[out]
[out] Final Model Output
[out] ------------------
[out]
[out] Model Summary
[out] ---------------------------------------------------------------
[out] R 0.562 RMSE 5.211
[out] R-Squared 0.316 MSE 27.154
[out] Adj. R-Squared 0.284 Coef. Var 5.015
[out] Pred R-Squared 0.227 AIC 439.757
[out] MAE 4.078 SBC 450.999
[out] ---------------------------------------------------------------
[out] RMSE: Root Mean Square Error
[out] MSE: Mean Square Error
[out] MAE: Mean Absolute Error
[out] AIC: Akaike Information Criteria
[out] SBC: Schwarz Bayesian Criteria
[out]
[out] ANOVA
[out] --------------------------------------------------------------------
[out] Sum of
[out] Squares DF Mean Square F Sig.
[out] --------------------------------------------------------------------
[out] Regression 876.442 3 292.147 10.144 0.0000
[out] Residual 1900.750 66 28.799
[out] Total 2777.193 69
[out] --------------------------------------------------------------------
[out]
[out] Parameter Estimates
[out] ---------------------------------------------------------------------------------------
[out] model Beta Std. Error Std. Beta t Sig lower upper
[out] ---------------------------------------------------------------------------------------
[out] (Intercept) 5.418 18.632 0.291 0.772 -31.782 42.618
[out] x4 2.612 0.704 0.383 3.710 0.000 1.206 4.017
[out] x3 2.045 0.659 0.317 3.103 0.003 0.729 3.361
[out] x2 1.852 0.669 0.285 2.769 0.007 0.516 3.187
[out] ---------------------------------------------------------------------------------------After three steps, we have obtained the right model.
To implement stepwise regression, we will use the following function (using p-value), found here.
def stepwise(data, target):
init_reg = data.columns.tolist()
best_reg = []
while init_reg:
remaining_reg = list(set(init_reg)-set(best_reg))
new_p = pd.Series(index=remaining_reg, dtype='float64')
for new_column in remaining_reg:
model = sm.OLS(target, sm.add_constant(data[best_reg+[new_column]])).fit()
new_p[new_column] = model.pvalues[new_column]
min_p_value = new_p.min()
if(min_p_value<0.05):
best_reg.append(new_p.idxmin())
while(len(best_reg)>0):
best_reg_with_constant = sm.add_constant(data[best_reg])
p_values = sm.OLS(target, best_reg_with_constant).fit().pvalues[1:]
max_p_value = p_values.max()
if(max_p_value >= 0.05):
excluded_reg = p_values.idxmax()
best_reg.remove(excluded_reg)
else:
break
else:
break
return best_reg
# run the stepwise regression
stepwise(df_py.iloc[:,1:], df_py.y)[out] ['x4', 'x3', 'x2']