2
$\begingroup$

I am using logistic regression to model habitat suitability for bighorn sheep. I've followed the "purposeful selection of covariates" method outlined in Hosmer, D. W., S. Lemeshow, and R. X. Sturdivant. 2013. Applied logistic regression. Hoboken, NJ: John Wiley & Sons, Inc. This went well, and I came up with a model with predictors that make sense ecologically:

FinalMod<-glm(Presence~GV5_10*DistEscp + Eastness + SEI + TWI + Slope + bio9 + prec6 + tmin12, family = binomial(link = "logit"),data=sheep)

enter image description here

Here's the data. This includes a lot of variables that I excluded from the model fitting based on colinearity.

I then went ahead and did some model evaluation stuff as follows:

1) Leave-one-out cross-validation:

kcv<-cv.glm(sheep, FinalMod)

This gives me delta values of 0.1.

2) ROC/AUC using the same data I used to fit the model:

p <- predict(FinalMod, newdata = sheep, type = 'response')
pr<- prediction(p, sheep$Presence)
prf<- performance(pr, measure = 'tpr', x.measure = 'fpr')
plot(prf)

#AUC
auc<- performance(pr, measure = 'auc')
auc09_10<- auc@y.values[[1]]

Which gave me an AUC of 0.94.

3) ROC/AUC of data from previous observation years that had not been used to fit the model, with the following AUC values: 1985 = 0.97, 1986 = 0.96.

4) An assessment of model overfitting from Frank Harrell, with a g^ of 0.94 indicating ~6% overfitting.

At this point, I'm very happy with my model, so I went to ArcGIS to use the raster calculator to make a probability surface. However, when I put in the equation to get probability, with the general form $\hat{p} = {exp(\hat{\beta_0} + x_1\hat{\beta_1} +...+ x_k\hat{\beta_k})\over{1+exp(\hat{\beta_0} + x_1\hat{\beta_1} +...+ x_k\hat{\beta_k}})}$, I came up with probabilities of 0.97-1.0 for the entire output raster. I also get this result when I apply the equation to the data in excel- $p\approx1$ for all records, including absences. I have entered and re-entered the formula, I'm pretty confident it isn't a typo that's giving me grief. Here is the formula that I've been using:

(exp(0.401 - 0.00823*GV5-10 - 0.0041*DistEscp + 0.7356*Eastness + 0.03364*SEI + 0.1812*TWI + 0.01635*Slope + 0.0497*bio9 + 0.4012*prec6 + 0.00002611*GV5-10*DistEscp)) 
/(1+(exp(0.401 - 0.00823*GV5-10 - 0.0041*DistEscp + 0.7356*Eastness + 0.03364*SEI + 0.1812*TWI + 0.01635*Slope + 0.0497*bio9 + 0.4012*prec6 + 0.00002611*GV5-10*DistEscp))

So, I am wondering why I'm getting ridiculous results when I try and apply the equation by hand, but good results from the assessments in R.

Final notes: There is one interaction term, which I was using in the formula as multiplication: ($\hat{\beta}_{12}(x_1x_2)$). I have a fairly small sample size- roughly 250 occurrence points, with a few more than that randomly selected 'background' points for a total of about 500 observations. Also, all explanatory variables are continuous.

Thanks for any help, and please let me know if I need to include more information. I'm trying to keep this as concise as possible.

$\endgroup$
  • $\begingroup$ Without the data it is hard to provide an answer specific to your situation. In general though, the AUC can be very high ('good'), while the actual predicted probabilities are way off compared to the observed ones. See this wiki (en.wikipedia.org/wiki/…). If you are concerned about the predicted probabilities specifically, look for calibration measures instead. Note: this is assuming your model and calculations are correct. $\endgroup$ – IWS May 29 '17 at 7:49
  • $\begingroup$ Thanks for the quick response, IWS. I've added a link to my data and the full equation I've been using. Do the low delta and adjusted delta from the cross validation fall into the same category as the AUC? I had been under the impression that those numbers were telling me that there was about a 10% missclassification rate with a cutoff of p=0.5. Is that incorrect? $\endgroup$ – Brody May 29 '17 at 16:41
  • $\begingroup$ You reference two different models above. I'm wondering if this has something to do with it? In your cross-validation, you reference "FinalMod1," but in you called the model "FinalMod" initially (no 1). When I ran the model for FinalMod, I didn't get coefficients anywhere near yours. Is it possible you have mixing up your models FinalMod and FinalMod1 perhaps? $\endgroup$ – StatsStudent May 29 '17 at 16:46
  • 1
    $\begingroup$ The typo @whuber found does fix the problem, thanks very much for catching that. I'll put that together as an answer. Thanks again for your help. $\endgroup$ – Brody May 29 '17 at 18:48
  • 1
    $\begingroup$ Thanks for letting us know about your success, Brody. I would assign a high probability to there still being errors in the ArcGIS calculation. That is because it should have choked in parsing the "GV5" expression. Because it didn't, I suspect there could be issues about differences in how that is coded in R and ArcGIS. You need to check your work by--at a minimum--computing predictions for all the original data both in R and in ArcGIS and checking that they agree to within ArcGIS's limited (single) precision. $\endgroup$ – whuber May 29 '17 at 18:52
0
$\begingroup$

I suspect the problems you are experiencing with your hand calculations in Excel (and ArcGIS) are result of Excel's well-known numerical accuracy issues, which make it ill-suited to handle regression problems. See here for more details on this: http://www.pucrs.br/famat/viali/tic_literatura/artigos/planilhas/pottel.pdf. When I perform "manual" calculations in R, I have no problem with this getting the same results as those from FinalMod in R.

Try this in R, for example:

    newsheep<-data.frame("ItcTerm"=1, sheep$GV5_10, sheep$DistEscp, sheep$Eastness, sheep$SEI, sheep$TWI, sheep$Slope, sheep$bio9, sheep$prec6, sheep$tmin12, sheep$GV5_10*sheep$DistEscp)
    betaHatTimesPredictor<-as.matrix(newsheep)%*%as.matrix(FinalMod$coef)
    ManualPhat<-exp(betaHatTimesPredictor)/(1+exp(betaHatTimesPredictor))

    #Now Compare Manual Calculations vs. GLM Calculations
    head(cbind(ManualPhat, FinalMod$fit))
             [,1]        [,2]
    1 0.002263072 0.002263072
    2 0.005236063 0.005236063
    3 0.003695555 0.003695555
    4 0.041600082 0.041600082
    5 0.001737018 0.001737018
    6 0.000795582 0.000795582

Now compare those calculations to ones you obtain from Excel. You might consider standardizing your calculations if you'd like to perform manual calculations in Excel to obtain your predicted values. Also, please be sure you are careful about what terms you have included in your model. As @whuber rightly pointed out, your post has terms appearing in parts of your model and then disappearing later in your post. Be sure to clean those up first and verify that the model you have built in R is indeed the one you have performed hand calculations on in excel.

$\endgroup$
  • $\begingroup$ Excel has no accuracy issues when it comes to basic arithmetic: it uses the native IEEE double-precision calculations. Indeed, it can accurately reproduce R calculations, which (for basic arithmetic) work exactly the same. There will be an issue with the Raster Calculator expressions, because ArcGIS likes to use single-precision floats: but those will still agree with R to six significant figures. $\endgroup$ – whuber May 29 '17 at 18:54
  • $\begingroup$ @whuber, correct. I understood the OP to be running regression in Excel or via one of its add-ons, -- not simply multiplying the terms he obtained in R and adding them up, which of course, should be fine. Then he was calculating $\hat{p}$ manually (i.e. paper and pen or with a calculator) and seeing discrepancies. Perhaps I was mistaken? Maybe OP can clarify? Either way, it looks with your help he solved the issue! $\endgroup$ – StatsStudent May 29 '17 at 18:57
  • $\begingroup$ The situation is a little different. ArcGIS is a platform that maintains raster data of all the regressors. R was used to fit a model for predicting the response based on a sample of the regressors and some observed responses for them. Its results needed to be applied by ArcGIS to the raster data in order to create a raster of the predicted response variables ("make a probability surface"). There were no new regressions being performed. $\endgroup$ – whuber May 29 '17 at 19:00
  • $\begingroup$ @whuber, I don't use ArcGIS, so I'm hoping to learn something here: can't one simply feed one's predicted probabilities into ArcGIS? Maybe that's what you are saying? $\endgroup$ – StatsStudent May 29 '17 at 19:02
  • 1
    $\begingroup$ Yes, there are various ways to do that: ArcGIS can be linked directly to R and it can export and import raster data. However, under some circumstances R will fail to perform a raster calculation, such as when there are huge datasets involved: the add-ons for R tend to be much slower and more cumbersome. (This circumstance is remarkably similar to a story I related at the end of an answer at stats.stackexchange.com/a/7933/919.) $\endgroup$ – whuber May 29 '17 at 19:05

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.