I have tried to find answers to my question by googling but haven't found anything relevant - if anyone could pass on some resources or relevant key words that would be brilliant, in absence of actually getting an answer to my question.

My model predicts deforestation at pixel resolution using both pixel- and lower-resolution predictors. The formula is as follows:

gam(formula = loss ~ ecozone + s(dist_river) + s(dist_prevdf) + s(allarea) + s(dist_road) + s(altitude) + s(slope) + s(aspect, bs = "cc") + s(pcARA, k = 3) + s(popdens) + s(pcpriorloss) + s(pcforestGFC), family = binomial, data = gfcpix, weights = weight/mean(weight), method = "REML", select = TRUE, knots = list(aspect = c(0, 360)))

Here, aspect is an angle and therefore given a cyclic spline, and the weights argument refers to loss having been oversampled by a factor of 5 in the data as it was very rare to begin with. However, as the maximum pixels sampled at the unit of the lower-resolution predictors was constrained, no-loss pixels were variably undersampled. Weights therefore differ between loss and no loss pixels across communities. Instead of including communities as random effects, as our main question relates to a few predictors at this resolution, I therefore added these predictors as-is. The total number of data rows in the model is 123288.

I used Faraway's logistic model instructions for evaluating the model and grouped predicted probabilities into intervals, comparing these groups into the corresponding proportion of 1 in the data. Predicted vs observed probabilities I multiplied probabilities by 5 to correspond to the sample rather than population level of 1s.

I am unsure why the model overpredicts at high probabilities. Plotting binned residuals against individual predictors did not suggest any clear issues, but I assume this could have been one way of chasing the issue? I attempted adding some plausible interaction structures into the model but they were not significant nor did they impact on the predictions. Of course, I may be missing some important predictors that would help in constraining the predictions. However, I am at a loss in terms of how to troubleshoot this issue. What might be the most likely cause of the problem, and how could I test this plus other possible causes? Simply saying "The model overpredicted loss at high probabilities" seems inadequate.

Any help would be much appreciated.

I have not provided data or code in this case as my question is more generic than my example, however can try to generate a similarly 'flawed' model with invented data if necessary.

UPDATE: To clarify the situation and partially respond to suggestions in previous comments:

  • I've plotted binned partial residuals as well for reference. I'm no expert on interpreting logistic model binned residuals but again following Faraway's book, I would say there is nothing SERIOUS in here (looks mostly like outlier departures rather than systematic issues?). Residuals vs predictors

  • I redrew the predicted vs observed bins for various versions of the model, and it did seem that something was awry with the weights argument. Varying this argument alone by switching it off changed the model from having the bias to not having the bias. Playing with using the probit vs logit link also had a slight effect (when keeping weights as is).Predicted vs observed multiple model versions

  • In a way it is a solution that the weights argument is causing problems, but as I don't understand WHY I've not shut this case as resolved. High variability and low mean weights (0.2 value for all loss pixels, variable values >=1 for no loss pixels) characterises the area of bias, so can I assume that the low weights and fewer observations at high probability together combine to make observations in this range count less in the model fitting, allowing for a wobble?Weights relationship with predictions

  • 1
    $\begingroup$ You might want to look at the summary output for your model and see if any of the s(.) terms have a model degrees of freedom close to the upper bound k, which defaults to 10 I believe. That could indicate that there's more wiggliness in the relationship between that individual predictor and the target than k is allowing to fit, and often this will manifest itself in a poorly fit tail. (Just an added suggestion to the ones your question already contains.) $\endgroup$ – jbowman Jan 16 '18 at 16:39
  • 1
    $\begingroup$ To echo @jbowman's comment - look at the output from gam.check() to confirm whether the defaults basis size is sufficiently large (you want the p values in the table to be larger than 0.05) and for any smoother indicated to have too-low a basis dimension, double k and refit the model. You can also plot residuals vs each predictor and or fit a Gaussian GAM to the residuals with one covariate at a time to diagnose if you have over smoothed. $\endgroup$ – Gavin Simpson Jan 16 '18 at 22:03
  • 1
    $\begingroup$ @jbowman, @gavin, thanks for your suggestions. I've added information in this regard to my question. Note that I had previously looked at summary() output for edf and ref.df without feeling there were major issues. But this brings up a question of why gam.check is throwing errors, preventing me to see anything there. The function works for the option without weights in. Error in dm[, i] <- sort(residuals(object, type = type)) : replacement has length zero In addition: Warning message: In qbinom(p, wt, mu) : NaNs produced $\endgroup$ – Emma Wiik Feb 7 '18 at 12:45

Your Answer

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

Browse other questions tagged or ask your own question.