[R-sig-ME] Spatial correlation in lme

[Ben Bolker]

 <Jeffrey_Warren at ...> writes:

> 
> 
> We've created 35 home ranges for individual birds that we tracked during
> the pre-breeding season. Each home range is based on an individual
> utilization distribution (UD) that we've used to predict relative at each
> 100 x 100 m pixel. Each pixel also has associated habitat attributes, such
> as water depth and percent emergent vegetation. We also have bird-level
> variables (age [categorical] and relative body condition [continuous]).
> We're now trying to relate habitat use (UD value) to habitat attributes and
> individual quality (with age and body condition as proxies). We've modeled
> the data with individual bird as a random effect to account for the related
> nature of pixels within an individual's home range. We then imposed a
> rational quadratic correlation structure to the lme model to account for
> the spatial autocorrelation among pixels (corRatio was best fit compared to
> other structures available in nlme). When we look at variograms for each
> model we don't see any improvement in residual correlation, but our beta
> estimates for the effects of habitat on selection change significantly.

When you say "look at the variograms", do you mean that you're doing
plot(Variogram(fittedmodel)) ?  If so, you might want to take a careful
look at ?Variogram.lme, which says among other things:

> resType: an optional character string specifying the type of residuals
          to be used. If ‘"response"’, the "raw" residuals (observed -
          fitted) are used; else, if ‘"pearson"’, the standardized
          residuals (raw residuals divided by the corresponding
          standard errors) are used; else, if ‘"normalized"’, the
          normalized residuals (standardized residuals pre-multiplied
          by the inverse square-root factor of the estimated error
          correlation matrix) are used. Partial matching of arguments
          is used, so only the first character needs to be provided.
          Defaults to ‘"pearson"’.

  That is, if you type plot(Variogram(fittedmodel)) it does *not*
plot the variogram of the residuals corrected for autocorrelation:
it plots the variogram of the residuals corrected *only* for
differential variance ("pearson").  I believe you would need
plot(Variogram(fittedmodel,resType="normalized")) to look
at the variogram of the autocorrelation-corrected values ...

> Here is what our data look like:
> 
> >head(hrdata)
> 
>     XMIN   XMAX    YMIN    YMAX  BirdID DEP SUB H2O EDGE    UDval Age
> BCIndex Year_  logUDval
> 1 430400 430500 4942900 4943000 1653957   0   0   0    0 0.012085   1
> -18.82793     0 -4.415790
> 2 430400 430500 4943000 4943100 1653957   0   0   0    0 0.012856   1
> -18.82793     0 -4.353945
> 3 430400 430500 4943100 4943200 1653957   0   0   0    0 0.014132   1
> -18.82793     0 -4.259314
> 4 430400 430500 4943200 4943300 1653957   0   0   0    0 0.014939   1
> -18.82793     0 -4.203780
> 5 430400 430500 4943300 4943400 1653957   0   0   0    0 0.016364   1
> -18.82793     0 -4.112671
> 6 430400 430500 4943400 4943500 1653957   0   0   0    0 0.017369   1
> -18.82793     0 -4.053068

  You might want to center/scale your X/Y values; I have had some trouble
with that in the past ...

> Model statements:
> Form <- formula(logUDval ~ DEP + EDGE + SUB:H2O + Year_ + Age + BCIndex)
> m1.lme<-lme(Form, random = ~ 1|BirdID, method="REML", data=hrdata)
> m1.lme.ratio<-lme(Form, random = ~ 1|BirdID, correlation = corRatio(form =
> ~XMIN + YMIN), method="REML", data=hrdata)
> 
> Model summaries:
> 
> Summary of LME with no spatial autocorrelation structure
> 
> Linear mixed-effects model fit by REML
>  Data: hrdata
>        AIC      BIC    logLik
>   50535.09 50605.97 -25258.54
> 
> Random effects:
>  Formula: ~1 | BirdID
>         (Intercept)  Residual
> StdDev:   0.6385548 0.8809287
> 
> Fixed effects: list(Form)
>                   Value             Std.Error         DF    t-value
> p-value
> (Intercept)       -1.8712232  0.20935222  19415  -8.938158        0.0000
> DEP               0.0014205   0.00028237  19415       5.030507    0.0000
> EDGE              0.0028367   0.00029345  19415       9.666943    0.0000
> Year_1      -0.3690599  0.23942499        31    -1.541443   0.1334
> Age1        -0.2452163  0.23423909        31    -1.046863   0.3033
> BCIndex           0.0018829   0.00283624        31    0.663880    0.5117
> SUB:H2O           0.0037672   0.00030599  19415       12.311448   0.0000

  [snip]
 
> Summary of LME with rational quadratic spatial autocorrelation structure
> 
> Linear mixed-effects model fit by REML
>  Data: hrdata
>         AIC       BIC   logLik
>   -74203.99 -74125.23 37111.99

  120000 AIC units is a pretty big difference, even with 
20000 data points.  Hmmm ...
 
> Random effects:
>  Formula: ~1 | BirdID
>         (Intercept)  Residual
> StdDev:   0.5266185 0.5505424
> 
> Correlation Structure: Rational quadratic spatial correlation
>  Formula: ~XMIN + YMIN | BirdID
>  Parameter estimate(s):
>    range
> 405.6316

  Is this range value reasonable?

> Fixed effects: list(Form)
>                   Value             Std.Error         DF
> t-value           p-value
> (Intercept)  -2.4777143 0.18583575  19415             -13.332819   0.0000
> DEP          -0.0000045       0.00001259  19415             -0.357142
> 0.7210
> EDGE          0.0000032       0.00000281  19415             1.150247
> 0.2501
> Year_1       -0.3041446       0.21134054        31          -1.439121
> 0.1601
> Age1              -0.2356964  0.20706315        31          -1.138283
> 0.2637
> BCIndex            0.0013983  0.00250764        31           0.557623
> 0.5811
> SUB:H2O     -0.0000001  0.00000827  19415             -0.006635   0.9947

  [snip]
 
> We don't understand why there appears to be no improvement in spatial
> autocorrelation but our results have changed so dramatically. Any ideas on
> why we're running into this issue would be greatly appreciated.