## ----include = FALSE---------------------------------------------------------- knitr::opts_chunk$set( comment = "#>" ) ## ----echo = FALSE, message = FALSE-------------------------------------------- library(occumb) set.seed(1) ## ----eval = FALSE------------------------------------------------------------- # data(fish) # summary(fish) ## ----------------------------------------------------------------------------- #> Sequence read counts: #> Number of species, I = 50 #> Number of sites, J = 50 #> Maximum number of replicates per site, K = 3 #> Number of missing observations = 6 #> Number of replicates per site: 2.88 (average), 0.33 (sd) #> Sequencing depth: 77910 (average), 98034.7 (sd) #> #> Species covariates: #> mismatch (continuous) #> Site covariates: #> riverbank (categorical) #> Replicate covariates: #> (None) #> #> Labels for species: #> Abbottina rivularis, Acanthogobius lactipes, Acheilognathus macropterus, Acheilognathus rhombeus, Anguilla japonica, Biwia zezera, Carassius cuvieri, Carassius spp., Channa argus, Ctenopharyngodon idella, Cyprinus carpio, Gambusia affinis, Gnathopogon spp., Gymnogobius castaneus, Gymnogobius petschiliensis, Gymnogobius urotaenia, Hemibarbus spp., Hypomesus nipponensis, Hypophthalmichthys spp., Hyporhamphus intermedius, Ictalurus punctatus, Ischikauia steenackeri, Lepomis macrochirus macrochirus, Leucopsarion petersii, Megalobrama amblycephala, Micropterus dolomieu dolomieu, Micropterus salmoides, Misgurnus spp., Monopterus albus, Mugil cephalus cephalus, Mylopharyngodon piceus, Nipponocypris sieboldii, Nipponocypris temminckii, Opsariichthys platypus, Opsariichthys uncirostris uncirostris, Oryzias latipes, Plecoglossus altivelis altivelis, Pseudogobio spp., Pseudorasbora parva, Rhinogobius spp., Rhodeus ocellatus ocellatus, Salangichthys microdon, Sarcocheilichthys variegatus microoculus, Silurus asotus, Squalidus chankaensis biwae, Tachysurus tokiensis, Tanakia lanceolata, Tribolodon brandtii maruta, Tribolodon hakonensis, Tridentiger spp. #> Labels for sites: #> 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, 48, 49, 50 #> Labels for replicates: #> L, C, R ## ----eval = FALSE------------------------------------------------------------- # fish # Result not shown ## ----eval = FALSE------------------------------------------------------------- # fit0 <- occumb(data = fish, parallel = TRUE) ## ----eval = FALSE------------------------------------------------------------- # fit1 <- occumb( # formula_psi = ~riverbank, # formula_phi_shared = ~mismatch, # data = fish, # parallel = TRUE # ) ## ----eval = FALSE------------------------------------------------------------- # fit1x <- occumb( # formula_psi = ~riverbank, # formula_phi_shared = ~mismatch, # data = fish, # n.thin = 20, # n.iter = 40000, # parallel = TRUE # ) ## ----eval = FALSE------------------------------------------------------------- # fit1 # Result not shown ## ----eval = FALSE------------------------------------------------------------- # plot(fit1) # Outputs many figures; enter "Esc" or "Ctrl + C" to quit ## ----eval = FALSE------------------------------------------------------------- # summary(fit1) ## ----------------------------------------------------------------------------- #> Summary for an occumbFit object #> #> Summary of data: #> Number of species, I = 50 #> Number of sites, J = 50 #> Maximum number of replicates per site, K = 3 #> Number of missing observations = 6 #> Number of replicates per site: 2.88 (average), 0.33 (sd) #> Sequencing depth: 77910 (average), 98034.7 (sd) #> #> Model specification: #> formula_phi: ~ 1 #> formula_theta: ~ 1 #> formula_psi: ~ riverbank #> formula_phi_shared: ~ mismatch #> formula_theta_shared: ~ 1 #> formula_psi_shared: ~ 1 #> prior_prec: 1e-04 #> prior_ulim: 10000 #> #> Saved parameters: #> Mu sigma rho alpha beta gamma alpha_shared phi theta psi z pi deviance #> #> MCMC ran in parallel for 43.593 minutes at time 2023-09-11 20:16:31.787026: #> For each of 4 chains: #> Adaptation: 100 iterations (sufficient) #> Burn-in: 10000 iterations #> Thin rate: 10 iterations #> Total chain length: 20100 iterations #> Posterior sample size: 1000 draws #> #> Summary of posterior samples: #> Mu: #> Number of parameters: 4 #> Rhat: 1.009 (min), 1.027 (median), 1.043 (mean), 1.109 (max) #> n.eff: 30 (min), 117.5 (median), 154.2 (mean), 352 (max) #> sigma: #> Number of parameters: 4 #> Rhat: 1.004 (min), 1.04 (median), 1.074 (mean), 1.213 (max) #> n.eff: 17 (min), 67.5 (median), 217 (mean), 716 (max) #> rho: #> Number of parameters: 6 #> Rhat: 1.042 (min), 1.107 (median), 1.106 (mean), 1.187 (max) #> n.eff: 19 (min), 33.5 (median), 40.7 (mean), 79 (max) #> alpha: #> Number of parameters: 50 #> Rhat: 1.027 (min), 1.093 (median), 1.157 (mean), 1.77 (max) #> n.eff: 7 (min), 34 (median), 40.6 (mean), 104 (max) #> beta: #> Number of parameters: 50 #> Rhat: 1.002 (min), 1.018 (median), 1.028 (mean), 1.237 (max) #> n.eff: 16 (min), 183.5 (median), 349.9 (mean), 4000 (max) #> gamma: #> Number of parameters: 100 #> Rhat: 1.001 (min), 1.019 (median), 1.029 (mean), 1.136 (max) #> n.eff: 24 (min), 173.5 (median), 363.5 (mean), 2756 (max) #> alpha_shared: #> Number of parameters: 1 #> Rhat: 1.1 #> n.eff: 31 #> phi: #> Number of parameters: 50 #> Rhat: 1.001 (min), 1.07 (median), 1.127 (mean), 1.664 (max) #> n.eff: 8 (min), 42.5 (median), 202.7 (mean), 3161 (max) #> theta: #> Number of parameters: 50 #> Rhat: 1.002 (min), 1.016 (median), 1.024 (mean), 1.219 (max) #> n.eff: 16 (min), 196 (median), 358.1 (mean), 4000 (max) #> psi: #> Number of parameters: 2500 #> Rhat: 1 (min), 1.008 (median), 1.013 (mean), 1.101 (max) #> n.eff: 36 (min), 396 (median), 679.3 (mean), 4000 (max) #> z: #> Number of parameters: 2500 #> Rhat: 1 (min), 1.006 (median), 1.009 (mean), 1.184 (max), 946 (Number of NAs) #> n.eff: 1 (min), 1233.5 (median), 1628.6 (mean), 4000 (max) #> pi: #> Number of parameters: 7500 #> Rhat: 1 (min), 1.078 (median), 1.102 (mean), 1.294 (max), 30 (Number of NAs) #> n.eff: 1 (min), 1 (median), 1009.1 (mean), 4000 (max) #> deviance: #> Rhat: 1.001 #> n.eff: 1332 ## ----eval = FALSE------------------------------------------------------------- # post_summary_psi <- get_post_summary(fit1, "psi") # post_summary_psi ## ----------------------------------------------------------------------------- #> mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff overlap0 f #> psi[1,1] 0.64607375 0.083034133 0.466158634 0.59382744 0.65078275 0.70478301 0.79318634 1.001697 1264 0 1 #> psi[2,1] 0.28250539 0.089275014 0.136970706 0.21789079 0.27480943 0.33672122 0.47971573 1.013796 317 0 1 #> psi[3,1] 0.38321762 0.088583160 0.218908875 0.32150169 0.38055944 0.44171514 0.56510758 1.005127 472 0 1 #> psi[4,1] 0.16899696 0.070833671 0.058363836 0.11684541 0.15943020 0.21075828 0.32821669 1.001143 3788 0 1 #> psi[5,1] 0.20435372 0.088327893 0.069729103 0.13950376 0.19315615 0.25286705 0.41681236 1.006308 596 0 1 #> psi[6,1] 0.10737782 0.061179251 0.024299550 0.06246265 0.09573034 0.13925811 0.25577838 1.002291 1418 0 1 #> psi[7,1] 0.84639119 0.053478041 0.730770215 0.81327035 0.85048238 0.88494478 0.93751979 1.003575 654 0 1 #> psi[8,1] 0.94296846 0.028288919 0.875373936 0.92739641 0.94770976 0.96354354 0.98337253 1.002581 1006 0 1 #> (Omitted the remaining) ## ----eval = FALSE------------------------------------------------------------- # attributes(post_summary_psi)$dimension ## ----------------------------------------------------------------------------- #> [1] "Species" "Site" ## ----eval = FALSE------------------------------------------------------------- # attributes(post_summary_psi)$label ## ----------------------------------------------------------------------------- #> $Species #> [1] "Abbottina rivularis" #> [2] "Acanthogobius lactipes" #> [3] "Acheilognathus macropterus" #> [4] "Acheilognathus rhombeus" #> [5] "Anguilla japonica" #> [6] "Biwia zezera" #> [7] "Carassius cuvieri" #> [8] "Carassius spp." #> [9] "Channa argus" #> [10] "Ctenopharyngodon idella" #> [11] "Cyprinus carpio" #> [12] "Gambusia affinis" #> [13] "Gnathopogon spp." #> [14] "Gymnogobius castaneus" #> [15] "Gymnogobius petschiliensis" #> [16] "Gymnogobius urotaenia" #> [17] "Hemibarbus spp." #> [18] "Hypomesus nipponensis" #> [19] "Hypophthalmichthys spp." #> [20] "Hyporhamphus intermedius" #> [21] "Ictalurus punctatus" #> [22] "Ischikauia steenackeri" #> [23] "Lepomis macrochirus macrochirus" #> [24] "Leucopsarion petersii" #> [25] "Megalobrama amblycephala" #> [26] "Micropterus dolomieu dolomieu" #> [27] "Micropterus salmoides" #> [28] "Misgurnus spp." #> [29] "Monopterus albus" #> [30] "Mugil cephalus cephalus" #> [31] "Mylopharyngodon piceus" #> [32] "Nipponocypris sieboldii" #> [33] "Nipponocypris temminckii" #> [34] "Opsariichthys platypus" #> [35] "Opsariichthys uncirostris uncirostris" #> [36] "Oryzias latipes" #> [37] "Plecoglossus altivelis altivelis" #> [38] "Pseudogobio spp." #> [39] "Pseudorasbora parva" #> [40] "Rhinogobius spp." #> [41] "Rhodeus ocellatus ocellatus" #> [42] "Salangichthys microdon" #> [43] "Sarcocheilichthys variegatus microoculus" #> [44] "Silurus asotus" #> [45] "Squalidus chankaensis biwae" #> [46] "Tachysurus tokiensis" #> [47] "Tanakia lanceolata" #> [48] "Tribolodon brandtii maruta" #> [49] "Tribolodon hakonensis" #> [50] "Tridentiger spp." #> #> $Site #> [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12" "13" "14" "15" #> [16] "16" "17" "18" "19" "20" "21" "22" "23" "24" "25" "26" "27" "28" "29" "30" #> [31] "31" "32" "33" "34" "35" "36" "37" "38" "39" "40" "41" "42" "43" "44" "45" #> [46] "46" "47" "48" "49" "50" ## ----eval = FALSE------------------------------------------------------------- # post_summary_gamma <- get_post_summary(fit1, "gamma") # attributes(post_summary_gamma)$dimension ## ----------------------------------------------------------------------------- #> [1] "Species" "Effects" ## ----eval = FALSE------------------------------------------------------------- # attributes(post_summary_gamma)$label$Effects ## ----------------------------------------------------------------------------- #> [1] "(Intercept)" "riverbankwithout_vegetation" ## ----eval = FALSE------------------------------------------------------------- # post_summary_gamma ## ----------------------------------------------------------------------------- #> mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff overlap0 f #> (Omitted the beginning) #> gamma[1,2] -0.68135951 0.4102181 -1.4106885 -0.964385346 -0.71198673 -0.41908665 0.197784615 1.0031831 720 1 0.93875 #> gamma[2,2] -0.57124068 0.4861383 -1.4368696 -0.906331725 -0.60645174 -0.27472631 0.480190462 1.0141750 414 1 0.87925 #> gamma[3,2] -1.07479024 0.4260032 -1.9550037 -1.336688784 -1.07308990 -0.80888005 -0.213626578 1.0040679 632 0 0.99250 #> gamma[4,2] -0.96563859 0.5007208 -1.9760743 -1.276917710 -0.95208744 -0.64960250 0.016721077 1.0020286 1835 1 0.97225 #> gamma[5,2] -1.01038376 0.5436033 -2.1340409 -1.339282006 -0.98783592 -0.65193007 0.006887964 1.0148153 190 1 0.97400 #> gamma[6,2] -0.96443244 0.5694898 -2.1291123 -1.307706597 -0.96394933 -0.62099851 0.205483803 1.0043719 989 1 0.95325 #> gamma[7,2] -1.44522243 0.4356167 -2.4011743 -1.721131931 -1.39965212 -1.14153735 -0.689367497 1.0021543 1065 0 0.99925 #> gamma[8,2] -1.36146823 0.4820477 -2.3661876 -1.654265660 -1.32819584 -1.04709249 -0.483307820 1.0014346 1791 0 0.99900 #> (Omitted the remaining) ## ----eval = FALSE------------------------------------------------------------- # get_post_summary(fit1, "Mu") ## ----------------------------------------------------------------------------- #> mean sd 2.5% 25% 50% 75% #> Mu[1] -1.18914419 0.1466486 -1.4845230 -1.2872349 -1.18205680 -1.0888031 #> Mu[2] 0.94327206 0.1945952 0.5640255 0.8102745 0.94803155 1.0691758 #> Mu[3] 0.02522906 0.3264013 -0.6034793 -0.1859071 0.01929096 0.2317505 #> Mu[4] -1.03746412 0.1609537 -1.3626390 -1.1414288 -1.03743346 -0.9292544 #> 97.5% Rhat n.eff overlap0 f #> Mu[1] -0.9119219 1.108793 30 0 1.0000 #> Mu[2] 1.3284452 1.031154 100 0 1.0000 #> Mu[3] 0.6951928 1.008576 352 1 0.5255 #> Mu[4] -0.7312586 1.022904 135 0 1.0000 #> attr(,"dimension") #> [1] "Effects" #> attr(,"label") #> attr(,"label")$Effects #> [1] "phi | (Intercept)" "theta | (Intercept)" #> [3] "psi | (Intercept)" "psi | riverbankwithout_vegetation" ## ----eval = FALSE------------------------------------------------------------- # post_sample_psi <- get_post_samples(fit1, "psi") # dim(post_sample_psi) ## ----------------------------------------------------------------------------- #> [1] 4000 50 50 ## ----eval = FALSE------------------------------------------------------------- # attributes(post_sample_psi) ## ----------------------------------------------------------------------------- #> $dim #> [1] 4000 50 50 #> #> $dimension #> [1] "Sample" "Species" "Site" #> #> $label #> $label$Sample #> NULL #> #> $label$Species #> [1] "Abbottina rivularis" #> [2] "Acanthogobius lactipes" #> [3] "Acheilognathus macropterus" #> [4] "Acheilognathus rhombeus" #> [5] "Anguilla japonica" #> [6] "Biwia zezera" #> [7] "Carassius cuvieri" #> [8] "Carassius spp." #> [9] "Channa argus" #> [10] "Ctenopharyngodon idella" #> [11] "Cyprinus carpio" #> [12] "Gambusia affinis" #> [13] "Gnathopogon spp." #> [14] "Gymnogobius castaneus" #> [15] "Gymnogobius petschiliensis" #> [16] "Gymnogobius urotaenia" #> [17] "Hemibarbus spp." #> [18] "Hypomesus nipponensis" #> [19] "Hypophthalmichthys spp." #> [20] "Hyporhamphus intermedius" #> [21] "Ictalurus punctatus" #> [22] "Ischikauia steenackeri" #> [23] "Lepomis macrochirus macrochirus" #> [24] "Leucopsarion petersii" #> [25] "Megalobrama amblycephala" #> [26] "Micropterus dolomieu dolomieu" #> [27] "Micropterus salmoides" #> [28] "Misgurnus spp." #> [29] "Monopterus albus" #> [30] "Mugil cephalus cephalus" #> [31] "Mylopharyngodon piceus" #> [32] "Nipponocypris sieboldii" #> [33] "Nipponocypris temminckii" #> [34] "Opsariichthys platypus" #> [35] "Opsariichthys uncirostris uncirostris" #> [36] "Oryzias latipes" #> [37] "Plecoglossus altivelis altivelis" #> [38] "Pseudogobio spp." #> [39] "Pseudorasbora parva" #> [40] "Rhinogobius spp." #> [41] "Rhodeus ocellatus ocellatus" #> [42] "Salangichthys microdon" #> [43] "Sarcocheilichthys variegatus microoculus" #> [44] "Silurus asotus" #> [45] "Squalidus chankaensis biwae" #> [46] "Tachysurus tokiensis" #> [47] "Tanakia lanceolata" #> [48] "Tribolodon brandtii maruta" #> [49] "Tribolodon hakonensis" #> [50] "Tridentiger spp." #> #> $label$Site #> [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12" "13" "14" "15" #> [16] "16" "17" "18" "19" "20" "21" "22" "23" "24" "25" "26" "27" "28" "29" "30" #> [31] "31" "32" "33" "34" "35" "36" "37" "38" "39" "40" "41" "42" "43" "44" "45" #> [46] "46" "47" "48" "49" "50" ## ----eval = FALSE------------------------------------------------------------- # post_sample_psi[, 1, 1] ## ----------------------------------------------------------------------------- #> [1] 0.6673664 0.7505735 0.6416083 0.7011978 0.5935338 0.6087940 0.6159829 #> [8] 0.7758206 0.6710466 0.7423708 0.6244308 0.6237291 0.6344345 0.6715390 #> [15] 0.6255393 0.7306071 0.6781538 0.7391794 0.6678208 0.6715868 0.7582302 #> [22] 0.5947020 0.4880676 0.6958859 0.6202520 0.5487584 0.6781225 0.7073072 #> (Omitted the remaining) ## ----eval = FALSE------------------------------------------------------------- # predict_psi <- predict(fit1, parameter = "psi", type = "quantiles") # dim(predict_psi) ## ----------------------------------------------------------------------------- #> [1] 3 50 50 ## ----eval = FALSE------------------------------------------------------------- # attributes(predict_psi) ## ----------------------------------------------------------------------------- #> $dim #> [1] 3 50 50 #> #> $dimnames #> $dimnames[[1]] #> [1] "50%" "2.5%" "97.5%" #> #> $dimnames[[2]] #> NULL #> #> $dimnames[[3]] #> NULL #> #> #> $parameter #> [1] "psi" #> #> $scale #> [1] "response" #> #> $dimension #> [1] "Statistics" "Species" "Sites" #> #> $label #> $label$Statistics #> [1] "50%" "2.5%" "97.5%" #> #> $label$Species #> [1] "Abbottina rivularis" #> [2] "Acanthogobius lactipes" #> [3] "Acheilognathus macropterus" #> [4] "Acheilognathus rhombeus" #> [5] "Anguilla japonica" #> [6] "Biwia zezera" #> [7] "Carassius cuvieri" #> [8] "Carassius spp." #> [9] "Channa argus" #> [10] "Ctenopharyngodon idella" #> [11] "Cyprinus carpio" #> [12] "Gambusia affinis" #> [13] "Gnathopogon spp." #> [14] "Gymnogobius castaneus" #> [15] "Gymnogobius petschiliensis" #> [16] "Gymnogobius urotaenia" #> [17] "Hemibarbus spp." #> [18] "Hypomesus nipponensis" #> [19] "Hypophthalmichthys spp." #> [20] "Hyporhamphus intermedius" #> [21] "Ictalurus punctatus" #> [22] "Ischikauia steenackeri" #> [23] "Lepomis macrochirus macrochirus" #> [24] "Leucopsarion petersii" #> [25] "Megalobrama amblycephala" #> [26] "Micropterus dolomieu dolomieu" #> [27] "Micropterus salmoides" #> [28] "Misgurnus spp." #> [29] "Monopterus albus" #> [30] "Mugil cephalus cephalus" #> [31] "Mylopharyngodon piceus" #> [32] "Nipponocypris sieboldii" #> [33] "Nipponocypris temminckii" #> [34] "Opsariichthys platypus" #> [35] "Opsariichthys uncirostris uncirostris" #> [36] "Oryzias latipes" #> [37] "Plecoglossus altivelis altivelis" #> [38] "Pseudogobio spp." #> [39] "Pseudorasbora parva" #> [40] "Rhinogobius spp." #> [41] "Rhodeus ocellatus ocellatus" #> [42] "Salangichthys microdon" #> [43] "Sarcocheilichthys variegatus microoculus" #> [44] "Silurus asotus" #> [45] "Squalidus chankaensis biwae" #> [46] "Tachysurus tokiensis" #> [47] "Tanakia lanceolata" #> [48] "Tribolodon brandtii maruta" #> [49] "Tribolodon hakonensis" #> [50] "Tridentiger spp." #> #> $label$Sites #> [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12" "13" "14" "15" #> [16] "16" "17" "18" "19" "20" "21" "22" "23" "24" "25" "26" "27" "28" "29" "30" #> [31] "31" "32" "33" "34" "35" "36" "37" "38" "39" "40" "41" "42" "43" "44" "45" #> [46] "46" "47" "48" "49" "50" ## ----eval = FALSE------------------------------------------------------------- # new_y <- array(1, dim = c(50, 2, 1)) # new_riverbank <- factor(levels(fish_raw$riverbank)) # dimnames(new_y)[[1]] <- dimnames(fish_raw$y)[[1]] # dimnames(new_y)[[2]] <- levels(fish_raw$riverbank) # newdata <- occumbData( # y = new_y, # spec_cov = list(mismatch = fish_raw$mismatch), # site_cov = list(riverbank = new_riverbank) # ) ## ----eval = FALSE------------------------------------------------------------- # summary(newdata) ## ----------------------------------------------------------------------------- #> Sequence read counts: #> Number of species, I = 50 #> Number of sites, J = 2 #> Maximum number of replicates per site, K = 1 #> Number of missing observations = 0 #> Number of replicates per site: 1 (average), 0 (sd) #> Sequencing depth: 50 (average), 0 (sd) #> #> Species covariates: #> mismatch (continuous) #> Site covariates: #> riverbank (categorical) #> Replicate covariates: #> (None) #> #> Labels for species: #> Abbottina rivularis, Acanthogobius lactipes, Acheilognathus macropterus, Acheilognathus rhombeus, Anguilla japonica, Biwia zezera, Carassius cuvieri, Carassius spp., Channa argus, Ctenopharyngodon idella, Cyprinus carpio, Gambusia affinis, Gnathopogon spp., Gymnogobius castaneus, Gymnogobius petschiliensis, Gymnogobius urotaenia, Hemibarbus spp., Hypomesus nipponensis, Hypophthalmichthys spp., Hyporhamphus intermedius, Ictalurus punctatus, Ischikauia steenackeri, Lepomis macrochirus macrochirus, Leucopsarion petersii, Megalobrama amblycephala, Micropterus dolomieu dolomieu, Micropterus salmoides, Misgurnus spp., Monopterus albus, Mugil cephalus cephalus, Mylopharyngodon piceus, Nipponocypris sieboldii, Nipponocypris temminckii, Opsariichthys platypus, Opsariichthys uncirostris uncirostris, Oryzias latipes, Plecoglossus altivelis altivelis, Pseudogobio spp., Pseudorasbora parva, Rhinogobius spp., Rhodeus ocellatus ocellatus, Salangichthys microdon, Sarcocheilichthys variegatus microoculus, Silurus asotus, Squalidus chankaensis biwae, Tachysurus tokiensis, Tanakia lanceolata, Tribolodon brandtii maruta, Tribolodon hakonensis, Tridentiger spp. #> Labels for sites: #> with_vegetation, without_vegetation #> Labels for replicates: #> (None) ## ----eval = FALSE------------------------------------------------------------- # predict(fit1, newdata = newdata, parameter = "psi", type = "quantiles") ## ----------------------------------------------------------------------------- #> , , 1 #> #> [,1] [,2] [,3] [,4] [,5] [,6] [,7] #> 50% 0.6615133 0.2751231 0.3827295 0.15224665 0.17718177 0.09182312 0.8490621 #> 2.5% 0.4802434 0.1300821 0.2306797 0.05507663 0.06783645 0.02441254 0.7212899 #> 97.5% 0.7915599 0.4723714 0.5649520 0.32693805 0.38346795 0.24695832 0.9368627 #> [,8] [,9] [,10] [,11] [,12] [,13] [,14] #> 50% 0.9491461 0.6094722 0.4620570 0.9901477 0.08998769 0.8800989 0.15458012 #> 2.5% 0.8789348 0.4101236 0.2996333 0.9557121 0.02421783 0.7735010 0.05625485 #> 97.5% 0.9844834 0.8144525 0.6430711 0.9986149 0.23751947 0.9596028 0.31844988 #> [,15] [,16] [,17] [,18] [,19] [,20] [,21] #> 50% 0.12125019 0.9695970 0.9028737 0.15509196 0.5713331 0.2300419 0.8826944 #> 2.5% 0.03724413 0.9110775 0.8015870 0.05106154 0.3814354 0.1092197 0.7664569 #> 97.5% 0.26923379 0.9924366 0.9631803 0.34754250 0.7741282 0.4051336 0.9550034 #> [,22] [,23] [,24] [,25] [,26] [,27] [,28] #> 50% 0.7916404 0.7469745 0.06964821 0.15088247 0.10296476 0.8278437 0.9750893 #> 2.5% 0.6445037 0.5788305 0.01384897 0.04985234 0.02831969 0.7060083 0.9284540 #> 97.5% 0.9020019 0.8627692 0.19692563 0.34110667 0.25821933 0.9216507 0.9943395 #> [,29] [,30] [,31] [,32] [,33] [,34] [,35] #> 50% 0.06220973 0.4863689 0.07515230 0.4353972 0.2299977 0.8332137 0.5024979 #> 2.5% 0.01098996 0.2994199 0.01757379 0.2667359 0.1056246 0.6977411 0.3344813 #> 97.5% 0.19309849 0.6536841 0.21486234 0.6188191 0.4107509 0.9326166 0.6992168 #> [,36] [,37] [,38] [,39] [,40] [,41] [,42] #> 50% 0.5793313 0.09849281 0.5462387 0.9190994 0.9970209 0.8173368 0.09874564 #> 2.5% 0.4008518 0.02705771 0.3669093 0.8250025 0.9815437 0.6718077 0.02554623 #> 97.5% 0.7541354 0.25569027 0.7265439 0.9755158 0.9997452 0.9114977 0.26906733 #> [,43] [,44] [,45] [,46] [,47] [,48] [,49] #> 50% 0.15160773 0.5472967 0.3177699 0.1336518 0.16453410 0.3119199 0.1010661 #> 2.5% 0.05407187 0.3702903 0.1750710 0.0401887 0.06069858 0.1646373 0.0270361 #> 97.5% 0.32445164 0.7429405 0.4983112 0.3168030 0.32694679 0.4919569 0.2257502 #> [,50] #> 50% 0.9885246 #> 2.5% 0.9521853 #> 97.5% 0.9984315 #> #> , , 2 #> #> [,1] [,2] [,3] [,4] [,5] [,6] #> 50% 0.4756873 0.16831019 0.17344502 0.06827157 0.08065680 0.03799725 #> 2.5% 0.3210251 0.07468112 0.08268459 0.02152790 0.02558621 0.00904520 #> 97.5% 0.6477671 0.31803370 0.30375314 0.16653407 0.20628815 0.10914663 #> [,7] [,8] [,9] [,10] [,11] [,12] [,13] #> 50% 0.5747380 0.8277408 0.4451700 0.2415741 0.9426904 0.04637033 0.6435848 #> 2.5% 0.4085796 0.6832309 0.2727587 0.1336498 0.8361653 0.01254011 0.4719155 #> 97.5% 0.7316505 0.9241055 0.6653103 0.3976390 0.9855532 0.12900086 0.7834428 #> [,14] [,15] [,16] [,17] [,18] [,19] [,20] #> 50% 0.07037047 0.06163691 0.8602899 0.6909106 0.08289891 0.3366483 0.1156426 #> 2.5% 0.02364689 0.01866216 0.7192849 0.5273816 0.02726444 0.2012909 0.0498471 #> 97.5% 0.16850684 0.15171748 0.9468762 0.8219908 0.20203910 0.5227359 0.2284146 #> [,21] [,22] [,23] [,24] [,25] [,26] [,27] #> 50% 0.7204299 0.5831452 0.5653456 0.030925710 0.06781351 0.05881219 0.6000844 #> 2.5% 0.5574411 0.4218203 0.4080076 0.006236538 0.02063876 0.01733949 0.4410841 #> 97.5% 0.8503492 0.7422631 0.7276702 0.094870829 0.17049638 0.16207666 0.7470483 #> [,28] [,29] [,30] [,31] [,32] [,33] #> 50% 0.9126102 0.035225441 0.3407391 0.04236923 0.16616006 0.08886842 #> 2.5% 0.8020100 0.007672009 0.2031182 0.01018349 0.07272355 0.03188771 #> 97.5% 0.9696633 0.110552419 0.5178231 0.13816851 0.29161845 0.18440940 #> [,34] [,35] [,36] [,37] [,38] [,39] [,40] #> 50% 0.5159626 0.2344728 0.3705463 0.04820736 0.2161039 0.7693729 0.9831972 #> 2.5% 0.3512136 0.1224432 0.2295261 0.01206818 0.1039605 0.6147098 0.9306534 #> 97.5% 0.6751210 0.3846941 0.5527593 0.13468654 0.3594184 0.8923270 0.9976299 #> [,41] [,42] [,43] [,44] [,45] [,46] #> 50% 0.6527351 0.06006526 0.06753850 0.2371226 0.15535572 0.05145973 #> 2.5% 0.4967540 0.01537025 0.02157825 0.1134572 0.07282467 0.01382803 #> 97.5% 0.7956659 0.16542028 0.16515022 0.3848913 0.27405126 0.13368309 #> [,47] [,48] [,49] [,50] #> 50% 0.08990812 0.13653424 0.05038579 0.9428926 #> 2.5% 0.03330197 0.05777936 0.01466287 0.8401821 #> 97.5% 0.20500919 0.25327092 0.12890174 0.9873862 #> #> attr(,"parameter") #> [1] "psi" #> attr(,"scale") #> [1] "response" #> attr(,"dimension") #> [1] "Statistics" "Species" "Sites" #> attr(,"label") #> attr(,"label")$Statistics #> [1] "50%" "2.5%" "97.5%" #> #> attr(,"label")$Species #> [1] "Abbottina rivularis" #> [2] "Acanthogobius lactipes" #> [3] "Acheilognathus macropterus" #> [4] "Acheilognathus rhombeus" #> [5] "Anguilla japonica" #> [6] "Biwia zezera" #> [7] "Carassius cuvieri" #> [8] "Carassius spp." #> [9] "Channa argus" #> [10] "Ctenopharyngodon idella" #> [11] "Cyprinus carpio" #> [12] "Gambusia affinis" #> [13] "Gnathopogon spp." #> [14] "Gymnogobius castaneus" #> [15] "Gymnogobius petschiliensis" #> [16] "Gymnogobius urotaenia" #> [17] "Hemibarbus spp." #> [18] "Hypomesus nipponensis" #> [19] "Hypophthalmichthys spp." #> [20] "Hyporhamphus intermedius" #> [21] "Ictalurus punctatus" #> [22] "Ischikauia steenackeri" #> [23] "Lepomis macrochirus macrochirus" #> [24] "Leucopsarion petersii" #> [25] "Megalobrama amblycephala" #> [26] "Micropterus dolomieu dolomieu" #> [27] "Micropterus salmoides" #> [28] "Misgurnus spp." #> [29] "Monopterus albus" #> [30] "Mugil cephalus cephalus" #> [31] "Mylopharyngodon piceus" #> [32] "Nipponocypris sieboldii" #> [33] "Nipponocypris temminckii" #> [34] "Opsariichthys platypus" #> [35] "Opsariichthys uncirostris uncirostris" #> [36] "Oryzias latipes" #> [37] "Plecoglossus altivelis altivelis" #> [38] "Pseudogobio spp." #> [39] "Pseudorasbora parva" #> [40] "Rhinogobius spp." #> [41] "Rhodeus ocellatus ocellatus" #> [42] "Salangichthys microdon" #> [43] "Sarcocheilichthys variegatus microoculus" #> [44] "Silurus asotus" #> [45] "Squalidus chankaensis biwae" #> [46] "Tachysurus tokiensis" #> [47] "Tanakia lanceolata" #> [48] "Tribolodon brandtii maruta" #> [49] "Tribolodon hakonensis" #> [50] "Tridentiger spp." #> #> attr(,"label")$Sites #> [1] "with_vegetation" "without_vegetation" ## ----eval = FALSE------------------------------------------------------------- # gof_result <- gof(fit1, cores = 4) ## ----eval = FALSE------------------------------------------------------------- # gof_result ## ----------------------------------------------------------------------------- #> Posterior predictive check for an occumbFit object: #> Statistics: Freeman-Tukey #> p-value: 0.2005 #> Discrepancy statistics for observed data: 542.81 (mean), 17.92 (sd) #> Discrepancy statistics for replicated data: 523.28 (mean), 16.47 (sd) ## ----eval = FALSE------------------------------------------------------------- # plot(gof_result) ## ----eval = FALSE------------------------------------------------------------- # utilL1 <- eval_util_L(expand.grid(K = 1:3, N = c(1E3, 1E4, 1E5)), # fit1, # cores = 4 # ) # utilL1 ## ----------------------------------------------------------------------------- #> K N Utility #> 1 1 1e+03 13.03225 #> 2 2 1e+03 16.72177 #> 3 3 1e+03 18.21783 #> 4 1 1e+04 14.34826 #> 5 2 1e+04 17.72304 #> 6 3 1e+04 18.92009 #> 7 1 1e+05 15.04591 #> 8 2 1e+05 18.20935 #> 9 3 1e+05 19.23478 ## ----eval = FALSE------------------------------------------------------------- # settings <- list_cond_L( # budget = 875 * 1E3, # lambda1 = 0.01, # lambda2 = 5000, # fit1 # ) # utilL2 <- eval_util_L(settings, fit1, cores = 4) # utilL2 ## ----------------------------------------------------------------------------- #> budget lambda1 lambda2 K N Utility #> 1 875000 0.01 5000 1 1250000.00 15.45671 #> 2 875000 0.01 5000 2 375000.00 18.37862 #> 3 875000 0.01 5000 3 83333.33 19.21953 ## ----eval = FALSE------------------------------------------------------------- # utilR1 <- eval_util_R(expand.grid(J = 1:3, K = 1:3, N = c(1E3, 1E4, 1E5)), # fit1, # cores = 4 # ) # utilR1 ## ----------------------------------------------------------------------------- #> 1 1 1 1e+03 13.06227 #> 2 2 1 1e+03 18.77537 #> 3 3 1 1e+03 22.21901 #> 4 1 2 1e+03 16.75208 #> 5 2 2 1e+03 22.80510 #> 6 3 2 1e+03 26.33776 #> 7 1 3 1e+03 18.18526 #> 8 2 3 1e+03 24.63089 #> 9 3 3 1e+03 28.13433 #> 10 1 1 1e+04 14.36287 #> 11 2 1 1e+04 20.30642 #> 12 3 1 1e+04 23.85725 #> 13 1 2 1e+04 17.65318 #> 14 2 2 1e+04 23.98433 #> 15 3 2 1e+04 27.55905 #> 16 1 3 1e+04 19.03370 #> 17 2 3 1e+04 25.30884 #> 18 3 3 1e+04 29.14037 #> 19 1 1 1e+05 14.97290 #> 20 2 1 1e+05 21.16216 #> 21 3 1 1e+05 24.76543 #> 22 1 2 1e+05 18.25193 #> 23 2 2 1e+05 24.64478 #> 24 3 2 1e+05 28.20150 #> 25 1 3 1e+05 19.33757 #> 26 2 3 1e+05 25.85813 #> 27 3 3 1e+05 29.68874 ## ----eval = FALSE------------------------------------------------------------- # settings <- list_cond_R( # budget = 1125 * 1E3, # lambda1 = 0.01, # lambda2 = 5000, # lambda3 = 5000, # J = seq(5, 50, 5), # K = 1:4 # ) # utilR2 <- eval_util_R(settings, fit1, cores = 4) # utilR2 ## ----------------------------------------------------------------------------- #> budget lambda1 lambda2 lambda3 J K N Utility #> 1 1125000 0.01 5000 5000 5 1 21500000.00 30.04136 #> 2 1125000 0.01 5000 5000 10 1 10250000.00 35.64544 #> 3 1125000 0.01 5000 5000 15 1 6500000.00 38.60797 #> 4 1125000 0.01 5000 5000 20 1 4625000.00 40.60944 #> 5 1125000 0.01 5000 5000 25 1 3500000.00 42.02112 #> 6 1125000 0.01 5000 5000 30 1 2750000.00 43.15096 #> 7 1125000 0.01 5000 5000 35 1 2214285.71 43.89833 #> 8 1125000 0.01 5000 5000 40 1 1812500.00 44.71556 #> 9 1125000 0.01 5000 5000 45 1 1500000.00 45.22666 #> 10 1125000 0.01 5000 5000 50 1 1250000.00 45.71167 #> 11 1125000 0.01 5000 5000 5 2 10500000.00 33.34409 #> 12 1125000 0.01 5000 5000 10 2 4875000.00 38.81906 #> 13 1125000 0.01 5000 5000 15 2 3000000.00 41.58322 #> 14 1125000 0.01 5000 5000 20 2 2062500.00 43.47530 #> 15 1125000 0.01 5000 5000 25 2 1500000.00 44.69525 #> 16 1125000 0.01 5000 5000 30 2 1125000.00 45.58398 #> 17 1125000 0.01 5000 5000 35 2 857142.86 46.27036 #> 18 1125000 0.01 5000 5000 40 2 656250.00 46.82490 #> 19 1125000 0.01 5000 5000 45 2 500000.00 47.24841 #> 20 1125000 0.01 5000 5000 50 2 375000.00 47.56585 #> 21 1125000 0.01 5000 5000 5 3 6833333.33 34.55557 #> 22 1125000 0.01 5000 5000 10 3 3083333.33 40.03808 #> 23 1125000 0.01 5000 5000 15 3 1833333.33 42.85560 #> 24 1125000 0.01 5000 5000 20 3 1208333.33 44.51759 #> 25 1125000 0.01 5000 5000 25 3 833333.33 45.67757 #> 26 1125000 0.01 5000 5000 30 3 583333.33 46.48685 #> 27 1125000 0.01 5000 5000 35 3 404761.90 47.14072 #> 28 1125000 0.01 5000 5000 40 3 270833.33 47.60570 #> 29 1125000 0.01 5000 5000 45 3 166666.67 47.93433 #> 30 1125000 0.01 5000 5000 50 3 83333.33 48.12801 #> 31 1125000 0.01 5000 5000 5 4 5000000.00 35.17767 #> 32 1125000 0.01 5000 5000 10 4 2187500.00 40.72769 #> 33 1125000 0.01 5000 5000 15 4 1250000.00 43.51095 #> 34 1125000 0.01 5000 5000 20 4 781250.00 45.16898 #> 35 1125000 0.01 5000 5000 25 4 500000.00 46.23139 #> 36 1125000 0.01 5000 5000 30 4 312500.00 46.98586 #> 37 1125000 0.01 5000 5000 35 4 178571.43 47.53393 #> 38 1125000 0.01 5000 5000 40 4 78125.00 47.89744