I have recently started working on a dataset, but my limited experience as a recent graduate is preventing me from making progress. After performing a rarefaction (a form of bootstrapping), which was necessary for the project, I now have a dataframe with 24 species diversity indices, each accompanied by a confidence interval generated from the random rarefaction process.
These 24 diversity indices belong to 2 different non-independent groups that I need to compare. If there were no confidence intervals and I could work solely with the means, I would use a paired t-test. I have 12 plots (and so 12 rarefied diversity indices) with the treatment 1 and 12 non-independent (spatially correlated) plots with the treatment 2. I want to test if diversity is different along treatment. I would treat it like a t-test but rarefaction/bootstrapping gives me a confidence interval alongside my indices.
To account for the confidence intervals obtained through rarefaction, I tried performing a weighted t-test using the variance derived from the interval. However, I am not confident in my method and cannot find any articles online that support this approach.
This is my dataframe:
size1 <- c(32.67256,30.59280,33.56214,30.15552,29.02073,24.92427,34.79967,34.26559,29.33457,25.75716,27.91638,33.87884)
size2 <-c(34.18847,30.94369,37.38462,22.96785,27.98805,31.39834,30.26401,33.59788,36.19856,31.19667,21.27245,28.87137)
size <-c(size1, size2)
size1CI05 <-c(29.50177,26.23491,29.60487,25.94388,24.48941,21.78924,29.24082,28.09738,25.20056,21.21051,24.40705,30.08490)
size2CI05 <-c(29.51654,23.75201,31.29203,17.60071,21.92296,26.38236,23.69115,26.15880,26.44932,26.72620,17.48761,23.76434)
sizeCI05 <-c(size1CI05,size2CI05)
size1CI95 <-c(35.84336,34.95070,37.51941,34.36716,33.55205,28.05929,40.35851,40.43380,33.46857,30.30381,31.42571,37.67278)
size2CI95 <-c(38.86039,38.13538,43.47722,28.33500,34.05315,36.41432,36.83687,41.03696,45.94780,35.66713,25.05730,33.97840)
sizeCI95 <-c(size1CI95,size2CI95)
group <- c(c(rep("1", 12), rep("2", 12))
df <- data.frame(size, sizeCI05, sizeCI95, group)
I also add the dput output as asked:
df <-structure(list(size = c(32.67256, 30.5928, 33.56214, 30.15552,
29.02073, 24.92427, 34.79967, 34.26559, 29.33457, 25.75716, 27.91638,
33.87884, 34.18847, 30.94369, 37.38462, 22.96785, 27.98805, 31.39834,
30.26401, 33.59788, 36.19856, 31.19667, 21.27245, 28.87137),
size05CI = c(29.5017662019491, 26.2349058385758, 29.6048735822698,
25.9438818737996, 24.4894140422372, 21.7892435074709, 29.2408246478003,
28.097376313181, 25.2005633026381, 21.2105115402011, 24.4070457401779,
30.0849045894848, 29.5165407340379, 23.7520071901393, 31.2920265619889,
17.6007091843407, 21.9229562124319, 26.3823600157115, 23.6911524798744,
26.158796800294, 26.4493176960902, 26.7262037524753, 17.4876088906701,
23.7643440166671), size95CI = c(35.8433590913617, 34.9506991382152,
37.5194133001315, 34.3671592174968, 33.5520534941053, 28.0592932279925,
40.3585126169722, 40.4338032764411, 33.4685713521404, 30.303809204322,
31.4257104335821, 37.6727751129716, 38.8603902871733, 38.1353758799793,
43.4772171945701, 28.3349990423006, 34.0531465846994, 36.4143170389483,
36.8368654619253, 41.0369594182707, 45.9478039221974, 35.6671344715295,
25.0572999101401, 33.9783955274194), group = c("1", "1",
"1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "2", "2",
"2", "2", "2", "2", "2", "2", "2", "2", "2", "2")), class = "data.frame", row.names = c(NA,
-24L))
This is what I tried to do:
install.packages("weights", dependencies=TRUE)
library(weights)
SE1 <- (size1CI95-size1CI05)/3.919928
SE2 <- (size2CI95-size2CI05)/3.919928
sizepaired <- size1-size2
combSE <- sqrt(SE1^2+SE2^2)
weights <- 1/combSE^2
wtd.t.test(x=sizepaired, y=0, weight=weights)
My approach was to take into account for the confidence interval (i.e., the standard deviation) originated due to rarefaction in my test by using a weighted test, where the weights are derived from the standard deviation. However, I am uncertain whether this is a valid approach, and I have not found any scientific articles addressing this method. My colleagues are somewhat skeptical and would prefer not to use 'experimental' methods.
Moreover, I am currently using a parametric test, but among my data (which I have not included here), some do not follow a normal distribution. Thus, I should consider how to conduct a weighted paired non-parametric test
As asked I'm editing with my raw data (how many individuals in each plot):
df <-structure(list(Group1_RH = c(16L, 0L, 1L, 0L, 4L, 8L, 0L, 1L, 6L,
0L, 4L, 2L, 0L, 1L, 0L, 0L, 1L, 5L, 1L, 6L, 0L, 4L, 2L, 0L, 0L,
2L, 0L, 0L, 3L, 14L, 16L, 2L, 1L, 0L, 4L, 1L, 3L, 8L, 2L, 0L,
6L, 15L, 3L, 6L, 0L, 0L, 0L, 4L, 2L, 0L, 5L, 0L, 0L, 0L, 0L,
0L, 9L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 3L, 1L, 1L, 0L, 0L, 0L),
Group2_RH = c(5L, 2L, 0L, 0L, 0L, 5L, 0L, 0L, 5L, 0L, 0L,
1L, 0L, 1L, 0L, 1L, 5L, 5L, 1L, 12L, 0L, 1L, 2L, 0L, 0L,
0L, 1L, 0L, 1L, 8L, 8L, 0L, 0L, 0L, 2L, 0L, 2L, 6L, 3L, 0L,
5L, 0L, 0L, 2L, 1L, 2L, 1L, 4L, 1L, 0L, 1L, 4L, 0L, 0L, 1L,
0L, 0L, 0L, 0L, 3L, 0L, 0L, 0L, 0L, 2L, 1L, 1L, 0L, 0L, 0L
), Group1_VA = c(9L, 0L, 1L, 0L, 31L, 1L, 0L, 1L, 1L, 0L,
0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 6L, 0L, 4L, 2L, 1L, 1L,
0L, 0L, 0L, 4L, 13L, 7L, 3L, 5L, 0L, 1L, 2L, 0L, 1L, 1L,
0L, 5L, 0L, 15L, 2L, 0L, 0L, 0L, 1L, 1L, 2L, 2L, 0L, 0L,
1L, 5L, 0L, 0L, 0L, 2L, 4L, 0L, 0L, 1L, 0L, 4L, 1L, 1L, 0L,
0L, 0L), Group2_VA = c(3L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 2L, 5L, 2L, 2L,
0L, 0L, 0L, 0L, 0L, 1L, 8L, 3L, 0L, 0L, 0L, 0L, 2L, 0L, 0L,
1L, 0L, 7L, 1L, 0L, 1L, 0L, 0L, 0L, 5L, 1L, 1L, 1L, 0L, 0L,
2L, 4L, 0L, 0L, 1L, 0L, 2L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L,
0L, 0L), Group1_FO = c(6L, 0L, 3L, 0L, 20L, 2L, 0L, 1L, 1L,
0L, 0L, 1L, 0L, 0L, 3L, 0L, 1L, 10L, 9L, 2L, 0L, 5L, 2L,
0L, 0L, 1L, 1L, 0L, 1L, 10L, 7L, 0L, 1L, 0L, 22L, 1L, 0L,
1L, 2L, 4L, 5L, 0L, 1L, 5L, 0L, 0L, 0L, 4L, 2L, 1L, 2L, 0L,
2L, 0L, 3L, 0L, 0L, 0L, 0L, 2L, 0L, 1L, 1L, 0L, 4L, 1L, 2L,
0L, 0L, 0L), Group2_FO = c(1L, 0L, 1L, 0L, 0L, 2L, 0L,
0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 2L, 3L, 4L, 0L, 1L,
1L, 0L, 0L, 0L, 1L, 0L, 1L, 13L, 14L, 3L, 0L, 0L, 1L, 1L,
0L, 2L, 5L, 0L, 4L, 1L, 0L, 6L, 0L, 2L, 1L, 3L, 1L, 0L, 1L,
5L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 3L, 2L, 1L, 1L, 0L, 3L, 1L,
3L, 1L, 0L, 0L), Group1_AI = c(6L, 0L, 0L, 0L, 10L, 1L,
0L, 0L, 1L, 0L, 0L, 0L, 3L, 0L, 0L, 1L, 0L, 0L, 0L, 6L, 0L,
2L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 14L, 4L, 1L, 0L, 0L, 1L,
1L, 0L, 0L, 2L, 0L, 7L, 0L, 0L, 5L, 0L, 0L, 5L, 4L, 2L, 1L,
2L, 0L, 0L, 3L, 8L, 0L, 1L, 1L, 0L, 6L, 0L, 0L, 2L, 0L, 4L,
2L, 1L, 0L, 0L, 0L), Group2_AI = c(1L, 0L, 0L, 0L, 0L,
2L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 3L,
0L, 2L, 2L, 0L, 0L, 0L, 0L, 0L, 0L, 10L, 8L, 0L, 0L, 0L,
0L, 0L, 1L, 2L, 0L, 0L, 5L, 0L, 0L, 0L, 0L, 0L, 0L, 3L, 1L,
2L, 1L, 0L, 0L, 1L, 7L, 0L, 0L, 0L, 0L, 3L, 0L, 0L, 1L, 0L,
4L, 0L, 1L, 0L, 0L, 0L), Group1_AZ = c(8L, 0L, 0L, 0L,
6L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 2L, 1L, 0L, 0L, 0L, 0L, 0L,
2L, 0L, 3L, 2L, 1L, 2L, 0L, 0L, 0L, 0L, 13L, 6L, 1L, 0L,
0L, 1L, 2L, 0L, 1L, 1L, 0L, 8L, 0L, 0L, 0L, 0L, 0L, 0L, 1L,
3L, 2L, 0L, 0L, 0L, 0L, 6L, 0L, 0L, 0L, 0L, 5L, 0L, 0L, 1L,
0L, 5L, 3L, 2L, 1L, 0L, 0L), Group2_AZ = c(6L, 0L, 0L,
0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 4L, 0L, 0L, 0L, 0L, 0L,
0L, 2L, 0L, 2L, 3L, 0L, 0L, 0L, 0L, 0L, 0L, 11L, 4L, 0L,
0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 8L, 0L, 0L, 1L, 0L, 0L, 0L,
1L, 1L, 3L, 0L, 0L, 0L, 2L, 4L, 0L, 0L, 0L, 0L, 3L, 0L, 1L,
1L, 0L, 2L, 1L, 2L, 1L, 0L, 0L), Group1_LU = c(3L, 0L, 2L,
8L, 11L, 5L, 1L, 0L, 6L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 2L,
0L, 0L, 13L, 0L, 2L, 1L, 0L, 0L, 0L, 0L, 0L, 9L, 7L, 5L,
0L, 2L, 0L, 30L, 0L, 0L, 2L, 1L, 0L, 4L, 0L, 4L, 0L, 0L,
0L, 0L, 0L, 2L, 1L, 0L, 0L, 3L, 0L, 0L, 0L, 0L, 0L, 0L, 5L,
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L), Group2_LU = c(2L,
0L, 0L, 0L, 0L, 2L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L,
0L, 1L, 1L, 6L, 0L, 1L, 2L, 0L, 0L, 0L, 0L, 0L, 1L, 13L,
4L, 2L, 0L, 0L, 0L, 2L, 0L, 2L, 2L, 0L, 6L, 0L, 0L, 0L, 0L,
0L, 0L, 7L, 0L, 3L, 2L, 1L, 0L, 2L, 4L, 0L, 0L, 1L, 0L, 4L,
0L, 0L, 1L, 0L, 2L, 1L, 1L, 0L, 0L, 0L), Group1_ME = c(7L,
0L, 1L, 0L, 16L, 2L, 2L, 1L, 2L, 0L, 0L, 0L, 1L, 0L, 1L,
1L, 0L, 0L, 0L, 2L, 0L, 2L, 3L, 0L, 0L, 0L, 1L, 1L, 1L, 8L,
13L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 3L, 0L, 5L, 0L, 3L, 1L,
0L, 1L, 0L, 1L, 1L, 4L, 1L, 0L, 0L, 1L, 5L, 0L, 0L, 1L, 0L,
6L, 0L, 2L, 0L, 0L, 4L, 0L, 1L, 0L, 0L, 0L), Group2_ME = c(12L,
0L, 1L, 0L, 4L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 2L, 0L, 0L, 0L,
0L, 0L, 4L, 4L, 0L, 4L, 3L, 0L, 0L, 0L, 0L, 0L, 0L, 10L,
3L, 0L, 0L, 0L, 1L, 1L, 0L, 1L, 1L, 0L, 8L, 0L, 0L, 0L, 0L,
1L, 0L, 0L, 1L, 5L, 1L, 0L, 0L, 1L, 5L, 0L, 0L, 1L, 0L, 3L,
0L, 0L, 0L, 0L, 6L, 2L, 3L, 1L, 0L, 0L), Group1_CA = c(6L,
0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 2L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 1L, 0L, 2L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 12L,
6L, 2L, 0L, 0L, 1L, 2L, 0L, 0L, 0L, 0L, 6L, 0L, 0L, 1L, 0L,
0L, 2L, 2L, 1L, 1L, 3L, 0L, 0L, 1L, 4L, 0L, 0L, 0L, 0L, 4L,
0L, 0L, 2L, 0L, 5L, 1L, 1L, 1L, 1L, 0L), Group2_CA = c(3L,
0L, 0L, 0L, 0L, 2L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L,
1L, 0L, 0L, 1L, 0L, 3L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 10L,
7L, 0L, 0L, 0L, 0L, 2L, 0L, 1L, 2L, 0L, 4L, 0L, 0L, 0L, 0L,
2L, 0L, 0L, 0L, 3L, 3L, 0L, 0L, 1L, 2L, 0L, 0L, 0L, 0L, 2L,
0L, 0L, 1L, 0L, 3L, 2L, 1L, 0L, 0L, 0L), Group1_CR = c(6L,
0L, 0L, 0L, 13L, 1L, 0L, 0L, 3L, 0L, 0L, 1L, 1L, 0L, 0L,
0L, 0L, 0L, 1L, 1L, 0L, 4L, 1L, 0L, 0L, 0L, 1L, 0L, 5L, 8L,
6L, 0L, 10L, 0L, 15L, 1L, 0L, 2L, 4L, 0L, 3L, 0L, 0L, 8L,
0L, 0L, 0L, 2L, 1L, 1L, 2L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L,
2L, 0L, 0L, 1L, 2L, 6L, 0L, 2L, 0L, 0L, 0L), Group2_CR = c(2L,
0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
1L, 5L, 0L, 2L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 7L, 4L,
2L, 0L, 0L, 1L, 1L, 0L, 0L, 1L, 0L, 2L, 0L, 0L, 1L, 0L, 0L,
0L, 2L, 2L, 1L, 3L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L,
0L, 1L, 0L, 4L, 0L, 1L, 0L, 0L, 0L), Group1_VG = c(11L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 2L, 1L, 0L, 0L, 0L,
0L, 0L, 1L, 3L, 0L, 1L, 2L, 0L, 0L, 0L, 1L, 0L, 0L, 12L,
5L, 1L, 0L, 0L, 0L, 3L, 0L, 3L, 0L, 0L, 5L, 0L, 0L, 0L, 0L,
0L, 0L, 5L, 3L, 2L, 1L, 0L, 0L, 1L, 4L, 0L, 0L, 0L, 0L, 6L,
0L, 0L, 0L, 0L, 4L, 0L, 2L, 0L, 0L, 0L), Group2_VG = c(3L,
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 2L, 0L, 0L, 0L,
1L, 0L, 0L, 3L, 0L, 2L, 3L, 0L, 0L, 0L, 0L, 0L, 1L, 11L,
11L, 1L, 0L, 0L, 0L, 2L, 0L, 1L, 2L, 0L, 6L, 0L, 0L, 0L,
0L, 0L, 0L, 6L, 4L, 3L, 2L, 0L, 0L, 1L, 6L, 0L, 0L, 2L, 0L,
4L, 0L, 0L, 1L, 0L, 2L, 1L, 2L, 1L, 2L, 0L), Group1_MZ = c(8L,
0L, 0L, 4L, 30L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L,
1L, 0L, 0L, 0L, 2L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 3L, 17L,
6L, 1L, 23L, 4L, 26L, 1L, 0L, 1L, 3L, 0L, 6L, 0L, 0L, 1L,
0L, 0L, 0L, 3L, 2L, 2L, 6L, 0L, 2L, 1L, 3L, 0L, 0L, 0L, 1L,
4L, 0L, 0L, 0L, 3L, 4L, 2L, 1L, 0L, 0L, 0L), Group2_MZ = c(0L,
13L, 0L, 0L, 0L, 19L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L,
0L, 0L, 0L, 0L, 13L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 3L,
5L, 6L, 2L, 0L, 0L, 0L, 0L, 1L, 12L, 0L, 0L, 2L, 0L, 0L,
0L, 0L, 1L, 0L, 2L, 1L, 0L, 1L, 0L, 0L, 0L, 1L, 3L, 0L, 0L,
0L, 8L, 0L, 5L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 20L), Group1_GO = c(3L,
0L, 2L, 1L, 5L, 2L, 0L, 0L, 3L, 0L, 0L, 0L, 2L, 0L, 0L, 0L,
1L, 8L, 1L, 8L, 0L, 4L, 1L, 0L, 0L, 0L, 1L, 0L, 4L, 16L,
2L, 1L, 5L, 0L, 8L, 1L, 0L, 2L, 1L, 0L, 11L, 0L, 3L, 2L,
0L, 0L, 0L, 3L, 1L, 2L, 2L, 0L, 1L, 3L, 2L, 0L, 0L, 0L, 0L,
5L, 0L, 0L, 0L, 1L, 2L, 0L, 0L, 0L, 0L, 0L), Group2_GO = c(10L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 3L, 0L, 0L, 1L,
0L, 5L, 1L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 14L,
11L, 2L, 0L, 0L, 0L, 1L, 0L, 0L, 2L, 0L, 7L, 0L, 0L, 0L,
0L, 0L, 0L, 3L, 2L, 1L, 5L, 0L, 0L, 1L, 5L, 0L, 0L, 1L, 0L,
3L, 0L, 0L, 0L, 0L, 5L, 2L, 2L, 0L, 0L, 0L)), class = "data.frame", row.names = c("Aegithalos_caudatus",
"Alauda_arvensis", "Alcedo_atthis", "Anas_crecca", "Anas_platyrhynchos",
"Anthus_pratensis", "Anthus_spinoletta", "Ardea_alba", "Ardea_cinerea",
"Botaurus_stellaris", "Bubulcus_ibis", "Carduelis_carduelis",
"Certhia_brachydactyla", "Chloris_chloris", "Chroicocephalus_ridibundus",
"Coccothraustes_coccothraustes", "Coloeus_monedula", "Columba_livia",
"Columba_palumbus", "Corvus_cornix", "Corvus_corone", "Cyanistes_caeruleus",
"Dendrocopos_major", "Dryobates_minor", "Dryocopus_martius",
"Egretta_garzetta", "Emberiza_cia", "Emberiza_citrinella", "Emberiza_schoeniclus",
"Erithacus_rubecula", "Fringilla_coelebs", "Fringilla_montifringilla",
"Fulica_atra", "Gallinago_gallinago", "Gallinula_chloropus",
"Garrulus_glandarius", "Linaria_cannabina", "Motacilla_alba",
"Motacilla_cinerea", "Netta_rufina", "Parus_major", "Passer_montanus",
"Phalacrocorax_carbo", "Phylloscopus_collybita", "Phasianus_colchicus",
"Phoenicurus_ochruros", "Periparus_ater", "Pica_pica", "Picus_viridis",
"Poecile_palustris", "Prunella_modularis", "Psittacula_krameri",
"Rallus_aquaticus", "Regulus_ignicapilla", "Regulus_regulus",
"Saxicola_rubicola", "Serinus_serinus", "Sitta_europaea", "Spatula_clypeata",
"Spinus_spinus", "Streptopelia_decaocto", "Sturnus_vulgaris",
"Sylvia_atricapilla", "Tachybaptus_ruficollis", "Troglodytes_troglodytes",
"Turdus_iliacus", "Turdus_merula", "Turdus_philomelos", "Turdus_pilaris",
"Vanellus_vanellus"))
library(iNEXT)
abbcove <-estimateD(df,q=0,datatype="abundance",base="coverage",nboot=1000)
In the abbcove object you can find qD (diversity index), qD.LCL (lower bound of confidence interval for the diversity index) and qd.UCL (upper bound of confidence interval for the diversity index)
