I want to perform an Ordinary Least Squares (OLS) regression (global) or even a Geographically Weighted (GWR) regression (local) to identify possible causes for specific phenomena, like crime.

I have around 60 independent variable candidates, that might be explanatory for the dependent variable.

To identify which of these variable candidates to use for in an OLS or GWR model, I use the tool "Exploratory Regression" within ArcGIS. The tool allows to Input multiple variable candidates and model criteria, tries all different combinations and lists the results, which can then be compared based on the AICc, R2 adjusted, VIF, etc.

I am aware how to use the tool, how to interpret the result and how the overall procedure works.

What I am unsure is, how to limit my huge number of variable candidates? The calculations are either too time consuming or not commutable at all.

One approach, suggested on the ESRI Website is to start with a low minimum and maximum number of candidates to be tested, e.g. starting with 1, then 2, and so on. When I run it with 1 candidate to be tested, I get a significance of 100% or 0% for each variable candidate. If it is done for 2, I get many different percentage values. If those lower e.g. 50% are removed, should I run the model again for 2 tested candidates or already for 3? There appear to be different results. Is there any threshold that could be used?

  • 4
    This is a case of model development or model fitting or variable selection for a linear model. There are many threads about the general problem on Cross Validated. One common theme there is that the ESRI suggestion is not a good approach. For OLS, regularization methods like the LASSO and Ridge Regression, coupled with cross-validation, are readily available and have been extremely well implemented in the R package glmnet. Some active members on Cross Validated, such as Frank Harrell, are published authors and well-known experts on this subject. – whuber Mar 3 '16 at 17:19
  • 2
    This is not a spatial regression! An OLS regression is more robust to iid assumptions so, is preferential for exploratory analysis where spatial autocorrelation may be present in the data. A spatial regression can take many forms but has a term that accounts for space in some way. A GWR would "more-or-less" be considered spatial but suffers from some serious issues and should never be used outside of an exploratory analysis. @whuber has a very good suggestion in specifying a LASSO regression with a cross-validation for parameter selection. This could even be a semi-parametric model. – Jeffrey Evans Mar 3 '16 at 19:15

It seems you are doing research. Here are some advices:

  1. Ask the domain experts: if you ask from an expert about the exploratory variables, he/she can easily filter many of the variables due to their irrelevancy to the response. If you don't have access to an expert, then review the literature.
  2. Use PCA (principal component analysis): PCA is a type of analysis that looks for correlation in between your exploratory variables. It then extracts few (2-3 or more) number of variables that account for the variance of all exploratory random variables. It also characterizes the amount of contribution of variables to a specific principal component. Note that in ArcGIS PCA is part of spatial analyst and can be leveraged with only raster datasets. However you can turn to other statistical software such as R or even python.

Best practice in machine learning and statistical analysis suggest to always avoid PCA in the first run of the inferential analysis. Just with the case of disk space and memory utilization optimization, PCA is recommended. Here is the design of a learning system with/without PCA:

  • Get training set {(X1,Y1),(X2,Y2)...}
  • Run PCA to reduce Xi in dimension to get Zi (optional step)
  • Train the (spatial) regression on {(Z1, Y1),(Z2, Y2),...}
  • Test on test set: Map Xi_test to Zi_test. Run the hypothesis function h(x) on {(Z1_test, Y),(Z2_test, Y2),....}

So here we do not inferr based on principal components (Zi s). We are just using PCA as an optimization technique to gain performance.

  • 1
    @whuber and I had a debate sometime ago regarding PCA. I both acquiesce and stand my ground on this issue. Whereas, I agree that in a pure data reduction context, PCA is perfectly valid. However, where I think that PCA can be dangerous is in inferential regression. If you have nice linear data in the x matrix then things are rosy but, if that x space is nonlinear then PCA can give very erroneous results because the hyperplane that defines the PCA space is effectively a first-order polynomial. If the multivariate space is curvilinear then you have serious issues in any inference. – Jeffrey Evans Apr 12 '16 at 15:37
  • @JeffreyEvans, PCA is most effective if a linear relation exists between predictors. Otherwise PCA will end up with the case of C (# components) = P (predictors). PCA is a lossy data compression technique and it depends on how much of the data variance is to be preserved e.g. 99%. In the worst case scenario we will have to keep all of the predictors. That being said, PCA doesn't impede us from using it in inferential regression rather it may not be that useful with some special case data. – Farid Cheraghi May 12 '16 at 9:16

Stepwise regression is generally an unwise approach, even though commonly used and published. For supported inference, evaluating model support and model fit, the canned ArcGIS tools are not a good platform for specifying a statistical model. A one-size fits all approach to fitting regression models is never a good idea.

The message regarding the foibles of stepwise regression is currently being pushed by some of the current greats in statistics/biometrics (Harrell, Bolker, Bivand, Cressie, Hastie, Friedman, …). I believe that these limitations apply equally to stepwise AIC because it is a forward selection procedure. Keep in mind that the original intent behind AIC was hypothesis testing, using competing models, not feature selection.

Some of the identified issues with stepwise approaches are:

  • The r-square values are inflated in relation to the sample population
  • The test-statistics do not retain the correct distribution (F, chi-squared) which results in a deflation of p-values and standard errors.
  • Bonferroni corrections are overly conservative
  • Regression coefficients are upwards biased and marginal variables are more likely to be included if coefficient is overestimated and excluded if underestimated.
  • Removing "insignificant" variables sets their coefficients to zero, which may be implausible.

The general advice is to fit a full model, only remove parameters if α > 0.5, if the sign (+/-) is not sensible or if a coefficient of 0 is plausible. One should account for colinearity by combining variables (via mathematical function or PCA reduction).

If one really needs to reduce parameter space (eg., for estimation purposes) then a step down (backwards selection) approach is highly recommend. Backwards elimination methods can be implemented in R using the “fastbw” function in the "rms" library. As recommend by @whuber, I would highly recommend exploring a LASSO regression. If you are comfortable letting the data fit your model, you could implement a Bayesian Model Averaging approach. A Metropolis Hastings algorithm is a good choice when the parameter space is large. The R packages "BMA" and "BMS" have fairly simple implementations for Bayesian Model Selection.

Here is a simple method, implemented in R, for selecting parameter(s) with a p < 0.05 using a generalized linear model. This could be used to remove spurious independent variables. It would be easy to adapt this for an Ordinary Least Squares but, generally, in this case, I would use a maximum likelihood. I would think of this as more of a grab-bag approach than stepwise. However, I would evaluate the parameter(s) coefficents before applying an automatic selection procedure such as this.

Create some example data with one highly significant variable (x4)

x1 <- rnorm(10); x2 <- rnorm(10); x3 <- rnorm(10)
y <- rnorm(10)
x4 <- y + 5 
( dat <- as.data.frame(cbind(x1,x2,x3,x4,y)) )

Specify a Generalized Linear Model (GLM)

model <- glm(y ~ x1 + x2 + x3 + x4, data = dat) 
  round(summary(model)$coeff, 4)

Use p-value from t-test to select supported parameter(s)

( sel.x <- summary(model)$coeff[-1,4] < 0.05 ) 
( sig.x <- names(sel.x)[sel.x == TRUE] ) 

Specify new GLM with selected parameters

sig.model <- glm(reformulate(sig.x, response="y"), data = dat)
  round(summary(sig.model)$coeff, 4)

round(anova(model, sig.model),4)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.