2
$\begingroup$

I am currently working on a dataset (count data) in which one observation corresponds to one day of monitoring at a site. The overall protocol is to monitor groups of sites along transects. Almost all transects (and by extension almost all sites) were monitored several times. Sometimes, several transects were monitored at on given date.

I would like to be able to predict the abundance of the studied species according to multiple landscape and meteorological metrics. To do so, I tried to construct GLMMs (with a negative binomial distribution) using the glmmTMB function:

Abundance ~ environmental variables + meteorological variables + Julian day (and its quadratic effect) + year (in factor) + (1 | transect / site) + (1 | date)

However, my models suffer from spatial autocorrelation. I tested it with the testSpatialAutocorrelation function from the DHARMa package:

groupLocations = aggregate(countData[, which(colnames(countData) %in% c("x","y"))]
                           , list(countData$xy), mean)
groupLocations$xy <- unique(countData$xy)
res <- simulateResiduals(myModel)
res2 <- recalculateResiduals(res, countData$xy)
all(unique(res2$group) == groupLocations$xy)
testSpatialAutocorrelation(res2,x=groupLocations$x,y=groupLocations$y)

I would thus like to correct it, ideally using a covariance structure like those presented in this vignette : https://cran.r-project.org/web/packages/glmmTMB/vignettes/covstruct.html So I tried to add exp(pos + 0 | group) to my model with (1) countData$pos <- numFactor(countData$x, countData$y), (2) with group being a dummy variable, my sites IDs or my date IDs... Whatever the solution I choose, it doesn't seem to work (R is working indefinitely).

I guess my problem might come from two main issues:

  • the factI don't fully understand what the "group" metric sould be,
  • the fact that I have several observations at the same location (but at different dates),

For information, latitude (y) and longitude (x) are in Lambert 93 in my dataset but I can transfrom them in WGS84 if needed.

*EDIT Here is an example of the type of data I use (but not the real ones, which are very heavy) and my code. Interestingly with this small dataset, exp(pos + 0 |group) seems to work whether I define group as a dummy variable or as my site IDs, but with different results. In both cases, the testSpatialAutocorrelation function still detect spatial autocorrelation.

library(DHARMa)
library(glmmTMB)

rm(list=ls())

####I. Data creation####

newData <- data.frame(transect=c(rep("Transect_A",8*8),rep("Transect_B",10*3),rep("Transect_C",12)))
newData$site <- c(rep(paste0("Site_A",c(1:8)),8),rep(paste0("Site_B",c(1:10)),3),paste0("Site_C",c(1:12)))
newData$date <- c(rep("2001-05-03",8),rep("2001-06-25",8),rep("2002-06-04",8),rep("2002-07-15",8)
  ,rep("2003-04-28",8),rep("2003-05-16",8),rep("2003-06-30",8),rep("2003-07-24",8)
  ,rep("2001-06-02",10),rep("2002-07-15",10),rep("2003-05-08",10)
  ,rep("2033-05-08",12))

coordX_TA <- c(524623.6,524379.1,524379.1,524614.5,524877.1,525040.1,525185.0,525112.5)
coordX_TB <- c(527023.3,526769.7,526615.8,526652.0,526941.8,527258.7,527430.8,527566.6,527620.9,527783.9) 
coordX_TC <- c(525329.9,525148.8,525230.3,525447.6,525710.2,525945.7,526126.8,526172.0,526090.5,525927.5,525728.3,525683.0)

coordY_TA <- c(6705186,6705041,6704842,6704715,6704806,6704942,6705132,6705358)
coordY_TB <- c(6703728,6703593,6703357,6703094,6702941,6703022,6703194,6703420,6703665,6703864)
coordY_TC <- c(6700966,6700767,6700486,6700369,6700323,6700360,6700559,6700803,6701057,6701220,6701365,6701102)

newData$X <- c(rep(coordX_TA,8),rep(coordX_TB,3),coordX_TC)
newData$Y <- c(rep(coordY_TA,8),rep(coordY_TB,3),coordY_TC)

newData$abundance <- c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,2,0,0,0,0,0,2,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
                       ,0,0,7,0,0,0,0,0,0,2,0,0,0,0,0,4,2,4,0,1,12,0,0,0,0,0,4,0,0,0,0,0,0,0,2,3,0,0,0,0,0,0,1
                       ,0,0,0,1,0,0,0,0,0,0,1,0,9,3,6,11,4,9,12)

newData$explanatoryVar <- c(rep(c(0.42899904,0.77110546,0.45891875,2.35145233,0.52373614,0.03604363,-0.20136060,3.23943293),8)
                          ,rep(c(0.21227248,1.31673835,3.63072002,1.29581033,3.55550595,1.36627219
                                 ,1.56734995,-0.11314501,-0.38866418,-0.18451185),3)
                          ,c(-0.39731842,-0.47077649,-0.46944891,-0.18973198,-0.53603923,-0.38094487,0.16941656,0.60947665
                             ,0.08765610,0.33079765,0.09083049,2.18221584))

####II. Model without correction####

newData$explanatoryVar_scaled <- scale(newData$explanatoryVar)

myModel <- glmmTMB(abundance ~ explanatoryVar + (1|transect/site) + (1|date)
                   , data = newData, family = "nbinom2",na.action="na.fail")

summary(myModel)

#test Spatial autocorrelation 

groupLocations = aggregate(newData[, which(colnames(newData) %in% c("X","Y"))]
                           , list(newData$site), mean)
groupLocations$site <- unique(newData$site)
res <- simulateResiduals(myModel)
res2 <- recalculateResiduals(res, newData$site)
all(unique(res2$group) == groupLocations$site)
testSpatialAutocorrelation(res2,x=groupLocations$X,y=groupLocations$Y)

####III. Model with correction####

newData$pos <- numFactor(scale(newData$X), scale(newData$Y))
# newData$group <- factor(rep(1, nrow(newData))) #group as a dummy variable
newData$group <- factor(newData$site)

myModelCorrection <- glmmTMB(abundance ~ explanatoryVar + (1|transect/site) + (1|date) + exp(pos + 0 | group)
                   ,data = newData, family = "nbinom2",na.action="na.fail")

summary(myModelCorrection)

#test Spatial autocorrelation 

groupLocations = aggregate(newData[, which(colnames(newData) %in% c("X","Y"))]
                           , list(newData$site), mean)
groupLocations$site <- unique(newData$site)
res <- simulateResiduals(myModelCorrection)
res2 <- recalculateResiduals(res, newData$site)
all(unique(res2$group) == groupLocations$site)
testSpatialAutocorrelation(res2,x=groupLocations$X,y=groupLocations$Y)
Improve this question
$\endgroup$
12
  • $\begingroup$ Have you seen if you can reproduce the volcano example "A spatial covariance example" in the covstruct link, to understand the required structures? $\endgroup$
    – Alex J
    Commented Apr 20, 2023 at 1:32
  • $\begingroup$ Yes, I tried to reproduce this example by adding the same structure : exp(pos + 0 | group). Interestingly, for R not to run indefinitely, I have to define "group" as a dummy variable equal to 1 for instance (even if I have several observations by sites) and for my model to converge I have to define "pos" as my x and y scaled (i.e. longitude and latitude scaled beforehand). However, if I use the testSpatialAutocorrelation function on the residuals (using my x and y not scaled), it seems that this structure does not correct the autocorrelation issue, or even makes it worse... $\endgroup$
    – Lea_M
    Commented Apr 20, 2023 at 10:07
  • $\begingroup$ What if you group by both day and site? $\endgroup$
    – Claire Richards
    Commented Apr 20, 2023 at 17:49
  • 2
    $\begingroup$ Because transect sampling tends to detect species with different probabilities depending on their distance from the transect, special techniques are usually needed to analyze the results. This might be of more importance for the accuracy of your results than any other consideration. Some good explanations and advice are in Steven Thompson's book Sampling. $\endgroup$
    – whuber ♦
    Commented Apr 21, 2023 at 14:46
  • 1
    $\begingroup$ I don't have any particular thoughts, although I agree with @whuber and also recommend Buckland's Introduction to Distance Sampling. For this particular data set in the reprex, I'm not sure you've got enough data to fit a spatial correlation structure + nested intercepts to a nbinom. I realise you've got a much larger set, but try to simplify until you get something working, and then go from there. $\endgroup$
    – Alex J
    Commented Apr 23, 2023 at 12:32

1 Answer 1

Reset to default
2
$\begingroup$

@Lea_M, I took a look at it, and I'm still learning about how to use these functions myself. The help info says you can have residual spatial autocorrelation even if the model takes care of it. I added the command rotation="estimated" during the recalculation and the autocorrelation is no longer statistically significant.

res2 <- recalculateResiduals(res, newData$site, rotation="estimated") 
testSpatialAutocorrelation(res2,x=groupLocations$X,y=groupLocations$Y)

data:  res2
observed = 0.072138, expected = -0.034483, sd = 0.056906, p-value = 0.06098
alternative hypothesis: Distance-based autocorrelation

[Quantile command help][2]

That said, the residual plot doesn't look great pretty bad. and I'm not certain why.

Another (non-exclusive) method suggested is to condition the simulations on the random effects, but it states here that for glmmTMB, it's not yet possible to do that.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thank you very much for this answer! I should have taken a closer look at the testSpatialAutocorrelation function help! With the command rotation="estimated", it does indeed appear that the spatial autocorrelation is corrected! The only problem that remains, however, is that when I use my real dataset, I cannot define "group" as my site IDs (R runs indefinitely or more exactly several hours at least, which is problematic because I need to perform a dredge). Finally, indeed, the residual graphs are not great for my example, fortunately they are OK for my real dataset! $\endgroup$
    – Lea_M
    Commented Apr 24, 2023 at 9:59
  • 1
    $\begingroup$ I'm glad it helped. It does seem like you should be able to account for the fact that the sites are clustered within the transects for the spatial correlation. I was able to write exp(pos + 0 | transect/site) but for this data, the results were the same, though it gave me a longer correlation structure. I'm stabbing in the dark here, I guess, but it seems to me like the random slopes of the spatial correlation should be nested, too. $\endgroup$
    – Claire Richards
    Commented Apr 24, 2023 at 22:39
  • $\begingroup$ Your're right... I'll test this on my dataset ! $\endgroup$
    – Lea_M
    Commented Apr 27, 2023 at 8:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.