最新消息:雨落星辰是一个专注网站SEO优化、网站SEO诊断、搜索引擎研究、网络营销推广、网站策划运营及站长类的自媒体原创博客

loops - R: Variogram per group - Stack Overflow

programmeradmin4浏览0评论

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")

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")
发布评论

评论列表(0)

  1. 暂无评论