1

I'm using a variogram to try and detect un-accounted for temporal auto-correlation in my model, but it shows the same pattern for all group levels (fSite). I allowed each fSite to have it's own trend, so I'm wondering if my code isn't creating a variogram for each fSite properly. CYR = year, fSeason = factor season ("WET", "DRY"), fSite = 1:47 locations, num = count data.

gam_szfs <- gam(num ~ 
                  s(CYR) +
                  s(CYR, fSeason, bs = "sz") +
                  
                  s(ToD) +
                  
                  s(DO) +
                  s(ave_tt) +
                  
                  s(CYR, fSite, bs="fs"),
                
                offset(log(area_sampled)), 
                
                data = toad2, 
                method = 'REML',
                select = FALSE,
                family = nb(link = "log"),
                control = list(trace = TRUE),
                drop.unused.levels=FALSE)

toad2$E1 <- residuals(gam_szfs, type = "deviance")  #' Residuals

#' I subset the database to make it easier to share
toad3 <- subset(toad2, fSite %in% c(1:5), 
                select = c(CYR, fSeason, fSite, num, E1))

MyData2 <- data.frame(E1   = toad3$E1, # Residuals
                      CYR = toad3$CYR,
                      fSite = toad3$fSite,
                      Ones = rep(1, nrow(toad3)))

coordinates(MyData2) <- c("CYR", "Ones")


# Initialize an empty list to hold variogram results
variogram_results <- list()

# Get unique fSites
unique_fSites <- unique(MyData2$fSite)

length(unique(toad3$CYR))

# Loop through each fSite, calculate the variogram, and store it
for (fSite in unique_fSites) {
  site_data <- subset(MyData2, fSite == fSite)
  
  V1 <- variogram(E1 ~ CYR, site_data, cutoff = 16, cressie = FALSE)

  # Add fSite information to the variogram result
  V1$fSite <- fSite
  
  # Append the result to the list
  variogram_results[[fSite]] <- V1
}


# Combine all the results into a single data frame
variogram_df <- do.call(rbind, variogram_results)

variogram_df$fSite <- as.numeric(variogram_df$fSite)


# Plot using ggplot2 with facet_wrap
ggplot(variogram_df, aes(x = dist, y = gamma)) + 
  geom_point() +
  geom_line() +
  facet_wrap(~ fSite, scales = "free") + 
  theme_minimal() +
  labs(title = "Variogram by fSite",
       x = "Distance",
       y = "Semivariance") +
  theme(strip.text = element_text(size = 12)) + 
  geom_smooth(method = "loess")

Variogram

Data (subset):

> dput(toad3)
structure(list(CYR = c(2008L, 2008L, 2008L, 2008L, 2008L, 2008L, 
2008L, 2008L, 2008L, 2008L, 2009L, 2009L, 2009L, 2009L, 2009L, 
2009L, 2009L, 2009L, 2009L, 2009L, 2010L, 2010L, 2010L, 2010L, 
2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2011L, 2011L, 2011L, 
2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2012L, 2012L, 
2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2013L, 
2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 
2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2015L, 2015L, 2015L, 
2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2015L, 2016L, 2016L, 
2016L, 2016L, 2016L, 2016L, 2016L, 2016L, 2016L, 2016L, 2017L, 
2017L, 2017L, 2017L, 2017L, 2017L, 2017L, 2017L, 2017L, 2017L, 
2018L, 2018L, 2018L, 2018L, 2018L, 2018L, 2018L, 2018L, 2018L, 
2018L, 2019L, 2019L, 2019L, 2019L, 2019L, 2019L, 2019L, 2019L, 
2019L, 2019L, 2020L, 2020L, 2020L, 2020L, 2020L, 2020L, 2020L, 
2020L, 2020L, 2020L, 2021L, 2021L, 2021L, 2021L, 2021L, 2022L, 
2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 
2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 
2023L, 2024L, 2024L, 2024L, 2024L, 2024L), fSeason = structure(c(1L, 
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 
2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 
1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 
2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 
2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 
2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 
1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L), levels = c("DRY", 
"WET"), class = "factor"), fSite = structure(c(1L, 2L, 3L, 4L, 
5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 
1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 
2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 
3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 1L, 2L, 
3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 
4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 
5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 
1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 
2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 
3L, 5L, 4L, 1L, 2L, 3L, 4L, 5L), levels = c("1", "2", "3", "4", 
"5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", 
"16", "17", "18", "19", "20", "21", "22", "23", "24", "25", "26", 
"27", "28", "29", "30", "31", "32", "33", "34", "35", "36", "37", 
"38", "39", "40", "41", "42", "43", "44", "45", "46", "47"), class = "factor"), 
    num = c(0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 
    0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 9L, 1L, 0L, 0L, 0L, 0L, 
    0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 2L, 0L, 0L, 
    0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
    0L, 2L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
    0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L), E1 = c(-0.247430445123603, 
    -0.22472807712087, -0.375974967991782, -0.249593492272773, 
    1.99017945741554, -0.801519670205052, -0.736040194113736, 
    -0.123161871252792, -0.895358458417623, 0.0494639074642558, 
    -0.249523011864716, -0.229733268283434, -0.305189512049232, 
    -0.345950930516627, -0.343089385209182, -0.631093065741142, 
    -0.583572354668092, -0.733598128382211, 0.84725333012146, 
    -0.763330641913427, -0.277668551241714, -0.193585056282638, 
    -0.432260633689169, -0.254316266229338, -0.248666488050697, 
    -1.06815739082575, -0.612826975973317, 0.100082933857372, 
    -0.80135441090358, -0.726054182378461, -0.425148377910311, 
    -0.148053095052621, -0.301389217660703, -0.292071269671274, 
    -0.239486807254894, 0.608813948659981, -0.701321582330218, 
    1.85085538983813, 0.604296487578992, -0.764765800267715, 
    -0.660765635587696, -0.264625481247797, -0.612658396105698, 
    -0.345263802121041, -0.318354950807924, -0.841343227715242, 
    -0.608530193895995, -1.19131742176282, 0.0710528702124216, 
    -0.595497263021366, -0.436596222887058, -0.3223713439725, 
    -0.71791434796134, -0.455175510100454, -0.367144943787426, 
    2.13816289256332, -0.347409136642035, -0.717334233994323, 
    -0.400052724630477, -0.372718053724584, -0.578478029873951, 
    -0.303460705957396, 0.740212211352622, -0.505096728402899, 
    -0.376361693421359, -0.41755717072046, -0.697744279470526, 
    -0.407916245163896, -0.59272948637825, -0.488056472560389, 
    -0.308905969831333, -0.597519263195212, -0.313561734339753, 
    -0.35994079627526, 2.50241206360777, -0.320189293258397, 
    0.830208225420878, -0.275218506548705, -0.524571213886048, 
    -0.434820227420392, 1.63138770054723, -0.648631480159925, 
    -0.198820453964549, -0.321511591826983, -0.697439631294457, 
    -0.324891498532058, -0.759245595025936, -0.405960581577857, 
    -0.508584847696458, -0.736064104089378, -0.531652121668709, 
    -0.498180862702523, -0.262224473535597, -0.339854643097809, 
    -0.396082824465216, -0.332226546716085, -0.657455020420682, 
    -0.253730346692792, -0.377288864370659, -0.462347770712867, 
    -0.356969880055938, -0.474332638603337, -0.220519611022768, 
    -0.457665352867075, -0.47508535461458, -0.341993532653451, 
    -0.441424307933895, -0.242274792268335, -0.355786532813558, 
    -0.402584979317751, -0.337554992649529, -0.45650232160878, 
    -0.24149419033347, -0.395035248029891, -0.438124422389352, 
    -0.379061581895076, -0.635182362577921, -0.24322209268536, 
    -0.261752528801338, -0.347877432982829, -0.263412127388931, 
    -0.532111539634991, -0.206747675294597, -0.251467216741288, 
    -0.398654457300848, -0.204833241746542, -0.400329864409977, 
    -0.120584180391436, -0.213059735948559, -0.293916885627306, 
    -0.21110175079172, -0.408456704811647, -0.169163650866059, 
    -0.23037288826383, -0.316055225876057, -0.211904110140088, 
    -0.260777881034904, -0.100062364875754, -0.166331729441553, 
    -0.191116746188383, -0.213469874700947, -0.619706455847356, 
    -0.159016176364313, -0.238631387346663, -0.345754981976852, 
    -0.228100665060894, -0.335259995479025, -0.0774414876508216, 
    -0.122961830664105, -0.135899888604339, -0.164185253978145, 
    1.42457872315356, -0.144794413408181, -0.22011605503903, 
    -0.283612468210839, -0.225082297174995)), row.names = c(10L, 
11L, 12L, 13L, 14L, 52L, 53L, 54L, 55L, 56L, 110L, 111L, 112L, 
113L, 114L, 126L, 127L, 128L, 129L, 130L, 190L, 191L, 192L, 193L, 
194L, 196L, 199L, 200L, 201L, 202L, 271L, 272L, 273L, 274L, 275L, 
323L, 328L, 329L, 330L, 331L, 369L, 370L, 371L, 372L, 373L, 424L, 
425L, 426L, 427L, 428L, 470L, 471L, 472L, 473L, 474L, 517L, 518L, 
519L, 520L, 521L, 554L, 555L, 556L, 557L, 558L, 592L, 642L, 643L, 
644L, 645L, 646L, 691L, 692L, 693L, 694L, 695L, 739L, 740L, 741L, 
742L, 743L, 780L, 781L, 782L, 783L, 784L, 798L, 801L, 802L, 806L, 
807L, 874L, 875L, 876L, 877L, 878L, 921L, 922L, 923L, 924L, 925L, 
938L, 939L, 940L, 941L, 942L, 982L, 983L, 984L, 985L, 986L, 1047L, 
1048L, 1049L, 1050L, 1051L, 1109L, 1110L, 1111L, 1112L, 1113L, 
1128L, 1129L, 1130L, 1131L, 1132L, 1204L, 1205L, 1206L, 1207L, 
1208L, 1220L, 1221L, 1222L, 1223L, 1224L, 1258L, 1261L, 1262L, 
1263L, 1264L, 1335L, 1336L, 1337L, 1338L, 1339L, 1374L, 1377L, 
1378L, 1379L, 1389L, 1403L, 1404L, 1405L, 1406L, 1407L), class = "data.frame")
Improve this question

1 Answer 1

Reset to default
3
+100

Non-standard evaluation is tripping you up. Your original code has subset(MyData2, fSite == fSite). R couldn't tell that the fSite on the left was supposed to refer to a column of MyData2 ... if you change the loop variable to fSite0, or use something like site_data[site_data$fSite == fSite, ], it should work ...

# Loop through each fSite, calculate the variogram, and store it
for (fSite0 in unique_fSites) {
  site_data <- subset(MyData2, fSite == fSite0)
  V1 <- gstat::variogram(E1 ~ CYR, site_data, cutoff = 16, cressie = FALSE)
  # Add fSite information to the variogram result
  V1$fSite <- fSite0
  # Append the result to the list
  variogram_results[[fSite0]] <- V1
}

facet plot of variograms by site

Improve this answer
Sign up to request clarification or add additional context in comments.

Comments

Your Answer

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

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.