# How to do a power analysis on density (count/area) data?

I would like to calculate the power to detect a hypothetical x% decrease from a mean value, but have a few qualifying questions. The data are strictly positive, right-skewed, continuous data containing many zeros (density data). A log(x+1) was insufficient in transforming the data to resemble normal data.

1. Is there a package in R that can run a power analysis on this data? For example, can pwr.t.test handle my data's departure from normally distributed data?

Seperate, but related questions...

1. If the power to detect a decline for 1 year's worth of data (n=47) is <80%, and since I don't have the money to increase my sample size each year, how many years of previously recorded data would I need to combine in order to detect an x% decrease in a mean (of ___ years of data), in the future? i.e. Would combing 2022 and 2023 data (n=94) result in a sufficient sample size to detect with 80%, or more, a hypothetical x% decline in 2025 (2 years ahead of 2023)?

Is combining datasets from previous years an "OK" method for increasing sample size? NOTE! The data are paired (same set of blocks sampled once a year, n=47). What are some things I need to keep in mind/account for when artificially increasing my sample size in this way?

"den" = density value per site, year, and season. "ln_den" = log(x+1) transformation of "den" values.

UPDATE: Distribution of full dataset (month x year, note - not always full 47):

    > table(data$$Month, data$$CYR)

2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 2022 2023
1    47   47   47   47   47    0    0    0    0    0    0    0    0    0    0    0
2     0    0    0    0    0    0    0    0    0    0   20    0    0    0    9    0
3     0    0    0    0    0   47   39   30   24   43   27   47   47    0   38   47
4     0    0    0    0    0    0    8   17   23    0    0    0    0    0    0    0
7    47   47   47   47    0    0    0    0    0    0    0    0    0    0    0    0
8     0    0    0    0   47    0   30    0    0    0    0    0    0    0    0    0
9     0    0    0    0    0   47   17   47   47    0   47   47   47   47   47   47
10    0    0    0    0    0    0    0    0    0   47    0    0    0    0    0    0


Data:

    > dput(gold)
structure(list(CYR = c(2022L, 2022L, 2022L, 2022L, 2022L, 2022L,
2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L,
2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L,
2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L,
2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L,
2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L,
2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L,
2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L,
2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L,
2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L,
2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2022L, 2023L, 2023L,
2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L,
2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L,
2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L,
2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L,
2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L,
2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L,
2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L,
2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L,
2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L,
2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L, 2023L,
2023L, 2023L), Season = c("DRY", "DRY", "DRY", "DRY", "DRY",
"DRY", "DRY", "DRY", "DRY", "DRY", "DRY", "DRY", "DRY", "DRY",
"DRY", "DRY", "DRY", "DRY", "DRY", "DRY", "DRY", "DRY", "DRY",
"DRY", "DRY", "DRY", "DRY", "DRY", "DRY", "DRY", "DRY", "DRY",
"DRY", "DRY", "DRY", "DRY", "DRY", "DRY", "DRY", "DRY", "DRY",
"DRY", "DRY", "DRY", "DRY", "DRY", "DRY", "WET", "WET", "WET",
"WET", "WET", "WET", "WET", "WET", "WET", "WET", "WET", "WET",
"WET", "WET", "WET", "WET", "WET", "WET", "WET", "WET", "WET",
"WET", "WET", "WET", "WET", "WET", "WET", "WET", "WET", "WET",
"WET", "WET", "WET", "WET", "WET", "WET", "WET", "WET", "WET",
"WET", "WET", "WET", "WET", "WET", "WET", "WET", "WET", "DRY",
"DRY", "DRY", "DRY", "DRY", "DRY", "DRY", "DRY", "DRY", "DRY",
"DRY", "DRY", "DRY", "DRY", "DRY", "DRY", "DRY", "DRY", "DRY",
"DRY", "DRY", "DRY", "DRY", "DRY", "DRY", "DRY", "DRY", "DRY",
"DRY", "DRY", "DRY", "DRY", "DRY", "DRY", "DRY", "DRY", "DRY",
"DRY", "DRY", "DRY", "DRY", "DRY", "DRY", "DRY", "DRY", "DRY",
"DRY", "WET", "WET", "WET", "WET", "WET", "WET", "WET", "WET",
"WET", "WET", "WET", "WET", "WET", "WET", "WET", "WET", "WET",
"WET", "WET", "WET", "WET", "WET", "WET", "WET", "WET", "WET",
"WET", "WET", "WET", "WET", "WET", "WET", "WET", "WET", "WET",
"WET", "WET", "WET", "WET", "WET", "WET", "WET", "WET", "WET",
"WET", "WET", "WET"), Month = c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L,
9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L,
9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L,
9L, 9L, 9L, 9L, 9L, 9L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L,
9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L,
9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L,
9L, 9L, 9L, 9L), Day = c(28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L,
28L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 18L, 18L, 18L, 18L, 18L,
18L, 18L, 18L, 18L, 18L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 12L, 12L, 12L,
12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 13L, 13L, 13L, 13L, 13L,
13L, 13L, 13L, 13L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L,
8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L,
10L, 10L, 10L, 10L, 10L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L,
6L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 8L, 8L,
8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 9L, 9L,
9L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 14L, 14L,
14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 15L, 15L, 15L,
15L, 15L, 15L, 15L, 15L, 15L, 18L, 18L, 18L, 18L, 18L, 18L, 18L,
18L, 18L, 18L, 19L, 19L, 19L, 19L, 19L, 19L), Site = c(39, 40,
41, 42, 43, 44, 45, 46, 47, 1, 2, 3, 4, 5, 6, 7, 10, 11, 12,
13, 14, 15, 16, 17, 8, 9, 18, 19, 20, 21, 22, 23, 24, 25, 26,
27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 1, 10, 11, 2,
3, 4, 5, 6, 7, 8, 9, 12, 13, 14, 15, 16, 17, 18, 19, 20, 40,
41, 42, 43, 44, 45, 46, 47, 29, 30, 31, 32, 33, 34, 35, 36, 37,
21, 22, 23, 24, 25, 26, 27, 28, 38, 39, 10, 11, 12, 13, 9, 38,
39, 40, 41, 42, 43, 44, 45, 46, 47, 18, 19, 20, 21, 22, 23, 24,
25, 26, 27, 28, 29, 30, 14, 15, 16, 17, 31, 32, 33, 34, 35, 36,
37, 1, 2, 3, 4, 5, 6, 7, 8, 38, 39, 40, 41, 42, 43, 44, 45, 46,
47, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 15, 30, 31,
32, 33, 34, 35, 36, 37, 1, 10, 11, 2, 3, 5, 6, 7, 8, 9, 12, 13,
14, 16, 17, 4), common_name = c("Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish",
"Goldspotted Killifish", "Goldspotted Killifish", "Goldspotted Killifish"
), num = c(0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
8L, 4L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 3L, 0L, 5L, 5L, 1L, 2L, 1L, 0L,
3L, 4L, 1L, 0L, 0L, 0L, 6L, 1L, 3L, 1L, 0L, 0L, 1L, 0L, 1L, 1L,
2L, 18L, 0L, 16L, 9L, 2L, 3L, 3L, 1L, 1L, 3L, 8L, 25L, 1L, 0L,
6L, 5L, 6L, 2L, 1L, 0L, 14L, 3L, 0L, 4L, 4L, 2L, 2L, 4L, 6L,
10L, 0L, 4L, 0L, 1L, 0L, 0L, 2L, 0L, 2L, 0L, 0L, 0L, 1L, 3L,
0L, 0L, 0L, 1L, 1L, 0L, 1L, 0L, 0L, 3L, 1L, 0L, 1L, 1L, 0L, 0L,
0L, 0L, 1L, 0L, 4L, 0L, 0L, 3L, 2L, 16L, 7L, 0L, 0L, 4L, 0L,
0L, 1L, 2L, 0L, 1L, 0L, 13L, 5L, 0L, 24L, 1L, 3L, 20L, 0L, 12L,
10L, 15L, 1L, 2L, 1L, 13L, 7L, 2L, 7L, 2L, 3L, 1L, 9L, 2L, 10L,
12L, 3L, 0L, 18L, 7L, 2L, 0L, 12L, 3L, 8L, 9L, 4L, 9L, 10L, 15L,
0L, 1L, 1L, 0L, 5L, 1L)), row.names = c(387311L, 387497L, 387683L,
387869L, 388055L, 388241L, 388427L, 388613L, 388799L, 388985L,
389171L, 389357L, 389543L, 389729L, 390101L, 390287L, 390473L,
390659L, 391031L, 391217L, 391589L, 391775L, 391961L, 392147L,
392519L, 392705L, 392891L, 393077L, 393263L, 393449L, 393635L,
393821L, 394007L, 394193L, 394379L, 394565L, 394751L, 394937L,
395123L, 395309L, 395495L, 395681L, 395867L, 396053L, 396239L,
396425L, 396797L, 396983L, 397169L, 397355L, 397541L, 397727L,
397913L, 398099L, 398471L, 398657L, 399029L, 399215L, 399587L,
399773L, 400145L, 400331L, 400517L, 400703L, 400889L, 401075L,
401261L, 401447L, 401633L, 401819L, 402005L, 402191L, 402377L,
402563L, 402749L, 402935L, 403121L, 403307L, 403493L, 403679L,
403865L, 404051L, 404237L, 404423L, 404609L, 404795L, 404981L,
405167L, 405353L, 405539L, 405725L, 405911L, 406283L, 406469L,
406655L, 406841L, 407213L, 407399L, 407585L, 407957L, 408143L,
408329L, 408515L, 408701L, 408887L, 409073L, 409259L, 409445L,
409631L, 409817L, 410003L, 410189L, 410375L, 410561L, 410747L,
410933L, 411119L, 411305L, 411491L, 411677L, 411863L, 412049L,
412421L, 412607L, 412793L, 412979L, 413165L, 413351L, 413537L,
413723L, 413909L, 414095L, 414281L, 414467L, 414653L, 414839L,
415025L, 415211L, 415583L, 415769L, 416141L, 416513L, 416699L,
416885L, 417071L, 417257L, 417443L, 417629L, 417815L, 418001L,
418187L, 418373L, 418559L, 418745L, 418931L, 419117L, 419303L,
419489L, 419675L, 419861L, 420047L, 420233L, 420419L, 420605L,
420791L, 420977L, 421163L, 421349L, 421535L, 421721L, 421907L,
422093L, 422279L, 422465L, 422651L, 422837L, 423023L, 423209L,
423581L, 423767L, 424139L, 424325L, 424697L, 424883L, 425255L,
425441L, 425627L, 425813L), class = "data.frame")

• Could you elaborate on what you mean by "density data"? These look like fractions and seem likely not to be the original data. To answer your question well, we need to understand better how the data arise.
– whuber
Commented Jul 8 at 18:34
• Certainly, density in this case is defined as fish counts per site (a site = 3m^2).
– Nate
Commented Jul 8 at 18:39
• Counts are natural numbers --these are not. What did you divide the counts by?
– whuber
Commented Jul 8 at 18:41
• I divided the counts by the area of water we searched for fish. So, a count of 1 fish in an area of 3m^2 results in a density of 0.333, for instance.
– Nate
Commented Jul 8 at 18:43
• If you have access to prior years' data, why not just use data going back as far as you can? That would seem to provide the most power of all. Also, it's best to model count data via the original count values with a Poisson or similar model; if sites have different areas, then you can include an offset term to get results in counts per unit area. Simulations based on your prior years' data might be more reliable than power calculations that depend on assumptions that your data might not meet.
– EdM
Commented Jul 8 at 18:44

Overview

Simulation is a good way to estimate power in situations that aren't covered by simple formulas. In your case, use your current data to estimate the characteristics of the data generating process, then choose hypothetical future values that follow the same general process except, say, for a change in mean value. Take a random sample from the hypothetical values, append to the current data, build the model on the combined data, and determine whether the difference between the future value and some appropriate combination of prior values is significant.

Do that sampling from the hypothetical values and model building a thousand times or more. The fraction of samples for which you found a "significant" difference is your estimate of the power.

This answer illustrates the process for a Poisson count model based on synthetic data similar to those in the question. It's possible that your data will require a negative binomial model instead, or even something more complicated. You would apply the same principles in that case for modeling of current data and simulation of future data.

The simulation shows the danger in working with individual years if there's noise in the data generating process, emphasizing the importance of intelligent application of your understanding of the subject matter to what the statistical analysis shows.

Current data

As the question doesn't contain complete data back to 2008, simulate some. The data provided for 2022 and 2023 are consistent with a random block effect with SD of about 1 in the log(mean) scale used by a Poisson model. Randomly sample 47 such blocks. As the blocks are the same from year to year, keep those values for all data.

Assume a slight downward trend in the log of the mean counts per block. Allow for some variability around the trend.

set.seed(20240714)
randBlocks <- rnorm(47, sd=1)
meanByYear <- log(6) - 0.05*c(0:15)  + rnorm(n=16,sd=0.4)


Set up a data structure to hold all the data, with the above meanByYear and randBlocks values. Despite their names, they are both in log(mean) scales, as a Poisson GLM models the log of the mean counts. Sum them to get the log(mean) for each block/CYR combination. Sample randomly from corresponding Poisson distributions; rpois() expects mean counts, so exponentiate log(mean) values.

currentData <- data.frame(CYR=rep(2008:2023,each=47),
block=rep(1:47,16),
yearMean=rep(meanByYear,each=47),
randBlock=rep(randBlocks,16))
## add random block effects to mean for year
currentData[,"meanForYrBlock"] <- currentData$yearMean + currentData$randBlock
## sample from corresponding Poisson distributions
currentData[,"count"] <- rpois(n=752, exp(currentData[,"meanForYrBlock"]))


Fit model to current data

Now fit a Poisson model. This uses a spline fit for CYR. One might consider treating CYR instead as a categorical predictor, depending on your understanding of the subject matter.

## extract data for modeling
currentCountData <- currentData[,c("CYR","block","count")]
##
## fit spline model
library(lme4)
library(splines)
currentSplineModel <- glmer(count~ns(CYR, df=3) + (1|block),
data=currentCountData, family="poisson")


Plot the current data and model estimates

plot(2008:2023, meanByYear, xlab="Year",
ylab="Log mean counts per block", bty="n",
xlim=c(2005,2025))
lines(2008:2023,log(6) - 0.05*c(0:15)) ## trend intended
currentDataPreds <- predict(currentSplineModel,
newdata=data.frame(CYR=2008:2023),
re.form=NA) ## main effects only
lines(2008:2023, currentDataPreds, col="red") ## model estimates
legend("topright", "black line, intended trend\ndots, trend plus error\nred line, mixed model estimates",bty="n")


Examine the plot. The spline does a reasonable job at capturing the underlying trend. Note that some individual years (2015, 2016) differ greatly from the underlying trend, however. The danger in putting too much emphasis on individual years is that you might just be capturing noise.

Post-modeling analysis, current data

Use emmeans to compare model scenarios. The at() parameter sets up a grid with predictions for each year, named CYR2008, etc. If CYR had been treated as a factor, the at() parameter wouldn't be needed as all factor levels would be included by default. Some functions require using indices instead of names.

The ref() argument to contrast() specifies a set of reference levels against which to compare other levels. The exclude() argument specifies levels to omit from the analysis. This example compares the 8th level (CYR2015) against all previous. These results are in the log(mean) scale; you can ask for the "response" (counts) scale instead.

library(emmeans)
emmCurrent <- emmeans(currentSplineModel, specs=c("CYR"),
at=list(CYR=2008:2023)) ## get each year
contrast_2015vsPrior <- data.frame(contrast(emmCurrent,
method="trt.vs.ctrl", ref=1:7,
exclude=9:16))
contrast_2015vsPrior[,1];contrast_2015vsPrior[,2:6]
## [1] "CYR2015 - avg(CYR2008,CYR2009,CYR2010,CYR2011,CYR2012,CYR2013,CYR2014)"

##     estimate         SE  df  z.ratio      p.value
## 1 -0.2318068 0.02210378 Inf -10.4872 9.891488e-26


This shows the danger in putting too much emphasis on a single year. The log(mean) value per block is highest of all in 2015 but the model, allowing for a trend (and taking advantage of later information about the data generating process), estimates that the underlying 2015 value was significantly lower than the average over previous years.

Generate hypothetical future data

Assume that the "random effects" are maintained the same for each block into the future. Write a function that allows for a hypothetical log(mean) value for CYR2024, adds in the block "random effects," and then samples from the corresponding Poisson distributions. Illustrated with log(mean) of 1.25.

## define simulation function
return(counts)
}
## hypothesize log(mean) of 1.25 for 2024, put into data frame
data2024 <- data.frame(CYR=rep(2024,47), block=1:47,


For power estimates you would do the above several thousand times to get a large set of hypothetical (future, e.g., 2024 or beyond) data from which to draw. For each hypothetical set of data, do something like the following:

fullData <- rbind(currentCountData, data2024)
model2024 <- glmer(count~ns(CYR,df=3) + (1|block),
data=fullData,family="poisson")
emm2024 <- emmeans(model2024,
specs=c("CYR"),
at=list(CYR=2008:2024))


You then run the comparison that you want for each of the several thousand hypothetical 2024 data sets. You can extract the estimated difference and the p-value, as shown below. The fraction of data sets for which you find a "significant" difference (at a specified false-positive rate) between 2024 and your reference years is the estimate of the power.

The problem is that such comparisons depend heavily upon your choice of reference years. With the overall trend I assumed, if you compare 2024 against all previous years it's quite significantly lower than the reference years.

contrast_2024vsAllPrev <- data.frame(contrast(emm2024,
method="trt.vs.ctrl", ref=c(1:16)))
contrast_2024vsAllPrev[,1];contrast_2024vsAllPrev[,2:6]
## [1] "CYR2024 - avg(CYR2008,CYR2009,CYR2010,CYR2011,CYR2012,CYR2013,CYR2014,CYR2015,CYR2016,CYR2017,CYR2018,CYR2019,CYR2020,CYR2021,CYR2022,CYR2023)"
##     estimate         SE  df   z.ratio      p.value
## 1 -0.4390201 0.0528057 Inf -8.313877 9.26253e-17
contrast_2024vsAllPrev\$p.value
## [1] 9.26253e-17


If you instead restrict reference values to the previous 3 years, it isn't.

contrast_2024vsPrev3 <- data.frame(contrast(emm2024,
method="trt.vs.ctrl",
exclude=1:13, ref=c(14:16)))
contrast_2024vsPrev3[,1];contrast_2024vsPrev3[,2:6]
## [1] "CYR2024 - avg(CYR2021,CYR2022,CYR2023)"
##      estimate         SE  df  z.ratio     p.value
## 1 -0.02897765 0.03422821 Inf -0.8466014 0.3972173


The choice of reference years should primarily be based on your understanding of the subject matter, informed by the statistical analysis of your current data.

The first step to start to do what you want is to fit an appropriate model to all the count data that you have for that species, over all years. Working with "densities" ignores the variability that comes from low numbers of counts. You would include calendar year as a predictor in some way, to check for prior trends. As you want to be able to estimate a decrease from a mean value over prior years, you need to document that such a mean value isn't already hiding a lot of variability.

The snippet of data that you show suggest that there's already a lot of year-to-year difference. As your "density" values are evidently over 3 square meters, I multiplied all the values by 3 and compared these two calendar years. There's a massive and highly significant difference: 180 counts in 2022, and 284 in 2023. That's over a 50% increase from 2022 to 2023. A simple Poisson count model with CYR as the predictor puts the p-value at 1.70e-06.

I'm thus skeptical that "detect[ing] a hypothetical x% decrease from a mean value" would be meaningful in this case, at least in the way that you propose. With large year-to-year variations, the "mean value" taken over some set of previous years will depend heavily on your choice of the previous years. That makes the "decrease from [the] mean value" essentially uninterpretable. You would need to document that the values in the coming year(s) went beyond the typical year-to-year variation in counts.

If you have enough years of prior data you might be able to build a Poisson or negative binomial model for the counts that could inform your work. To account for the repeated measures on the same "blocks" over time, use a mixed model with the "blocks" treated as random effects. The R lme4 package provides tools for both Poisson and negative binomial models.

You should work with an experienced local statistician to decide how best to model the CYR values. At first thought I can see arguments for modeling them continuously with a regression spline, treating them as fixed factors, or including them as random factors. Much would depend on the number of years for which you have data and on your understanding of the subject matter.

Once the fundamental characteristics of your data have been modeled, you (or your statistician) can simulate large numbers of data values to represent what you think might occur in a subsequent year or years. For example, if a simple Poisson model is adequate, you could assume a new (lower) mean count value per block for a coming year and the observed inter-block variability from your mixed model to simulate those future data, by drawing a very large number of random values with a function like rpois() in R. Then take multiple samples from your simulated data, run the proper comparison against the historical data (probably not just against the mean value) at a specified false-positive rate (typically 0.05), and see how frequently among those samples that you can detect your hypothesized lower mean count as being "significant." That fraction of samples in which you find a "significant" difference is the power estimate.

First, plot the data that you have to see what the annual changes are. That might give you a gut-level sense of how much a value in some future year needs to be below what there already is in order to say that there's a believable difference. I suspect that it might have to be a very, very low value.

Simulating data isn't so hard. For Poisson data, all you need is your estimate of the mean number of counts (lambda) for each hypothesized combination of year, season, and block. Then in R rpois(N,lambda) will return N values sampled from a Poisson distribution having mean value lambda. Use a very large value of N (many thousands) so that you can sample repeatedly when you do the power estimate.*

You pretty much get to choose the values of lambda, to correspond to your hypothesized low number of counts in future years. Although you are modeling block as a random effect, with the same blocks tested each year you could add, effectively as fixed effects, the individual random effects found in your model of data to date to the fixed-effect estimates. Be careful, as the coefficients of the linear model (including the individual random effects) will be in the log scale, and you need to do that addition in the log scale before you convert to a value of lambda.

The correct test, after you model the data to date together with a single sample of your hypothesized future data, would be a post-hoc test comparing the model estimate for future year (or years) to some appropriate combination of the estimates for the prior years. The R emmeans package has tools that can perform tests on arbitrary sets or linear combinations of model predictions. You could test the future year against each of the prior years, or against the average of some pre-specified set of prior years.

Then do that for a few thousand samples from your simulated future data plus the known data to date. The fraction of samples that give you a "significant" difference is the power estimate you seek.

The "appropriate combination of prior years" will depend on your understanding of the subject matter and, given the big difference between 2022 and 2023, is probably going to be difficult to justify if there are already trends or large variations. Testing against each of the previous years might be of interest, but that will involve large multiple-comparison corrections that will certainly lower the power.

*If you need a negative binomial model for your data to date then there are two parameters for rnbinom(). You'd probably want to use the "alternative parameterization" noted in its help page, keeping the dispersion parameter that was found in your model and simulating for your choice of mean value.

• thank you (+1)! I especially appreciate your last paragraph. I do know how to set-up both GLM(M)s and GAM(M)s and your suggestions sound like a great approach. Unfortunately (for me) I don't know any statisticians (and my advisor has a "linear regression or nothing" outlook on stats), I also don't know how to simulate new data from said model or run that correct comparison you mentioned lastly.
– Nate
Commented Jul 12 at 17:25
• I just added a table of my full dataset, and although there are a lot of gaps, I would argue for a spline of year and a 2 level factor for season (instead of month - we only have 2 seasons) in a model.
– Nate
Commented Jul 12 at 17:52
• would a quick example be possible in R?
– Nate
Commented Jul 12 at 22:00
• That would take a few days; I do have other things to do. While you're waiting, could you please edit the question to show a plot of the counts per year summed over all blocks? I recognize that there are a couple of missing values in 2017 and 2021, but seeing the full count data would make the situation a lot more concrete.
– EdM
Commented Jul 13 at 2:15
• Absolutely, thank you for your time.
– Nate
Commented Jul 13 at 12:34