Tables
Table 1: Poisson regression model output: summary statistics of the posterior samples of diversity measurements. The posterior means and medians are used as the estimated diversity measurements and the 95% credible intervals (95% CI) give ranges of the diversity measurements associated with 0.95 inclusion probability.The diversity measurements computed based on the observed counts of different species of butterflies in different locations is also listed to compare with the model estimated diversity measurements. The indicator variable “violation” shows whether the credible interval covers the observed value, 1 stands for failed to cover and 0 means that the observed value is covered. The potential scale reproduction factor, R-hat, measures the convergence of the posterior samples. And R-hat=1.0 indicates that convergence is achieved and the resulting Bayesian estimates are reliable. βA is the additive model of species richness and βM is the multiplicative model.
MEAN / MEDIAN / 95 % CI / OBSERVED / VIOLATION / R-hatαpatch / 19.43 / 19.42 / (18.12, 20.77) / 16.69 / 1 / 1.0
αcluster / 33.56 / 33.50 / (31.17, 36.00) / 31.67 / 0 / 1.0
γ / 47.13 / 47.00 / (44.00, 49.00) / 49.00 / 0 / 1.0
βApatch / 14.14 / 14.13 / (12.35, 15.96) / 14.93 / 0 / 1.0
βAcluster / 13.57 / 13.67 / (10.83, 16.17) / 17.33 / 1 / 1.0
βMpatch / 1.73 / 1.73 / (1.63, 1.83) / 1.90 / 1 / 1.0
βMcluster / 1.41 / 1.41 / (1.31, 1.51) / 1.55 / 1 / 1.0
DIC = 2840.750
Table 2: ZIP model (species-specific priors for the mixture proportion) output: summary statistics of the posterior samples of diversity measurements. The posterior means and medians are both considered as the estimated diversity measurements and the 95% credible intervals (95% CI) give ranges of the diversity measurements associated with 0.95 inclusion probability.The diversity measurements computed based on the observed counts of different species of butterflies in different locations is also listed to compare with the model estimated diversity measurements. The indicator variable “violation” shows whether the credible interval covers the observed value, 1 stands for failed to cover and 0 means that the observed value is covered. The potential scale reproduction factor, R-hat, measures the convergence of the posterior samples. And R-hat=1.0 indicates that convergence is achieved and the resulting Bayesian estimates are reliable. βA is the additive model of species richness and βM is the multiplicative model.
MEAN / MEDIAN / 2.5% / 97.5% / OBSERVED / VIOLATION / R-hatαcluster / 31.59 / 31.67 / 29.50 / 33.83 / 31.67 / 0 / 1.0
αpatch / 16.76 / 16.77 / 15.69 / 17.85 / 16.69 / 0 / 1.0
γ / 46.87 / 47.00 / 44.00 / 49.00 / 49.00 / 0 / 1.0
βApatch / 14.83 / 14.83 / 13.17 / 16.47 / 14.93 / 0 / 1.0
βAcluster / 15.29 / 15.33 / 12.67 / 17.83 / 17.33 / 0 / 1.0
Mpatch / 1.89 / 1.88 / 1.79 / 1.99 / 1.90 / 0 / 1.0
βMcluster / 1.48 / 1.48 / 1.38 / 1.59 / 1.55 / 0 / 1.0
DIC = 2940.5
Table 3: ZIP model with environmental predictors (species-specific priors for the mixture proportion) output: summary statistics of the posterior samples of diversity measurements. The posterior means and medians are both considered as the estimated diversity measurements and the 95% credible intervals (95% CI) give ranges of the diversity measurements associated with 0.95 inclusion probability.The diversity measurements computed based on the observed counts of different species of butterflies in different locations is also listed to compare with the model estimated diversity measurements. The indicator variable “violation” shows whether the credible interval covers the observed value, 1 stands for failed to cover and 0 means that the observed value is covered. The potential scale reproduction factor, R-hat, measures the convergence of the posterior samples. And R-hat=1.0 indicates that convergence is achieved and the resulting Bayesian estimates are reliable. βA is the additive model of species richness and βM is the multiplicative model.
MEAN / MEDIAN / 2.5% / 97.5% / OBSERVED / VIOLATION / R-hatαcluster / 31.86 / 31.83 / 29.67 / 34.00 / 31.67 / 0 / 1.0
αpatch / 16.89 / 16.88 / 15.85 / 18.00 / 16.69 / 0 / 1.0
γ / 46.99 / 47.00 / 44.00 / 49.00 / 49.00 / 0 / 1.0
βApatch / 14.97 / 14.97 / 13.29 / 16.62 / 14.93 / 0 / 1.0
βAcluster / 15.13 / 15.17 / 12.50 / 17.67 / 17.33 / 0 / 1.0
Mpatch / 1.89 / 1.89 / 1.79 / 1.99 / 1.90 / 0 / 1.0
βMcluster / 1.48 / 1.48 / 1.38 / 1.58 / 1.55 / 0 / 1.0
DIC = 3078.1
Table 4: ZIP model with environmental predictors (species-specific priors for the mixture proportion) output: summary statistics of the posterior samples of regression coefficients
MEAN / MEDIAN / SD / 2.5% / 97.5%b1 / 0.18 / 0.18 / 0.04 / 0.10 / 0.26
b2 / 0.17 / 0.17 / 0.08 / 0.01 / 0.32
b3 / -0.05 / -0.04 / 0.07 / -0.20 / 0.08
b4 / 0.02 / 0.02 / 0.01 / 0.01 / 0.04
Code
------
# ZIP model without Environmental information:
library("R2WinBUGS")
########################################################################
#######################data set########################################
# “y” is the count of butterflies, “s” represents the species, “c” represents the cluster, “p”
# represents the patch.
#########################################################################
y <-c(
0,0,0,0,0,0,0,1,0,1,0,0,0,1,0,0,1,1,2,0,1,1,0,1,5,0,
0,0,0,0,0,0,0,1,1,0,0,0,0,2,0,0,0,0,0,1,0,0,0,0,0,0,
0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,1,3,
2,2,1,7,2,2,5,8,4,0,1,0,0,0,0,2,3,0,4,7,1,3,0,1,13,0,
4,3,0,2,1,3,2,4,3,2,4,2,1,2,7,6,1,3,9,4,3,11,0,3,3,1,
3,1,0,5,0,0,5,1,0,1,2,0,0,1,15,3,0,5,10,5,1,6,1,1,2,2,
6,2,3,3,2,0,1,2,2,2,2,0,2,0,7,3,1,0,0,1,0,2,1,4,1,1,
7,0,0,5,1,0,4,0,5,3,1,0,1,0,2,0,0,0,0,0,0,3,0,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,1,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,1,
1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
0,0,0,0,0,0,1,0,0,0,1,0,0,1,1,0,0,0,1,0,0,2,0,0,0,3,
1,0,2,1,0,1,1,0,5,0,2,0,1,0,1,1,0,0,0,0,0,4,0,1,0,0,
0,6,2,2,3,0,3,0,12,1,1,2,8,1,5,2,12,7,4,2,12,24,3,5,0,15,
0,0,0,0,1,0,0,0,0,0,0,1,0,1,1,0,4,1,1,1,3,0,0,0,1,0,
0,4,0,0,0,0,6,1,4,2,1,0,0,2,4,1,1,1,2,1,4,3,0,2,1,1,
0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,3,0,5,0,1,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,3,0,0,0,0,
0,1,2,0,0,0,1,1,1,1,2,0,0,0,0,0,2,1,0,1,1,0,0,1,0,2,
23,4,6,10,9,8,15,6,15,2,23,8,1,17,13,12,11,9,3,3,9,22,13,13,2,12,
2,0,2,1,1,1,10,3,0,1,2,0,1,6,11,16,3,7,13,14,6,16,5,5,14,22,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,0,4,0,0,0,
0,0,0,1,0,0,0,0,0,0,2,0,0,1,0,0,0,0,0,0,0,0,0,0,0,2,
0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,2,
2,1,0,1,1,2,3,0,1,0,0,0,0,1,0,0,0,0,1,4,0,1,0,0,0,2,
1,1,0,0,0,0,2,0,1,0,0,0,0,0,1,0,1,0,0,2,0,1,0,0,0,0,
4,1,1,4,0,0,1,3,1,2,0,0,2,0,1,0,1,0,2,2,1,7,1,3,0,0,
0,1,3,4,1,0,0,0,0,0,0,0,0,0,2,0,3,2,1,2,3,3,0,3,3,0,
0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,
0,1,0,0,2,0,1,0,1,0,0,0,1,1,6,0,0,0,3,0,0,0,0,1,0,0,
0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,3,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,
1,0,1,0,0,0,2,2,0,0,0,0,0,0,0,1,1,1,0,0,0,1,1,1,2,0,
0,0,0,0,0,0,0,1,1,0,2,0,1,1,0,0,1,0,0,0,0,0,1,0,0,0,
0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
2,0,0,0,0,1,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,
4,0,1,0,0,0,0,0,3,2,0,0,0,0,0,2,0,2,2,0,0,1,0,1,1,2,
4,0,0,2,0,0,0,3,1,1,1,0,1,2,1,0,0,0,0,0,0,0,0,3,1,0,
0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,
0,0,1,0,0,0,0,1,1,1,0,0,1,0,0,0,0,0,0,0,0,1,0,1,0,0,
0,0,1,0,3,0,0,0,2,0,2,2,4,4,3,1,2,1,2,1,1,3,2,2,0,2,
3,1,0,0,1,1,0,0,0,0,1,1,0,2,0,0,0,0,3,0,0,2,0,6,0,7,
0,0,0,1,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,
0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
3,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,0,3,0,0,1,0,4,1,0)
s=c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,
3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,
4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,
5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,
6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,
9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9,
10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,
11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,
12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,12,
13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,
14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,
15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,
16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,16,
17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,17,
18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,
19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,
20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,20,
21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,21,
22,22,22,22,22,22,22,22,22,22,22,22,22,22,22,22,22,22,22,22,22,22,22,22,22,22,
23,23,23,23,23,23,23,23,23,23,23,23,23,23,23,23,23,23,23,23,23,23,23,23,23,23,
24,24,24,24,24,24,24,24,24,24,24,24,24,24,24,24,24,24,24,24,24,24,24,24,24,24,
25,25,25,25,25,25,25,25,25,25,25,25,25,25,25,25,25,25,25,25,25,25,25,25,25,25,
26,26,26,26,26,26,26,26,26,26,26,26,26,26,26,26,26,26,26,26,26,26,26,26,26,26,
27,27,27,27,27,27,27,27,27,27,27,27,27,27,27,27,27,27,27,27,27,27,27,27,27,27,
28,28,28,28,28,28,28,28,28,28,28,28,28,28,28,28,28,28,28,28,28,28,28,28,28,28,
29,29,29,29,29,29,29,29,29,29,29,29,29,29,29,29,29,29,29,29,29,29,29,29,29,29,
30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,30,
31,31,31,31,31,31,31,31,31,31,31,31,31,31,31,31,31,31,31,31,31,31,31,31,31,31,
32,32,32,32,32,32,32,32,32,32,32,32,32,32,32,32,32,32,32,32,32,32,32,32,32,32,
33,33,33,33,33,33,33,33,33,33,33,33,33,33,33,33,33,33,33,33,33,33,33,33,33,33,
34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,
35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,35,
36,36,36,36,36,36,36,36,36,36,36,36,36,36,36,36,36,36,36,36,36,36,36,36,36,36,
37,37,37,37,37,37,37,37,37,37,37,37,37,37,37,37,37,37,37,37,37,37,37,37,37,37,
38,38,38,38,38,38,38,38,38,38,38,38,38,38,38,38,38,38,38,38,38,38,38,38,38,38,
39,39,39,39,39,39,39,39,39,39,39,39,39,39,39,39,39,39,39,39,39,39,39,39,39,39,
40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,40,
41,41,41,41,41,41,41,41,41,41,41,41,41,41,41,41,41,41,41,41,41,41,41,41,41,41,
42,42,42,42,42,42,42,42,42,42,42,42,42,42,42,42,42,42,42,42,42,42,42,42,42,42,
43,43,43,43,43,43,43,43,43,43,43,43,43,43,43,43,43,43,43,43,43,43,43,43,43,43,
44,44,44,44,44,44,44,44,44,44,44,44,44,44,44,44,44,44,44,44,44,44,44,44,44,44,
45,45,45,45,45,45,45,45,45,45,45,45,45,45,45,45,45,45,45,45,45,45,45,45,45,45,
46,46,46,46,46,46,46,46,46,46,46,46,46,46,46,46,46,46,46,46,46,46,46,46,46,46,
47,47,47,47,47,47,47,47,47,47,47,47,47,47,47,47,47,47,47,47,47,47,47,47,47,47,
48,48,48,48,48,48,48,48,48,48,48,48,48,48,48,48,48,48,48,48,48,48,48,48,48,48,
49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49,49)
c=c(
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6,
1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,5,5,5,5,5,6,6,6,6,6)
p=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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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)
###################Load the saved workspace file#################################
setwd("C:/Jing_miami/PJHou/paper_draft_Crist/code_07032012/")
modfile <- "bugs_code.txt"
bugsdir <- "C:/WinBUGS14"
############Write the Winbugs model file directly in R#############################
mod <- function()
{
mu ~ dnorm(0.0,1.0)
for(i in 1:49) {
species[i] ~ dnorm(0.0,tau.s)
}
for(i in 1:6) {
cluster[i] ~ dnorm(0.0,tau.c)
}
for(i in 1:26) {
patch[i] ~ dnorm(0.0,tau.p)
}
sigma.s ~ dunif(0,5)
sigma.c ~ dunif(0,5)
sigma.p ~ dunif(0,5)
tau.s <- pow(sigma.s,-2)
tau.c <- pow(sigma.c,-2)
tau.p <- pow(sigma.p,-2)
for(i in 1:1274) {
y[i] ~ dpois(m[i])
m[i]<-(1-eta[i])*lambda[i]
eta[i] ~dbern(p0[s[i]])
log(lambda[i]) <- mu+species[s[i]]+cluster[c[i]]+patch[p[i]]
fd[i]<- exp(-lambda[i]+y[i]*log(lambda[i])-loggam(y[i]+1))
l[i] <- log(p0[s[i]]*equals(y[i],0)+(1-p0[s[i]])*fd[i])
y.pred[i]~dpois(m[i])
}
for (j in 1:49)
{
p0[j]~dbeta(1,1)
}
# compute alpha and beta(cluster) and beta(patch) here
for(i in 1:26){
for(j in 1:49){
yp_matrix[i,j]<-y.pred[(j-1)*26+i]
}
}
totabun<-sum(yp_matrix[,])
for(i in 1:49){
s_sum[i]<-sum(yp_matrix[,i])
ind1[i]<-step(s_sum[i]-0.5)
}
gamma<-sum(ind1[1:49])
# calculate richness by patch and cluster
for(i in 1:26){ for(j in 1:49){
ind2[i,j]<-step(yp_matrix[i,j]-0.5)
}
patch.rich[i]<-sum(ind2[i,])
}
for(j in 1:49){
cluster.mat[1,j]<-sum(yp_matrix[1:4,j])
cluster.mat[2,j]<-sum(yp_matrix[5:7,j])
cluster.mat[3,j]<-sum(yp_matrix[8:10,j])
cluster.mat[4,j]<-sum(yp_matrix[11:13,j])
cluster.mat[5,j]<-sum(yp_matrix[14:21,j])
cluster.mat[6,j]<-sum(yp_matrix[22:26,j])
}
for(i in 1:6){ for(j in 1:49){
ind3[i,j]<-step(cluster.mat[i,j]-0.5)
}
crich[i]<-sum(ind3[i,])
cabun[i]<-sum(cluster.mat[i,])
}
# calculate additive alphas and betas for patch and cluster
alpha.patch<-sum(patch.rich[1:26])/26
alpha.cluster<-sum(crich[1:6])/6
beta.patch<-alpha.cluster-alpha.patch
beta.cluster<-gamma-alpha.cluster
# calculate multiplicative betas
beta.patch.mult<-alpha.cluster/alpha.patch
beta.cluster.mult<-gamma/alpha.cluster
}
#############################################################################
write.model(mod,modfile)
#############################################################################
############################################################################
###########R+WINBUGS#########################################################
#############################################################################
library("R2WinBUGS")
butterfly.data <-list("y","c","s","p")
butterfly.parameters <- c("gamma" , "alpha.cluster" , "beta.cluster" , "alpha.patch" , "beta.patch" , "beta.patch.mult" , "beta.cluster.mult" , "y.pred","p0")
inits <- NULL
butterfly.sim <-bugs(butterfly.data,inits,butterfly.parameters,modfile,n.chains=3,n.iter=40000,n.burnin=10000,n.thin=5,bugs.directory="C:/WinBUGS14/")
#############################################################################
alpha.cluster<-butterfly.sim$sims.list$alpha.cluster
gamma<-butterfly.sim$sims.list$gamma
beta.cluster<-butterfly.sim$sims.list$beta.cluster
alpha.patch<-butterfly.sim$sims.list$alpha.patch
beta.patch<-butterfly.sim$sims.list$beta.patch
beta.patch.mult<-butterfly.sim$sims.list$beta.patch.mult
beta.cluster.mult<-butterfly.sim$sims.list$beta.cluster.mult
p0<-butterfly.sim$sims.list$p0
############################################################################
p0m<-matrix(,49,6)
for(i in 1:49){
p0m[i,1] <- mean(p0[,i])
p0m[i,2] <- median(p0[,i])
p0m[i,3]<-quantile(p0[,i],0.025)
p0m[i,4]<-quantile(p0[,i],0.975)
}
plot(p0m[,2],ylim=c(0,1),xlab='species: 1 to 49',ylab='Posterior quantiles of the mixture proportion')
for(i in 1:49){
xx<-c(i-0.05,i+0.05)
yy<-c(p0m[i,3],p0m[i,3])
yy2<-c(p0m[i,4],p0m[i,4])
zz<-c(i-0.05,i-0.05)
zz2<-c(i+0.05,i+0.05)
ss<-c(p0m[i,3],p0m[i,4])
lines(xx,yy,col='blue')
lines(xx,yy2,col='blue')
lines(zz,ss,col='blue')
lines(zz2,ss,col='blue')
}
summary<-matrix(,7,6)
summary[1,1] <- mean(alpha.cluster)
summary[1,2]<-median(alpha.cluster)
summary[1,3]<-quantile(alpha.cluster,0.025)
summary[1,4]<-quantile(alpha.cluster,0.975)
summary[1,5]<-31.67
summary[2,1] <- mean(alpha.patch)
summary[2,2]<-median(alpha.patch)
summary[2,3]<-quantile(alpha.patch,0.025)
summary[2,4]<-quantile(alpha.patch,0.975)
summary[2,5]<-16.69
summary[3,1] <- mean(beta.cluster)
summary[3,2]<-median(beta.cluster)
summary[3,3]<-quantile(beta.cluster,0.025)
summary[3,4]<-quantile(beta.cluster,0.975)
summary[3,5]<-17.33
summary[4,1] <- mean(beta.patch)
summary[4,2]<-median(beta.patch)
summary[4,3]<-quantile(beta.patch,0.025)
summary[4,4]<-quantile(beta.patch,0.975)
summary[4,5]<-14.93
summary[5,1] <- mean(gamma)
summary[5,2]<-median(gamma)
summary[5,3]<-quantile(gamma,0.025)
summary[5,4]<-quantile(gamma,0.975)
summary[5,5]<-49
summary[6,1] <- mean(beta.cluster.mult)
summary[6,2]<-median(beta.cluster.mult)
summary[6,3]<-quantile(beta.cluster.mult,0.025)
summary[6,4]<-quantile(beta.cluster.mult,0.975)
summary[6,5]<-1.55
summary[7,1] <- mean(beta.patch.mult)
summary[7,2]<-median(beta.patch.mult)
summary[7,3]<-quantile(beta.patch.mult,0.025)
summary[7,4]<-quantile(beta.patch.mult,0.975)
summary[7,5]<-1.9
for (i in 1:7)
{
if ((summary[i,5]>summary[i,4])||(summary[i,5]<summary[i,3])) (summary[i,6]<-1)
else (summary[i,6]<-0)
}
rownames(summary) <- c("alpha.cluster","alpha.patch","beta.cluster","beta.patch","gamma","beta.cluster.mult","beta.patch.mult")
colnames(summary) <- c("mean","median","2.5%","97.5","observed","violation")
summary