title: How to deal with spatially autocorrelated residuals in GLMM
body: <p>I am conducting an analysis of where on the landscape a predator encounters potential prey. My response data is binary with an Encounter location = 1 and a Random location = 0 and my independent variables are continuous but have been rescaled. </p>
  
  <p>I originally used a GLM structure </p>
  
  <pre><code>glm_global &lt;- glm(Encounter ~ Dist_water_cs+coverMN_cs+I(coverMN_cs^2)+
                    Prey_bio_stand_cs+Prey_freq_stand_cs+Dist_centre_cs,
                    data=Data_scaled, family=binomial)
  </code></pre>
  
  <p>but realized that this failed to account for potential spatial-autocorrelation in the data (a spline correlogram showed high residual correlation up to ~1000m). </p>
  
  <pre><code>Correlog_glm_global &lt;- spline.correlog (x = Data_scaled[, "Y"],
                                          y = Data_scaled[, "X"],
                                          z = residuals(glm_global,
                                          type = "pearson"), xmax = 1000)
  </code></pre>
  
  <p>I attempted to account for this by implementing a GLMM (in lme4) with the predator group as the random effect. </p>
  
  <pre><code>glmm_global &lt;- glmer(Encounter ~ Dist_water_cs+coverMN_cs+I(coverMN_cs^2)+
                       Prey_bio_stand_cs+Prey_freq_stand_cs+Dist_centre_cs+(1|Group),
                       data=Data_scaled, family=binomial)
  </code></pre>
  
  <p>When comparing AIC of the global GLMM (1144.7) to the global GLM (1149.2) I get a Delta AIC value >2 which suggests that the GLMM fits the data better. However I am still getting essentially the same correlation in the residuals, as shown on the spline correlogram for the GLMM model).</p>
  
  <pre><code>Correlog_glmm_global &lt;- spline.correlog (x = Data_scaled[, "Y"],
                                           y = Data_scaled[, "X"],
                                           z = residuals(glmm_global,
                                           type = "pearson"), xmax = 10000)
  </code></pre>
  
  <p>I also tried explicitly including the Lat*Long of all the locations as an independent variable but results are the same. 
  After reading up on options, I tried running Generalized Estimating Equations (GEEs) in “geepack” thinking this would allow me more flexibility with regards to explicitly defining the correlation structure (as in GLS models for normally distributed response data) instead of being limited to compound symmetry (which is what we get with GLMM). However I realized that my data still demanded the use of compound symmetry (or “exchangeable” in geepack) since I didn’t have temporal sequence in the data. When I ran the global model</p>
  
  <pre><code>gee_global &lt;- geeglm(Encounter ~ Dist_water_cs+coverMN_cs+I(coverMN_cs^2)+
                       Prey_bio_stand_cs+Prey_freq_stand_cs+Dist_centre_cs,
                       id=Pride, corstr="exchangeable", data=Data_scaled, family=binomial)
  </code></pre>
  
  <p>(using scaled or unscaled data made no difference so this is with scaled data for consistency)
  suddenly none of my covariates were significant. However, being a novice with GEE modelling I don’t know a) if this is a valid approach for this data or b) whether this has even accounted for the residual autocorrelation that has been evident throughout. </p>
  
  <p>I would be most appreciative for some constructive feedback as to 1) which direction to go once I realized that the GLMM model (with predator group as a random effect) still showed spatially autocorrelated Pearson residuals (up to ~1000m), 2) if indeed GEE models make sense at this point and 3) if I have missed something in my GEE modelling. Many thanks.  </p>