Q:

Exploratory Spatial Regression: How to limit 60 variable candidates?

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?

A:

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

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.
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.

A:

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)

