0

I am trying to plot a semivariogram of my model residuals for a generalised mixed effect model in R. Doing this for a mixed effect model with normal distribution is straightforward with the nlme package, and using the quakes dataset as an example.

library(nlme)
data(quakes)
head(quakes)

model1 <- lme(mag ~ depth  , random = ~1|stations, data = quakes)
summary(model1)
semivario <- Variogram(model1, form = ~long+lat,resType = "normalized")
plot(semivario, smooth = TRUE)

I want to create a model with a non-normal distribution, which I can't do with nlme, so I have tried glmer and glmmPQL. I have turned the 'mag' into a binomial variable, then try to reapply the Variogram function to make a plot with models.

quakes$thresh <- ifelse(quakes$mag > "5", 0, 1)

library(MASS)
model2 <- glmmPQL(as.factor(thresh) ~ depth  , random = ~1|stations, family = binomial, data = quakes)
summary(model2)
semivario <- Variogram(model2, form = ~long+lat,resType = "normalized")
plot(semivario, smooth = TRUE)

library(lme4)
model3 <-glmer(as.factor(thresh) ~ depth + (1|stations), data = quakes, family = binomial)
summary(model3)
semivario <- Variogram(model3, form = ~long+lat,resType = "normalized")
plot(semivario, smooth = TRUE)

Neither of these appear to work for plotting the variogram. The glmmPQL says that lat and long isn't found, and the glmer says distance isn't specified.

How can I code a plot of semivariogram of these models? Is the Variogram function from the nlme package unusable for them? And if so what alternatives can I use?

willbutdontwanna
  • 3
  • 2
  • best way to do this is probably to get the residuals from the NB fit, fit a trivial (intercept + grouping only) `lme` model, then apply `Variogram()` to the result ... – Ben Bolker Sep 02 '21 at 16:18
  • Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking. – Community Sep 02 '21 at 16:37
  • if you can provide a [mcve] that would be great (unfortunately I'm not immediately aware of publicly available/handy data sets that would support a spatial approach to a negative binomial mixed model, although something could obviously be simulated ... – Ben Bolker Sep 02 '21 at 19:16
  • There are two ways to do this, I think. (1) use `MASS::glmmPQL`. (2) Use my example above. If you can provide a [mcve] I can take a shot at showing an answer ... – Ben Bolker Sep 03 '21 at 22:22
  • thanks again. for some reason the glmPQL didn't work either. I have updated the question with an MRE. let me know if this is more helpful/clear – willbutdontwanna Sep 05 '21 at 13:29

0 Answers0