Libraries used:

library(tidyverse)
library(lme4) # for mixed models
library(emmeans) # for marginal means

Phenology

Data:

phenology2019 <- read.csv("phenology_2019.csv")
phenology2020 <- read.csv("phenology_2020.csv")
phenology2019 |> 
  head(10) |> 
  knitr::kable()

Date	Day	Week	Location	Site	Treatment	Plot	PlotU	Plot4	Species	Anthesis	Floral.density	Group
5/29/2019	149	1	Toolik	Dry	C	5	69	C	AA	1	0	DS
5/29/2019	149	1	Toolik	Dry	C	6	70	C	AA	1	0	DS
5/29/2019	149	1	Toolik	Dry	C	8	72	C	AA	1	0	DS
5/29/2019	149	1	Toolik	Dry	C	10	74	C	AA	1	0	DS
5/29/2019	149	1	Toolik	Dry	C	11	75	C	AA	1	7	DS
5/29/2019	149	1	Toolik	Dry	C	13	77	C	AA	1	2	DS
5/29/2019	149	1	Toolik	Dry	C	16	80	C	AA	1	3	DS
5/29/2019	149	1	Toolik	Dry	E	2	82	y27	AA	1	0	DS
5/29/2019	149	1	Toolik	Dry	E	7	87	y27	AA	1	0	DS
5/29/2019	149	1	Toolik	Dry	E	8	88	y27	AA	1	0	DS

Auxiliary functions:

f_mean <- function(x){
  if (sum(!is.na(x)) == 0) {
    return(NA)
  } else {
    return(mean(x, na.rm = T))
  }
}
f_sum <- function(x){
  if (sum(!is.na(x)) == 0) {
    return(NA)
  } else {
    return(sum(x, na.rm = T))
  }
}

# place stars on p-values:
f_star <- function(x) {
  if (is.na(x)) {
    return(NA)
  } else if (x < 0.05 & x >= 0.01) {
    return("*")
  } else if (x < 0.01 & x >= 0.001) {
    return("**")
  } else if (x < 0.001) {
    return("***")
  } else {
    return(NA)
  }
}

Quasi-Poisson adjustment to glmer:

f_qp <- function(mod, coefs = coef(summary(mod))) {
  phi <- sum(residuals(mod, type = "pearson")^2)/df.residual(mod)
  tab <- within(as.data.frame(coefs), {
    `Std. Error` <- `Std. Error`*sqrt(phi)
    `z value` <- Estimate/`Std. Error`
    `Pr(>|z|)` <- 2*pnorm(abs(`z value`), lower.tail = FALSE)
  })
  return(as_tibble(tab))
}

Number of species that bloomed in each plot and mean species floral density in each plot:

# All species together, 2019
plot_summaries_2019 <- phenology2019 |> 
  group_by(PlotU, Plot4, Location, Site) |> 
summarize(
    n_sp = length(levels(factor(Species))),
    fd_mean = f_mean(Floral.density), 
    fd_sum = f_sum(Floral.density),
    n = n()
    ) |> 
  filter(!is.na(PlotU)) |> 
  ungroup() |> 
  mutate(PlotU = factor(PlotU))

# By growth forms, 2019
plot_summaries_GF_2019 <- phenology2019 |> 
  group_by(PlotU, Plot4, Location, Site, Group) |> 
summarize(
    n_sp = length(levels(factor(Species))),
    fd_mean = f_mean(Floral.density), 
    fd_sum = f_sum(Floral.density), 
    n = n()
    ) |> 
  filter(!is.na(PlotU)) |> 
  ungroup() |> 
  mutate(PlotU = factor(PlotU))

# All species together, 2020
plot_summaries_2020 <- phenology2020 |> 
  group_by(PlotU, Plot4, Location, Site) |> 
summarize(
    n_sp = length(levels(factor(Species)))
    ) |> 
  filter(!is.na(PlotU)) |> 
  ungroup() |> 
  mutate(PlotU = factor(PlotU))

# By growth forms, 2020
plot_summaries_GF_2020 <- phenology2020 |> 
  group_by(PlotU, Plot4, Location, Site, Group) |> 
summarize(
    n_sp = length(levels(factor(Species)))
    ) |> 
  filter(!is.na(PlotU)) |> 
  ungroup() |> 
  mutate(PlotU = factor(PlotU))

Number of species in bloom per plot per week and mean weekly floral density in each plot per week:

# All species together, 2019
weekly_summaries_2019 <- phenology2019 |> 
  group_by(Week, PlotU, Plot4, Location, Site) |> 
summarize(
    n_sp = length(levels(factor(Species))),
    fd_mean = f_mean(Floral.density),
    fd_sum = f_sum(Floral.density),
    n = n()
    ) |> 
  filter(!is.na(PlotU)) |> 
  ungroup() |> 
  mutate(PlotU = factor(PlotU))

# By growth forms, 2019
weekly_summaries_GF_2019 <- phenology2019 |> 
  group_by(Week, PlotU, Plot4, Location, Site, Group) |> 
summarize(
    n_sp = length(levels(factor(Species))),
    fd_mean = f_mean(Floral.density),
    fd_sum = f_sum(Floral.density),
    n = n()
    ) |> 
  filter(!is.na(PlotU)) |> 
  ungroup() |> 
  mutate(PlotU = factor(PlotU))

# All species together, 2020
weekly_summaries_2020 <- phenology2020 |> 
  group_by(Week, PlotU, Plot4, Location, Site) |> 
summarize(
    n_sp = length(levels(factor(Species)))
    ) |> 
  filter(!is.na(PlotU)) |> 
  ungroup() |> 
  mutate(PlotU = factor(PlotU))

# By growth forms, 2020
weekly_summaries_GF_2020 <- phenology2020 |> 
  group_by(Week, PlotU, Plot4, Location, Site, Group) |> 
summarize(
    n_sp = length(levels(factor(Species)))
    ) |> 
  filter(!is.na(PlotU)) |> 
  ungroup() |> 
  mutate(PlotU = factor(PlotU))

Floral density (2019 only)

Checking for over/under dispersion:

m <- glmer(fd_sum ~ Plot4 + offset(log(n)) + (1 | Location/Site), family = poisson, data = plot_summaries_2019)
deviance(m)/df.residual(m)

## [1] 17.68742

There’s overdispersion. Will use a negative binomial family.

Model selection:

m1 <- glmer.nb(fd_sum ~ Plot4 + offset(log(n)) + (1 | Location/Site), data = plot_summaries_2019)
m2 <- glmer.nb(fd_sum ~ Plot4 + offset(log(n)) + (1 | Site), data = plot_summaries_2019)
m3 <- glmer.nb(fd_sum ~ Plot4 + offset(log(n)) + (1 | Location), data = plot_summaries_2019)
anova(m1, m2)

## Data: plot_summaries_2019
## Models:
## m2: fd_sum ~ Plot4 + offset(log(n)) + (1 | Site)
## m1: fd_sum ~ Plot4 + offset(log(n)) + (1 | Location/Site)
##    npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## m2    6 983.08 1000.2 -485.54   971.08                    
## m1    7 989.02 1009.0 -487.51   975.02     0  1          1

anova(m1, m3)

## Data: plot_summaries_2019
## Models:
## m3: fd_sum ~ Plot4 + offset(log(n)) + (1 | Location)
## m1: fd_sum ~ Plot4 + offset(log(n)) + (1 | Location/Site)
##    npar     AIC  BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## m3    6 1000.88 1018 -494.44   988.88                         
## m1    7  989.02 1009 -487.51   975.02 13.855  1  0.0001975 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(m2, m3)

## Data: plot_summaries_2019
## Models:
## m2: fd_sum ~ Plot4 + offset(log(n)) + (1 | Site)
## m3: fd_sum ~ Plot4 + offset(log(n)) + (1 | Location)
##    npar     AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## m2    6  983.08 1000.2 -485.54   971.08                    
## m3    6 1000.88 1018.0 -494.44   988.88     0  0

A model with only community type (Site) as a random effect gives the best AIC. Nesting Site within Location does not improve the model. Will proceed with only community type as a random effect.

All species:

# compute EMMs, SEs, and p-values
m <- glmer.nb(fd_sum ~ Plot4 + offset(log(n)) + (1 | Location/Site), data = plot_summaries_2019)
emms <- as_tibble(emmeans(m, ~ Plot4, type = "response", offset = log(1)))
s <- summary(m)
all <- bind_cols(emms[,1:3], as_tibble(coef(s)[,3:4]))
knitr::kable(all, digits = 4)

Plot4	response	SE	z value	Pr(>\|z\|)
C	0.5056	0.1335	-2.5831	0.0098
y02	0.7325	0.2153	1.6601	0.0969
y05	0.4118	0.1482	-0.6711	0.5021
y27	1.7266	0.6254	3.9522	0.0001

By growth forms:

# auxiliary function to compute EMMs, SEs, and p-values
f <- function(gf){
  m <- glmer.nb(fd_sum ~ Plot4 + offset(log(n)) + (1 | Site), data = filter(plot_summaries_GF_2019, Group == gf))
  emms <- as_tibble(emmeans(m, ~ Plot4, type = "response", offset = log(1)))
  s <- summary(m)
  return(bind_cols(emms[,1:3], as_tibble(coef(s)[,3:4], 4)))
}

# p-values
comps <- map_df(c("F", "ES", "DS", "GS"), f) |> 
  mutate(Group = rep(c("F", "ES", "DS", "G"), each = 4), .before = Plot4)
knitr::kable(comps, digits = 4)

Group	Plot4	response	SE	z value	Pr(>\|z\|)
F	C	0.3658	0.1170	-3.1444	0.0017
F	y02	0.1000	0.0527	-2.0540	0.0400
F	y05	0.1293	0.0768	-1.5273	0.1267
F	y27	0.0322	0.0401	-1.9291	0.0537
ES	C	0.5582	0.3194	-1.0189	0.3083
ES	y02	0.6393	0.3806	0.4419	0.6585
ES	y05	0.2882	0.1922	-1.5415	0.1232
ES	y27	0.7220	0.4701	0.6325	0.5270
DS	C	1.0315	0.5290	0.0604	0.9518
DS	y02	1.7333	0.9444	1.5384	0.1239
DS	y05	0.5606	0.3691	-1.2457	0.2129
DS	y27	3.3053	2.1205	2.3802	0.0173
G	C	0.2259	0.0901	-3.7308	0.0002
G	y02	0.1428	0.0858	-0.6364	0.5245
G	y05	0.9443	0.6848	1.7286	0.0839
G	y27	0.3474	0.4035	0.3504	0.7260

Combining the floral density data for plotting:

fd <- all |> 
  mutate(Group = rep("All", 4), .before = Plot4) |> 
  bind_rows(comps) |> 
  mutate(Lower = response - SE, 
         Upper = response + SE,
         Group = factor(Group, levels = c("All", "F", "ES", "DS", "G")))
fd$`z value`[fd$Plot4 == "C"] <- NA
fd$`Pr(>|z|)`[fd$Plot4 == "C"] <- 1
fd$sig <- map_vec(fd$`Pr(>|z|)`, f_star)
levels(fd$Plot4) <- list(C  = "C", `W (1 yr)` = "y02", `W (4 yrs)` = "y05", `W (26 yrs)` = "y27")
names(fd)[1:2] <- c("Growth form", "Plot treatment")

Plot:

p <- fd |> 
  ggplot(aes(x = `Growth form`, y = response, color = `Plot treatment`)) + 
  ylab(as.expression(bquote("Floral density (flowers per 30.5 x 30.5" ~ cm^2 ~ ")"))) + 
  geom_pointrange(aes(ymin = Lower, ymax = Upper), position = position_dodge(width = 0.3)) +
  geom_text(aes(label = sig, x = `Growth form`, y = Upper + 0.17, group = `Plot treatment`), angle = 90, position = position_dodge(width = 0.34), vjust = 0.5, show.legend = FALSE) +
  theme_bw()
p

Weekly, all species:

# auxiliary function to compute EMMs, SEs, and p-values
f <- function(w){
  m <- glmer.nb(fd_sum ~ Plot4 + offset(log(n)) + (1 | Site), data = filter(weekly_summaries_2019, Week == w))
  emms <- as_tibble(emmeans(m, ~ Plot4, type = "response", offset = log(1)))
  s <- summary(m)
  return(bind_cols(emms[, 1:3], as_tibble(coef(s)[,3:4], 4)))
}

# p-values
all <- map_df(1:7, f) |> 
  mutate(Week = rep(1:7, each = 4), .before = Plot4) |> 
  mutate(Group = "All", .before = Plot4)
knitr::kable(all, digits = 4)

Week	Group	Plot4	response	SE	z value	Pr(>\|z\|)
1	All	C	0.8936	0.3923	-0.2563	0.7977
1	All	y02	1.5779	0.7596	1.6023	0.1091
1	All	y05	1.4788	0.8344	1.0683	0.2854
1	All	y27	1.8889	1.0944	1.4972	0.1343
2	All	C	0.6082	0.2224	-1.3596	0.1740
2	All	y02	1.0082	0.4111	1.4590	0.1446
2	All	y05	0.2068	0.1258	-1.9019	0.0572
2	All	y27	3.0250	1.4475	3.7271	0.0002
3	All	C	0.4848	0.1552	-2.2613	0.0237
3	All	y02	0.6229	0.2427	0.6134	0.5396
3	All	y05	0.1771	0.0968	-1.7789	0.0753
3	All	y27	1.8786	0.9287	2.6060	0.0092
4	All	C	0.4226	0.1872	-1.9440	0.0519
4	All	y02	0.6549	0.3267	1.1125	0.2659
4	All	y05	0.2475	0.1531	-0.9956	0.3194
4	All	y27	1.3517	0.7817	2.3673	0.0179
5	All	C	0.4595	0.1782	-2.0051	0.0449
5	All	y02	0.6149	0.2796	0.7334	0.4633
5	All	y05	0.4084	0.2271	-0.2305	0.8177
5	All	y27	1.2233	0.6847	1.8519	0.0640
6	All	C	0.2635	0.0713	-4.9259	0.0000
6	All	y02	0.2406	0.0917	-0.1940	0.8461
6	All	y05	0.2823	0.1460	0.1186	0.9056
6	All	y27	1.5312	0.7001	3.3122	0.0009
7	All	C	0.0703	0.0402	-4.6403	0.0000
7	All	y02	0.1048	0.0798	0.4194	0.6749
7	All	y05	0.0000	0.0000	0.0000	1.0000
7	All	y27	0.7500	0.7141	2.1311	0.0331

Weekly, by growth forms:

# Means and standard errors
df1 <- weekly_summaries_GF_2019 |> 
  filter(!is.na(Group)) |> 
  group_by(Week, Group, Plot4) |> 
  summarize(
    response = round(mean(fd_mean, na.rm = T), 2),
    SE = round(sd(fd_mean, na.rm = T) / sqrt(sum(!is.na(fd_mean))), 2)
    )

# auxiliary function to compute EMMs, SEs, and p-values
f <- function(w, gf){
  m <- glmer.nb(fd_sum ~ Plot4 + offset(log(n)) + (1 | Site), data = filter(weekly_summaries_GF_2019, Week == w, Group == gf))
  emms <- as_tibble(emmeans(m, ~ Plot4, type = "response", offset = log(1)))
  s <- summary(m)
  return(bind_cols(emms[, 1:3], as_tibble(coef(s)[,3:4], 4)))
}

# p-values
comps <- map_df(1:6, \(x)f(x, "F")) # weeks 7 and 8 not enough obs
wk <- rep(1:6, each = 4); wk <- wk[c(-2, -24)];
comps_F <- comps |> 
  mutate(Group = "F", .before = Plot4) |> 
  mutate(Week = wk, .before = Group)
df_F <- df1 |> 
  filter(Group == "F", Week == 7) |> 
  mutate(`z value` = NA, `Pr(>|z|)` = NA)
comps_F <- comps_F |> 
  bind_rows(df_F)

comps <- map_df(1:7, \(x)f(x, "ES")) # week 8 not enough obs
wk <- rep(1:7, each = 4); wk <- wk[-28];
comps_ES <- comps |> 
  mutate(Group = "ES", .before = Plot4) |> 
  mutate(Week = wk, .before = Group)

comps <- map_df(1:6, \(x)f(x, "DS")) # weeks 7 and 8 not enough obs
wk <- rep(1:6, each = 4)
comps_DS <- comps |> 
  mutate(Group = "DS", .before = Plot4) |> 
  mutate(Week = wk, .before = Group)
df_DS <- df1 |> 
  filter(Group == "DS", Week == 7) |> 
  mutate(`z value` = NA, `Pr(>|z|)` = NA)
comps_DS <- comps_DS |> 
  bind_rows(df_DS)

comps <- map_df(1:3, \(x)f(x, "GS")) # weeks 4-8 not enough obs
wk <- rep(1:3, each = 4)
comps_G <- comps |> 
  mutate(Group = "G", .before = Plot4) |> 
  mutate(Week = wk, .before = Group)
df_G <- df1 |> 
  filter(Group == "GS", Week %in% 4:7) |> 
  mutate(`z value` = NA, `Pr(>|z|)` = NA)
df_G$Group[df_G$Group == "GS"] <- "G"
comps_G <- comps_G |> 
  bind_rows(df_G)

df_weekly <- all |> 
  bind_rows(comps_F) |> 
  bind_rows(comps_ES) |> 
  bind_rows(comps_DS) |> 
  bind_rows(comps_G)
df_weekly$`z value`[df_weekly$Plot4 == "C"] <- NA
df_weekly$`Pr(>|z|)`[df_weekly$Plot4 == "C"] <- 1
knitr::kable(df_weekly)

Week	Group	Plot4	response	SE	z value	Pr(>\|z\|)
1	All	C	0.8935823	0.3923476	NA	1.0000000
1	All	y02	1.5778630	0.7595641	1.6023322	0.1090822
1	All	y05	1.4787743	0.8343856	1.0683427	0.2853660
1	All	y27	1.8888779	1.0943633	1.4972195	0.1343362
2	All	C	0.6081992	0.2224435	NA	1.0000000
2	All	y02	1.0082112	0.4111369	1.4589668	0.1445742
2	All	y05	0.2067515	0.1257931	-1.9018629	0.0571891
2	All	y27	3.0249933	1.4474927	3.7270504	0.0001937
3	All	C	0.4847671	0.1552286	NA	1.0000000
3	All	y02	0.6229257	0.2426944	0.6133628	0.5396365
3	All	y05	0.1770807	0.0967597	-1.7789068	0.0752550
3	All	y27	1.8786407	0.9286749	2.6059660	0.0091616
4	All	C	0.4226431	0.1872410	NA	1.0000000
4	All	y02	0.6549094	0.3266899	1.1124526	0.2659436
4	All	y05	0.2474874	0.1531240	-0.9956452	0.3194226
4	All	y27	1.3516976	0.7816569	2.3672690	0.0179199
5	All	C	0.4594750	0.1782016	NA	1.0000000
5	All	y02	0.6148759	0.2795550	0.7333929	0.4633188
5	All	y05	0.4083776	0.2270963	-0.2304923	0.8177093
5	All	y27	1.2233057	0.6846505	1.8519456	0.0640336
6	All	C	0.2634550	0.0713399	NA	1.0000000
6	All	y02	0.2406145	0.0916540	-0.1940415	0.8461434
6	All	y05	0.2823486	0.1460238	0.1186406	0.9055601
6	All	y27	1.5312494	0.7000600	3.3121904	0.0009257
7	All	C	0.0703125	0.0402270	NA	1.0000000
7	All	y02	0.1048338	0.0798231	0.4193847	0.6749350
7	All	y05	0.0000000	0.0000000	-0.0000052	0.9999959
7	All	y27	0.7500000	0.7140762	2.1310601	0.0330842
1	F	C	0.9525069	0.5360405	NA	1.0000000
1	F	y05	0.0000000	0.0000000	-0.0000012	0.9999990
1	F	y27	0.0000000	0.0000000	-0.0000036	0.9999971
2	F	C	0.6614338	0.5565597	NA	1.0000000
2	F	y02	0.4000000	0.4782849	-0.3439867	0.7308563
2	F	y05	0.0000000	0.0000000	-0.0000011	0.9999991
2	F	y27	0.1000000	0.1389147	-1.1632410	0.2447317
3	F	C	0.4313820	0.2947420	NA	1.0000000
3	F	y02	0.0874430	0.0882401	-1.3096353	0.1903192
3	F	y05	0.3637644	0.3401233	-0.1472197	0.8829586
3	F	y27	0.0000000	0.0000000	-0.0000033	0.9999974
4	F	C	0.2614209	0.1327915	NA	1.0000000
4	F	y02	0.1395554	0.1029979	-0.7005618	0.4835765
4	F	y05	0.0365136	0.0452095	-1.4708528	0.1413309
4	F	y27	0.0000000	0.0000000	-0.0000023	0.9999982
5	F	C	0.2099626	0.0744726	NA	1.0000000
5	F	y02	0.0934104	0.0599799	-1.1040988	0.2695503
5	F	y05	0.0735056	0.0581625	-1.2103953	0.2261272
5	F	y27	0.0000000	0.0000000	-0.0000029	0.9999977
6	F	C	0.1654441	0.1022623	NA	1.0000000
6	F	y02	0.1231631	0.1338409	-0.2566322	0.7974627
6	F	y05	0.1860649	0.1166005	0.1459425	0.8839668
7	F	C	0.0000000	0.0000000	NA	1.0000000
7	F	y02	0.0000000	NA	NA	NA
7	F	y05	0.0000000	0.0000000	NA	NA
1	ES	C	0.2598106	0.1298104	NA	1.0000000
1	ES	y02	0.3109820	0.2109119	0.2134211	0.8309985
1	ES	y05	0.3620275	0.4702113	0.2384048	0.8115671
1	ES	y27	1.4913255	1.1946514	1.8509197	0.0641811
2	ES	C	0.8377186	0.5631953	NA	1.0000000
2	ES	y02	1.0714638	0.7611039	0.4267691	0.6695475
2	ES	y05	0.3982378	0.5629115	-0.5436610	0.5866748
2	ES	y27	1.5196373	1.2452829	0.8846264	0.3763582
3	ES	C	0.4759027	0.3552931	NA	1.0000000
3	ES	y02	0.4752883	0.3805364	-0.0020796	0.9983408
3	ES	y05	0.0692212	0.0976584	-1.4565830	0.1452315
3	ES	y27	0.7564644	0.6761801	0.6514532	0.5147540
4	ES	C	0.5856812	0.4224746	NA	1.0000000
4	ES	y02	0.7362194	0.5582987	0.5241442	0.6001783
4	ES	y05	0.2699076	0.2326038	-1.2979940	0.1942894
4	ES	y27	0.4188003	0.3544879	-0.5826463	0.5601314
5	ES	C	0.5991436	0.3148544	NA	1.0000000
5	ES	y02	0.8777375	0.5035568	0.9691896	0.3324506
5	ES	y05	0.5086366	0.3329274	-0.3213867	0.7479174
5	ES	y27	0.0675011	0.0830553	-1.8854622	0.0593675
6	ES	C	0.5581433	0.1298812	NA	1.0000000
6	ES	y02	0.5431968	0.1621329	-0.0717202	0.9428246
6	ES	y05	0.2753500	0.1318402	-1.3272422	0.1844286
6	ES	y27	0.0000000	0.0000000	-0.0000007	0.9999994
7	ES	C	0.2586476	0.1725402	NA	1.0000000
7	ES	y02	0.5132518	0.2805168	0.7946511	0.4268165
7	ES	y05	0.0000000	0.0000000	-0.0000011	0.9999991
1	DS	C	1.1955432	0.6345120	NA	1.0000000
1	DS	y02	2.1604227	1.2215049	1.5810710	0.1138618
1	DS	y05	0.7002072	0.4832633	-1.0134621	0.3108395
1	DS	y27	2.6740970	1.8460733	1.4234650	0.1546014
2	DS	C	0.5295073	0.1419459	NA	1.0000000
2	DS	y02	1.1351536	0.3455672	1.8799774	0.0601112
2	DS	y05	0.3627801	0.2080858	-0.5972629	0.5503319
2	DS	y27	2.8398618	1.2743920	3.2130983	0.0013131
3	DS	C	0.4557638	0.1447556	NA	1.0000000
3	DS	y02	1.1284687	0.4339891	1.8177229	0.0691065
3	DS	y05	0.3166212	0.2268675	-0.4647676	0.6420979
3	DS	y27	2.2735752	1.2109308	2.5916430	0.0095519
4	DS	C	0.2832942	0.1224305	NA	1.0000000
4	DS	y02	0.1418516	0.1613611	-0.5684331	0.5697409
4	DS	y05	0.5404668	0.4067283	0.7443389	0.4566715
4	DS	y27	2.8042269	1.1878920	3.7881185	0.0001518
5	DS	C	0.4503208	0.2202864	NA	1.0000000
5	DS	y02	0.0000000	0.0000000	-0.0000036	0.9999972
5	DS	y05	1.5000000	1.3222594	1.1935442	0.2326563
5	DS	y27	3.7330067	1.7821295	3.0942804	0.0019729
6	DS	C	0.2288253	0.1351041	NA	1.0000000
6	DS	y02	0.1905523	0.2412042	-0.1310421	0.8957420
6	DS	y05	3.0000046	4.1313286	1.7175057	0.0858868
6	DS	y27	3.4999967	1.8102650	3.4749012	0.0005110
7	DS	C	0.5000000	0.2900000	NA	1.0000000
7	DS	y27	2.0000000	1.3200000	NA	NA
1	G	C	0.3192735	0.1423014	NA	1.0000000
1	G	y02	0.2102387	0.1408601	-0.5192016	0.6036202
1	G	y05	1.5394615	1.2083052	1.7428822	0.0813542
1	G	y27	1.2500052	2.7664492	0.6045643	0.5454685
2	G	C	0.1019134	0.0626354	NA	1.0000000
2	G	y02	0.1267690	0.0963975	0.2232134	0.8233694
2	G	y05	0.1428573	0.1539945	0.2721690	0.7854921
2	G	y27	1.2593231	0.8773174	2.7063309	0.0068031
3	G	C	0.1943292	0.1502488	NA	1.0000000
3	G	y02	0.1263586	0.1704571	-0.2768292	0.7819112
3	G	y05	0.1007669	0.1734028	-0.3481200	0.7277500
3	G	y27	0.0000000	0.0000000	-0.0000009	0.9999993
4	G	C	0.1100000	0.0700000	NA	1.0000000
4	G	y02	0.0000000	0.0000000	NA	NA
4	G	y05	0.0000000	0.0000000	NA	NA
4	G	y27	0.0000000	0.0000000	NA	NA
5	G	C	0.1000000	0.0700000	NA	1.0000000
5	G	y02	0.0000000	0.0000000	NA	NA
5	G	y05	0.5000000	NA	NA	NA
5	G	y27	0.0000000	0.0000000	NA	NA
6	G	C	0.0000000	0.0000000	NA	1.0000000
6	G	y02	0.0000000	0.0000000	NA	NA
6	G	y27	0.0000000	NA	NA	NA
7	G	C	3.0000000	NA	NA	1.0000000
7	G	y02	0.0000000	NA	NA	NA
7	G	y05	0.0000000	NA	NA	NA

Weekly plot by growth forms:

df <- df_weekly |> 
  mutate(Plot4 = factor(Plot4), Group = factor(Group))

df$sig <- map_vec(df$`Pr(>|z|)`, f_star)

# Lower and upper bounds for Mean +- SE
df$Lower <- df$response - df$SE
df$Upper <- df$response + df$SE

# Plot groups
df$paired <- rep(NA, nrow(df))
for (i in 1:nrow(df)) {
  df$paired[i] <- paste0(df$Plot4[i], "_", df$Group[i])
}

# Position of stars
df$star_pos <- rep(NA, nrow(df))
df$star_pos[df$Group == "F"] <- 1
df$star_pos[df$Group == "ES"] <- 2.5
df$star_pos[df$Group == "DS"] <- 6.5
df$star_pos[df$Group == "G"] <- 3

levels(df$Plot4) <- list(C  = "C", `W (1 yr)` = "y02", `W (4 yrs)` = "y05", `W (26 yrs)` = "y27")
levels(df$Group) <- list(All = "All", F  = "F", ES = "ES", DS = "DS", G = "G")
names(df)[2] <- "Growth form"
names(df)[3] <- "Plot treatment"

p <- df %>%
  ggplot(aes(x = Week, y = response, group = `Plot treatment`, color = `Plot treatment`)) + 
  ylab(as.expression(bquote("Floral density (flowers per 30.5 x 30.5" ~ cm^2 ~ ")"))) + 
  geom_point(data = df, mapping = aes(x = Week, y = star_pos + 0.5), color = "white", alpha = 0) +
  geom_point(position = position_dodge(width = 0.1)) +
  geom_line(aes(group = paired, color = `Plot treatment`)) +
  geom_errorbar(aes(ymin = Lower, ymax = Upper), width = 0, 
                position = position_dodge(width = 0.1)) +
  geom_text(aes(label = sig, y = star_pos), angle = 90, position = position_dodge(width = 0.15), vjust = 0.6, show.legend = FALSE) + facet_grid(rows = vars(`Growth form`), scales = "free_y") +
  theme_bw() 
p

# Auxiliary function to look at floral density at the species level
f_gf <- function(gf) {
  # Weekly summaries
  ws <- phenology2019 |> 
    group_by(Week, PlotU, Plot4, Species) |> 
    summarize(
      Group = Group[1],
      fd_mean = f_mean(Floral.density),
      fd_sum = f_sum(Floral.density),
      n = n()
    ) |> 
    filter(!is.na(PlotU), Group == gf) |> 
    ungroup() |> 
    mutate(PlotU = factor(PlotU))
  
  # Means and standard errors
  df <- ws |> 
    group_by(Week, Plot4, Species) |> 
    summarize(Mean = round(mean(fd_mean, na.rm = T), 2),
              SD = round(sd(fd_mean, na.rm = T), 2),
              SE = round(sd(fd_mean, na.rm = T)/sqrt(sum(!is.na(fd_mean))), 2))
  
  # Visualization
  p <- df |> 
    ggplot(aes(y = Mean, x = Week, color = Plot4)) +
    geom_line() +
    geom_point() +
    facet_wrap(vars(Species)) +
    theme_bw()
  return(p)
}

By forbs species:

f_gf("F")

By evergreen shrubs species:

f_gf("ES")

By deciduous shrubs species:

f_gf("DS")

By graminoid species:

f_gf("GS")

Number of species

2019

Checking for under/over dispersion:

m <- glmer(n_sp ~ Plot4 + (1 | Location/Site) , family = poisson, data = plot_summaries_2019)
deviance(m)/df.residual(m)

## [1] 0.5143948

There’s underdispersion. will use a quasi-poisson adjustment.

Model selection:

m1 <- glmer(n_sp ~ Plot4 + (1 | Location/Site), family = poisson, data = plot_summaries_2019)
m2 <- glmer(n_sp ~ Plot4 + (1 | Site) , family = poisson, data = plot_summaries_2019)
m3 <- glmer(n_sp ~ Plot4 + (1 | Location) , family = poisson, data = plot_summaries_2019)
anova(m1, m2)

## Data: plot_summaries_2019
## Models:
## m2: n_sp ~ Plot4 + (1 | Site)
## m1: n_sp ~ Plot4 + (1 | Location/Site)
##    npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## m2    5 516.72 530.98 -253.36   506.72                    
## m1    6 519.72 536.83 -253.86   507.72     0  1          1

anova(m1, m3)

## Data: plot_summaries_2019
## Models:
## m3: n_sp ~ Plot4 + (1 | Location)
## m1: n_sp ~ Plot4 + (1 | Location/Site)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## m3    5 524.28 538.54 -257.14   514.28                       
## m1    6 519.72 536.83 -253.86   507.72 6.5593  1    0.01043 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(m2, m3)

## Data: plot_summaries_2019
## Models:
## m2: n_sp ~ Plot4 + (1 | Site)
## m3: n_sp ~ Plot4 + (1 | Location)
##    npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## m2    5 516.72 530.98 -253.36   506.72                    
## m3    5 524.28 538.54 -257.14   514.28     0  0

Nesting Site within Location as a random effect does not improve a model with only Site (community type) as a random effect. Will only include community type as a random effect.

# p-values
m <- glmer(n_sp ~ Plot4 + (1 | Site), family = poisson, data = plot_summaries_2019)
emms <- as_tibble(emmeans(m, ~ Plot4, type = "response"))
s <- f_qp(m) # quasi-poisson adjustment
all <- bind_cols(emms[,1:3], s[,3:4])
knitr::kable(all, digits = 4)

Plot4	rate	SE	z value	Pr(>\|z\|)
C	4.9088	0.5423	19.8231	0.0000
y02	5.0482	0.6220	0.4001	0.6891
y05	5.6985	0.8041	1.7348	0.0828
y27	3.9642	0.6225	-2.1479	0.0317

Fill in for plots with no species in bloom:

df0 <- expand_grid(PlotU = factor(1:128), Group = c("F", "ES", "DS", "GS")) |>
  mutate(Plot4 = factor(if_else(
    PlotU %in% c(1:16, 33:48, 65:80, 97:112), "C", 
    if_else(PlotU %in% c(17:24, 49:56), "y05", 
            if_else(PlotU %in% c(25:32, 57:64, 89:96, 121:128), "y02",
                    if_else(PlotU %in% c(81:88, 113:120), "y27", NA))))),
    Location = if_else(PlotU %in% 1:64, "Imnavait", "Toolik"),
    Site = if_else(PlotU %in% c(1:32, 65:96), "Dry", "Moist"))
df1 <- plot_summaries_GF_2019 |> 
  right_join(df0) |> 
  mutate(n_sp = replace(n_sp, is.na(n_sp), 0)) 
df2 <- df1 |> 
  group_by(Group, Plot4) |> 
  summarize(Mean = round(mean(n_sp, na.rm = T), 2),
            SD = round(sd(n_sp, na.rm = T), 2))

# auxiliary function to compute EMMs, SEs, and p-values
f <- function(gf){
  m <- glmer(n_sp ~ Plot4 + (1 | Site), family = poisson, data = filter(df1, Group == gf))
  emms <- as_tibble(emmeans(m, ~ Plot4, type = "response"))
  s <- f_qp(m) # quasi-poisson adjustment
  return(bind_cols(emms[,1:3], s[,3:4]))
}

# p-values
comps <- map_df(c("F", "ES", "DS", "GS"), f) |> 
  mutate(Group = rep(c("F", "ES", "DS", "G"), each = 4), .before = Plot4)
knitr::kable(comps, digits = 4)

Group	Plot4	rate	SE	z value	Pr(>\|z\|)
F	C	0.8686	0.2148	-0.5858	0.5580
F	y02	0.8686	0.2431	0.0000	1.0000
F	y05	1.6174	0.4599	2.7528	0.0059
F	y27	0.4792	0.1966	-1.6268	0.1038
ES	C	1.9357	0.3069	6.0827	0.0000
ES	y02	1.9664	0.3562	0.1502	0.8806
ES	y05	1.8435	0.4142	-0.3512	0.7255
ES	y27	1.5977	0.3767	-1.3028	0.1926
DS	C	1.4063	0.1482	4.5081	0.0000
DS	y02	1.5000	0.2165	0.5033	0.6147
DS	y05	1.5000	0.3062	0.3916	0.6954
DS	y27	1.5625	0.3125	0.6496	0.5160
G	C	0.6686	0.1126	-2.9976	0.0027
G	y02	0.6842	0.1538	0.1100	0.9124
G	y05	0.6842	0.2119	0.0854	0.9319
G	y27	0.3110	0.1408	-2.0326	0.0421

Combining the data for plotting:

ns <- all |> 
  mutate(Group = rep("All", 4), .before = Plot4) |> 
  bind_rows(comps) |> 
  mutate(Lower = rate - SE, 
         Upper = rate + SE,
         Group = factor(Group, levels = c("All", "F", "ES", "DS", "G")))
ns$`z value`[ns$Plot4 == "C"] <- NA
ns$`Pr(>|z|)`[ns$Plot4 == "C"] <- 1
ns$sig <- map_vec(ns$`Pr(>|z|)`, f_star)
levels(ns$Plot4) <- list(C  = "C", `W (1 yr)` = "y02", `W (4 yrs)` = "y05", `W (26 yrs)` = "y27")
names(ns)[1:2] <- c("Growth form", "Plot treatment")

Plot:

p <- ns |> 
  ggplot(aes(x = `Growth form`, y = rate, color = `Plot treatment`)) + 
  ylab("Number of species per plot") + 
  geom_pointrange(aes(ymin = Lower, ymax = Upper), position = position_dodge(width = 0.3)) +
  geom_text(aes(label = sig, x = `Growth form`, y = Upper + 0.17, group = `Plot treatment`), angle = 90, position = position_dodge(width = 0.34), vjust = 0.5, show.legend = FALSE) +
  theme_bw()
p

Weekly, all species:

df0 <- expand_grid(PlotU = factor(1:128), Week = 1:8) |> 
  mutate(Plot4 = factor(if_else(
    PlotU %in% c(1:16, 33:48, 65:80, 97:112), "C", 
    if_else(PlotU %in% c(17:24, 49:56), "y05", 
            if_else(PlotU %in% c(25:32, 57:64, 89:96, 121:128), "y02",
                    if_else(PlotU %in% c(81:88, 113:120), "y27", NA))))),
    Location = if_else(PlotU %in% 1:64, "Imnavait", "Toolik"),
    Site = if_else(PlotU %in% c(1:32, 65:96), "Dry", "Moist"))
df1 <- weekly_summaries_2019 |> 
  right_join(df0) |> 
  mutate(n_sp = replace(n_sp, is.na(n_sp), 0))

f <- function(w){
  m <- glmer(n_sp ~ Plot4 + (1 | Site), family = poisson, data = filter(df1, Week == w))
  emms <- as_tibble(emmeans(m, ~ Plot4, type = "response"))
  s <- f_qp(m) # quasi-poisson adjustment
  return(bind_cols(emms[,1:3], s[,3:4]))
}

# p-values
all <- map_df(1:8, f) |> 
  mutate(Week = rep(1:8, each = 4), .before = Plot4) |> 
  mutate(Group = "All", .before = Plot4)
knitr::kable(all, digits = 4)

Week	Group	Plot4	rate	SE	z value	Pr(>\|z\|)
1	All	C	2.2031	0.1855	13.7833	0.0000
1	All	y02	2.0625	0.2539	-0.6499	0.5158
1	All	y05	2.3125	0.3802	0.3855	0.6999
1	All	y27	2.0625	0.3590	-0.5012	0.6162
2	All	C	1.9507	0.2561	7.1599	0.0000
2	All	y02	2.4152	0.3591	2.0873	0.0369
2	All	y05	2.0436	0.4065	0.3350	0.7376
2	All	y27	2.4152	0.4512	1.6411	0.1008
3	All	C	2.0614	0.4038	4.8656	0.0000
3	All	y02	2.3948	0.5004	1.3984	0.1620
3	All	y05	2.7283	0.6293	2.1504	0.0315
3	All	y27	2.4857	0.5849	1.3863	0.1657
4	All	C	2.2301	0.4282	5.7405	0.0000
4	All	y02	2.1240	0.4471	-0.4622	0.6439
4	All	y05	2.6701	0.6132	1.4417	0.1494
4	All	y27	2.1846	0.5251	-0.1525	0.8788
5	All	C	1.8747	0.2975	5.0491	0.0000
5	All	y02	1.5366	0.2954	-1.5119	0.1306
5	All	y05	2.3971	0.4946	1.7063	0.0880
5	All	y27	1.2293	0.3179	-2.2333	0.0255
6	All	C	0.9844	0.1240	-0.1518	0.8794
6	All	y02	1.0938	0.1849	0.6068	0.5440
6	All	y05	1.5000	0.3062	2.1320	0.0330
6	All	y27	0.5625	0.1875	-1.9067	0.0566
7	All	C	0.2804	0.0683	-5.7613	0.0000
7	All	y02	0.3738	0.1103	0.8522	0.3941
7	All	y05	0.4984	0.1789	1.4949	0.1349
7	All	y27	0.3738	0.1543	0.6737	0.5005
8	All	C	0.0000	0.0000	0.0000	1.0000
8	All	y02	0.0000	0.0000	0.0000	1.0000
8	All	y05	0.0625	0.0625	0.0000	1.0000
8	All	y27	0.0000	0.0000	-0.0011	0.9992

Weekly, by growth forms:

df0 <- expand_grid(PlotU = factor(1:128), Week = 1:8, Group = c("F", "ES", "DS", "GS")) |> 
  mutate(Plot4 = factor(if_else(
    PlotU %in% c(1:16, 33:48, 65:80, 97:112), "C", 
    if_else(PlotU %in% c(17:24, 49:56), "y05", 
            if_else(PlotU %in% c(25:32, 57:64, 89:96, 121:128), "y02",
                    if_else(PlotU %in% c(81:88, 113:120), "y27", NA))))),
    Location = if_else(PlotU %in% 1:64, "Imnavait", "Toolik"),
    Site = if_else(PlotU %in% c(1:32, 65:96), "Dry", "Moist"))
df1 <- weekly_summaries_GF_2019 |> 
  right_join(df0) |> 
  mutate(n_sp = replace(n_sp, is.na(n_sp), 0))

# auxiliary function to compute EMMs, SEs, and p-values
f <- function(w, gf){
  m <- glmer(n_sp ~ Plot4 + (1 | Site), family = poisson, data = filter(df1, Week == w, Group == gf))
  emms <- as_tibble(emmeans(m, ~ Plot4, type = "response"))
  s <- f_qp(m)
  return(bind_cols(emms[,1:3], s[,3:4]))
}

# Forbs
comps <- map_df(1:8, \(x)f(x, "F")) 
wk <- rep(1:8, each = 4)
comps_F <- comps |> 
  mutate(Group = "F", .before = Plot4) |> 
  mutate(Week = wk, .before = Group)

# Evergreen shrubs
comps <- map_df(1:7, \(x)f(x, "ES")) # week 8 not enough obs
wk <- rep(1:7, each = 4)
comps_ES <- comps |> 
  mutate(Group = "ES", .before = Plot4) |> 
  mutate(Week = wk, .before = Group)
aux <- tibble(Week = 8, Group = "ES", Plot4 = c("C", "y02", "y05", "y27"), rate = numeric(4), SE = numeric(4), `z value` = rep(NA, 4), `Pr(>|z|)` = rep(NA, 4))
comps_ES <- bind_rows(comps_ES, aux)

# Deciduous shrubs
comps <- map_df(1:7, \(x)f(x, "DS")) # week 8 not enough obs
wk <- rep(1:7, each = 4)
comps_DS <- comps |> 
  mutate(Group = "DS", .before = Plot4) |> 
  mutate(Week = wk, .before = Group)
aux <- tibble(Week = 8, Group = "DS", Plot4 = c("C", "y02", "y05", "y27"), rate = numeric(4), SE = numeric(4), `z value` = rep(NA, 4), `Pr(>|z|)` = rep(NA, 4))
comps_DS <- bind_rows(comps_DS, aux)
comps_DS$`Pr(>|z|)`[comps_DS$Plot4 == "y27" & comps_DS$Week == 7] <- 1 # method did not converge for week 7

# Graminoid
comps <- map_df(1:7, \(x)f(x, "GS")) # week 8 not enough obs
wk <- rep(1:7, each = 4)
comps_G <- comps |> 
  mutate(Group = "G", .before = Plot4) |> 
  mutate(Week = wk, .before = Group)
aux <- tibble(Week = 8, Group = "G", Plot4 = c("C", "y02", "y05", "y27"), rate = numeric(4), SE = numeric(4), `z value` = rep(NA, 4), `Pr(>|z|)` = rep(NA, 4))
comps_G <- bind_rows(comps_G, aux)

nsp <- all |> 
  bind_rows(comps_F) |> 
  bind_rows(comps_ES) |> 
  bind_rows(comps_DS) |>
  bind_rows(comps_G)
nsp$`z value`[nsp$Plot4 == "C"] <- 0
nsp$`Pr(>|z|)`[nsp$Plot4 == "C"] <- 1
knitr::kable(nsp)

Week	Group	Plot4	rate	SE	z value	Pr(>\|z\|)
1	All	C	2.2031250	0.1855366	0.0000000	1.0000000
1	All	y02	2.0625000	0.2538762	-0.6499028	0.5157550
1	All	y05	2.3125000	0.3801727	0.3854786	0.6998829
1	All	y27	2.0625000	0.3590353	-0.5012377	0.6162038
2	All	C	1.9507362	0.2560986	0.0000000	1.0000000
2	All	y02	2.4152049	0.3591283	2.0872940	0.0368616
2	All	y05	2.0436460	0.4064749	0.3350242	0.7376069
2	All	y27	2.4151902	0.4512399	1.6410740	0.1007821
3	All	C	2.0613510	0.4037716	0.0000000	1.0000000
3	All	y02	2.3948082	0.5003885	1.3984035	0.1619919
3	All	y05	2.7282754	0.6292638	2.1504018	0.0315234
3	All	y27	2.4857487	0.5848623	1.3862954	0.1656567
4	All	C	2.2301378	0.4281538	0.0000000	1.0000000
4	All	y02	2.1239588	0.4471483	-0.4621955	0.6439411
4	All	y05	2.6701179	0.6131977	1.4417010	0.1493868
4	All	y27	2.1846295	0.5251164	-0.1525466	0.8787559
5	All	C	1.8746693	0.2974949	0.0000000	1.0000000
5	All	y02	1.5366103	0.2953909	-1.5119339	0.1305507
5	All	y05	2.3971194	0.4946443	1.7062538	0.0879608
5	All	y27	1.2292950	0.3178930	-2.2333396	0.0255266
6	All	C	0.9843750	0.1240196	0.0000000	1.0000000
6	All	y02	1.0937500	0.1848775	0.6067915	0.5439893
6	All	y05	1.5000000	0.3061862	2.1320082	0.0330062
6	All	y27	0.5625000	0.1875002	-1.9067123	0.0565578
7	All	C	0.2803651	0.0683081	0.0000000	1.0000000
7	All	y02	0.3738174	0.1103423	0.8521928	0.3941071
7	All	y05	0.4984237	0.1788655	1.4948684	0.1349488
7	All	y27	0.3738207	0.1543267	0.6737409	0.5004761
8	All	C	0.0000000	0.0000000	0.0000000	1.0000000
8	All	y02	0.0000000	0.0000000	0.0000004	0.9999997
8	All	y05	0.0625000	0.0625000	0.0000118	0.9999906
8	All	y27	0.0000000	0.0000000	-0.0010514	0.9991611
1	F	C	0.1221660	0.0995865	0.0000000	1.0000000
1	F	y02	0.0000000	0.0000000	-0.1751009	0.8610004
1	F	y05	0.1221661	0.1162670	0.0000002	0.9999998
1	F	y27	0.2036101	0.1794884	1.1646791	0.2441489
2	F	C	0.1562500	0.0494105	0.0000000	1.0000000
2	F	y02	0.1562500	0.0698771	0.0000000	1.0000000
2	F	y05	0.2500000	0.1250001	0.7739575	0.4389559
2	F	y27	0.3125000	0.1397543	1.2328638	0.2176266
3	F	C	0.2350826	0.1191705	0.0000000	1.0000000
3	F	y02	0.3918047	0.2028033	1.5479318	0.1216387
3	F	y05	0.7313670	0.3818657	3.3741642	0.0007404
3	F	y27	0.2612022	0.1650497	0.2207882	0.8252574
4	F	C	0.4065915	0.1767089	0.0000000	1.0000000
4	F	y02	0.5963414	0.2672721	1.6012021	0.1093322
4	F	y05	0.8674052	0.4047407	2.8721670	0.0040767
4	F	y27	0.1626356	0.1133245	-1.7757075	0.0757811
5	F	C	0.4800132	0.1565493	0.0000000	1.0000000
5	F	y02	0.5527411	0.1981233	0.4838108	0.6285202
5	F	y05	0.9891123	0.3628979	2.3917746	0.0167671
5	F	y27	0.0581830	0.0600944	-2.0531441	0.0400586
6	F	C	0.2031250	0.0563367	0.0000000	1.0000000
6	F	y02	0.2187500	0.0826797	0.1452177	0.8845390
6	F	y05	0.5625000	0.1875000	2.1578489	0.0309396
6	F	y27	0.0000000	0.0000000	-0.0000130	0.9999896
7	F	C	0.0312500	0.0220971	0.0000000	1.0000000
7	F	y02	0.0312500	0.0312500	-0.0000001	0.9999999
7	F	y05	0.3125000	0.1397542	2.9049234	0.0036734
7	F	y27	0.0000000	0.0000000	-0.0000042	0.9999966
8	F	C	0.0000000	0.0000000	0.0000000	1.0000000
8	F	y02	0.0000000	0.0000000	0.0000007	0.9999994
8	F	y05	0.0625000	0.0625000	0.0000131	0.9999896
8	F	y27	0.0000000	0.0000000	-0.0014699	0.9988272
1	ES	C	0.4687500	0.0855817	0.0000000	1.0000000
1	ES	y02	0.4687500	0.1210308	0.0000000	1.0000000
1	ES	y05	0.2500000	0.1250002	-1.4296045	0.1528306
1	ES	y27	0.6250000	0.1976423	0.9537375	0.3402165
2	ES	C	0.5156250	0.0897587	0.0000000	1.0000000
2	ES	y02	0.8125000	0.1593443	2.0816730	0.0373723
2	ES	y05	0.3125000	0.1397541	-1.2526469	0.2103342
2	ES	y27	0.9375000	0.2420614	2.3046237	0.0211877
3	ES	C	0.8401041	0.1711001	0.0000000	1.0000000
3	ES	y02	0.9775735	0.2281971	0.8225988	0.4107362
3	ES	y05	0.9775731	0.2860196	0.6438757	0.5196561
3	ES	y27	1.1608723	0.3195635	1.4666486	0.1424717
4	ES	C	1.2308128	0.2490531	0.0000000	1.0000000
4	ES	y02	1.1852371	0.2758955	-0.2514820	0.8014415
4	ES	y05	1.2764026	0.3520467	0.1929382	0.8470074
4	ES	y27	1.0940504	0.3170410	-0.5871833	0.5570806
5	ES	C	0.9861780	0.1705044	0.0000000	1.0000000
5	ES	y02	0.8012774	0.1838653	-1.3985042	0.1619617
5	ES	y05	1.1710795	0.3026258	1.0303388	0.3028510
5	ES	y27	0.6163624	0.2081010	-2.1651719	0.0303745
6	ES	C	0.5156250	0.0897588	0.0000000	1.0000000
6	ES	y02	0.6875000	0.1465755	1.4471218	0.1478628
6	ES	y05	0.8750000	0.2338536	2.2956455	0.0216962
6	ES	y27	0.0625000	0.0624997	-2.8783998	0.0039970
7	ES	C	0.1718750	0.0518223	0.0000000	1.0000000
7	ES	y02	0.3125000	0.0988212	1.6085219	0.1077209
7	ES	y05	0.1250000	0.0883883	-0.4870166	0.6262466
7	ES	y27	0.0000000	0.0000000	-0.0000050	0.9999960
8	ES	C	0.0000000	0.0000000	0.0000000	1.0000000
8	ES	y02	0.0000000	0.0000000	NA	NA
8	ES	y05	0.0000000	0.0000000	NA	NA
8	ES	y27	0.0000000	0.0000000	NA	NA
1	DS	C	1.0780460	0.1303359	0.0000000	1.0000000
1	DS	y02	1.1249181	0.1879045	0.3250744	0.7451248
1	DS	y05	1.3748969	0.2935291	1.5600130	0.1187568
1	DS	y27	0.9999268	0.2502299	-0.4257216	0.6703107
2	DS	C	0.8286372	0.1958048	0.0000000	1.0000000
2	DS	y02	0.9642221	0.2529744	0.8379075	0.4020826
2	DS	y05	1.0244989	0.3178520	0.9399665	0.3472347
2	DS	y27	0.9039487	0.2915955	0.3671296	0.7135223
3	DS	C	0.6591588	0.1978246	0.0000000	1.0000000
3	DS	y02	0.7616904	0.2482070	0.6541014	0.5130465
3	DS	y05	0.6445088	0.2563810	-0.0744755	0.9406320
3	DS	y27	0.8202858	0.3057205	0.7964459	0.4257729
4	DS	C	0.3534621	0.1146554	0.0000000	1.0000000
4	DS	y02	0.1767324	0.0845368	-1.4848308	0.1375886
4	DS	y05	0.2945501	0.1506989	-0.3626408	0.7168733
4	DS	y27	0.6480110	0.2542202	1.6276702	0.1035948
5	DS	C	0.1632794	0.0807957	0.0000000	1.0000000
5	DS	y02	0.0544277	0.0439357	-1.6021444	0.1091237
5	DS	y05	0.1088540	0.0878709	-0.5912961	0.5543220
5	DS	y27	0.3809824	0.2097027	1.9843396	0.0472180
6	DS	C	0.1137931	0.0998467	0.0000000	1.0000000
6	DS	y02	0.0568961	0.0572333	-1.0746264	0.2825420
6	DS	y05	0.0379305	0.0488298	-1.0562711	0.2908444
6	DS	y27	0.2655156	0.2415338	1.7828626	0.0746087
7	DS	C	0.0146536	0.0270795	0.0000000	1.0000000
7	DS	y02	0.0000000	0.0000001	-0.0007658	0.9993890
7	DS	y05	0.0000000	0.0000004	-0.0011419	0.9990889
7	DS	y27	0.0879219	0.1604557	5.4420645	1.0000000
8	DS	C	0.0000000	0.0000000	0.0000000	1.0000000
8	DS	y02	0.0000000	0.0000000	NA	NA
8	DS	y05	0.0000000	0.0000000	NA	NA
8	DS	y27	0.0000000	0.0000000	NA	NA
1	G	C	0.4068090	0.1772990	0.0000000	1.0000000
1	G	y02	0.4068108	0.1920328	0.0000154	0.9999877
1	G	y05	0.4339252	0.2297616	0.1730569	0.8626067
1	G	y27	0.1084842	0.0874796	-1.9314113	0.0534322
2	G	C	0.4347706	0.0896535	0.0000000	1.0000000
2	G	y02	0.4658280	0.1262366	0.2566563	0.7974441
2	G	y05	0.4347639	0.1681252	-0.0000434	0.9999654
2	G	y27	0.2484453	0.1258390	-1.2460617	0.2127417
3	G	C	0.2314888	0.1231876	0.0000000	1.0000000
3	G	y02	0.1543256	0.0966040	-0.8961438	0.3701760
3	G	y05	0.2057700	0.1415685	-0.2219812	0.8243285
3	G	y27	0.1543266	0.1149533	-0.6774081	0.4981471
4	G	C	0.0207459	0.0635616	0.0000000	1.0000000
4	G	y02	0.0103729	0.0322012	-1.6970411	0.0896889
4	G	y05	0.0138306	0.0433046	-0.8389845	0.4014780
4	G	y27	0.0276612	0.0854977	0.7874734	0.4310048
5	G	C	0.0248652	0.0735755	0.0000000	1.0000000
5	G	y02	0.0114762	0.0344515	-1.8751083	0.0607778
5	G	y05	0.0076508	0.0238019	-1.7642528	0.0776894
5	G	y27	0.0153017	0.0463582	-0.9928933	0.3207619
6	G	C	0.0595071	0.0439893	0.0000000	1.0000000
6	G	y02	0.0714080	0.0590039	0.3023629	0.7623754
6	G	y05	0.0000000	0.0000151	-0.0061999	0.9950532
6	G	y27	0.0476057	0.0556041	-0.2467203	0.8051247
7	G	C	0.0156250	0.0156250	0.0000000	1.0000000
7	G	y02	0.0312500	0.0312500	0.5206544	0.6026075
7	G	y05	0.0625000	0.0625000	1.0413093	0.2977320
7	G	y27	0.0000000	0.0000000	-0.0000149	0.9999881
8	G	C	0.0000000	0.0000000	0.0000000	1.0000000
8	G	y02	0.0000000	0.0000000	NA	NA
8	G	y05	0.0000000	0.0000000	NA	NA
8	G	y27	0.0000000	0.0000000	NA	NA

Weekly plot by growth forms:

df <- nsp |> 
  mutate(Plot4 = factor(Plot4), Group = factor(Group))

df$sig <- map_vec(df$`Pr(>|z|)`, f_star)

# Lower and upper bounds for Mean +- SE
df$Lower <- df$rate - df$SE
df$Upper <- df$rate + df$SE

# Plot groups
df$paired <- rep(NA, nrow(df))
for (i in 1:nrow(df)) {
  df$paired[i] <- paste0(df$Plot4[i], "_", df$Group[i])
}

# Position of stars
df$star_pos <- rep(NA, nrow(df))
df$star_pos[df$Group == "All"] <- 3.5
df$star_pos[df$Group == "F"] <- 1.5
df$star_pos[df$Group == "ES"] <- 1.6
df$star_pos[df$Group == "DS"] <- 1.5
df$star_pos[df$Group == "G"] <- 0.5

levels(df$Plot4) <- list(C  = "C", `W (1 yr)` = "y02", `W (4 yrs)` = "y05", `W (26 yrs)` = "y27")
levels(df$Group) <- list(All = "All", F  = "F", ES = "ES", DS = "DS", G = "G")
names(df)[2] <- "Growth form"
names(df)[3] <- "Plot treatment"

p <- df %>%
  ggplot(aes(x = Week, y = rate, group = `Plot treatment`, color = `Plot treatment`)) + 
  ylab(as.expression("Number of species per plot")) + 
  geom_point(mapping = aes(x = Week, y = star_pos + 0.25), color = "white", alpha = 0) +
  geom_point(position = position_dodge(width = 0.1)) +
  geom_line(aes(group = paired, color = `Plot treatment`)) +
  geom_errorbar(aes(ymin = Lower, ymax = Upper), width = 0, 
                position = position_dodge(width = 0.1)) +
  geom_text(aes(label = sig, y = star_pos), angle = 90, position = position_dodge(width = 0.15), vjust = 0.6, show.legend = FALSE) +
  facet_grid(rows = vars(`Growth form`), scales = "free_y") +
  theme_bw() 
p

2020

Model selection:

m1 <- glmer(n_sp ~ Plot4 + (1 | Site) + (1 | Location) , family = poisson, data = plot_summaries_2020)
m2 <- glmer(n_sp ~ Plot4 + (1 | Site) , family = poisson, data = plot_summaries_2020)
m3 <- glmer(n_sp ~ Plot4 + (1 | Location) , family = poisson, data = plot_summaries_2020)
anova(m1, m2)

## Data: plot_summaries_2020
## Models:
## m2: n_sp ~ Plot4 + (1 | Site)
## m1: n_sp ~ Plot4 + (1 | Site) + (1 | Location)
##    npar   AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## m2    4 103.2 107.74 -47.601   95.201                    
## m1    5 105.2 110.88 -47.601   95.201     0  1          1

anova(m1, m3)

## Data: plot_summaries_2020
## Models:
## m3: n_sp ~ Plot4 + (1 | Location)
## m1: n_sp ~ Plot4 + (1 | Site) + (1 | Location)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## m3    4 103.47 108.01 -47.734   95.467                     
## m1    5 105.20 110.88 -47.601   95.201 0.2661  1      0.606

Adding Location as a random effect does not improve a model with Site (community type) as a random effect. Will only include community type as a random effect.

# p-values
m <- glmer(n_sp ~ Plot4 + (1 | Site), family = poisson, data = plot_summaries_2020)
emms <- as_tibble(emmeans(m, ~ Plot4, type = "response"))
s <- f_qp(m) # quasi-poisson adjustment
all <- bind_cols(emms[,1:3], s[,3:4])
knitr::kable(all, digits = 4)

Plot4	rate	SE	z value	Pr(>\|z\|)
C	4.2736	0.7041	9.8066	0.0000
y05	4.8051	0.9672	0.5520	0.5810
y27	4.6394	0.9479	0.3825	0.7021

Fill in for plots with no species in bloom:

df0 <- expand_grid(PlotU = factor(plot_summaries_2020$PlotU), Group = c("F", "ES", "DS", "GS")) |>
  mutate(Plot4 = factor(if_else(
    PlotU %in% c(1:16, 33:48, 65:80, 97:112), "C", 
    if_else(PlotU %in% c(17:24, 49:56), "y05", 
            if_else(PlotU %in% c(25:32, 57:64, 89:96, 121:128), "y02",
                    if_else(PlotU %in% c(81:88, 113:120), "y27", NA))))),
    Location = if_else(PlotU %in% 1:64, "Imnavait", "Toolik"),
    Site = if_else(PlotU %in% c(1:32, 65:96), "Dry", "Moist"))
df1 <- plot_summaries_GF_2020 |> 
  right_join(df0) |> 
  mutate(n_sp = replace(n_sp, is.na(n_sp), 0)) 
df2 <- df1 |> 
  group_by(Group, Plot4) |> 
  summarize(Mean = round(mean(n_sp, na.rm = T), 2),
            SD = round(sd(n_sp, na.rm = T), 2))

# auxiliary function to compute EMMs, SEs, and p-values
f <- function(gf){
  m <- glmer(n_sp ~ Plot4 + (1 | Site), family = poisson, data = filter(df1, Group == gf))
  emms <- as_tibble(emmeans(m, ~ Plot4, type = "response"))
  s <- f_qp(m) # quasi-poisson adjustment
  return(bind_cols(emms[,1:3], s[,3:4]))
}

# p-values
comps <- map_df(c("F", "ES", "DS", "GS"), f) |> 
  mutate(Group = rep(c("F", "ES", "DS", "G"), each = 3), .before = Plot4)
knitr::kable(comps, digits = 4)

Group	Plot4	rate	SE	z value	Pr(>\|z\|)
F	C	0.7070	0.2965	-0.9394	0.3475
F	y05	1.2727	0.5401	1.3417	0.1797
F	y27	0.1591	0.1626	-1.6068	0.1081
ES	C	2.0000	0.4264	4.3453	0.0000
ES	y05	1.8333	0.5528	-0.3149	0.7528
ES	y27	2.8333	0.6872	1.4416	0.1494
DS	C	1.0909	0.3149	0.3567	0.7213
DS	y05	0.8333	0.3727	-0.5988	0.5493
DS	y27	1.6667	0.5270	1.1714	0.2414
G	C	0.4545	0.2033	-2.2761	0.0228
G	y05	0.8333	0.3727	1.2373	0.2160
G	y27	0.0000	0.0000	0.0000	1.0000

Combining the data for plotting:

ns <- all |> 
  mutate(Group = rep("All", 3), .before = Plot4) |> 
  bind_rows(comps) |> 
  mutate(Lower = rate - SE, 
         Upper = rate + SE,
         Group = factor(Group, levels = c("All", "F", "ES", "DS", "G")))
ns$`z value`[ns$Plot4 == "C"] <- NA
ns$`Pr(>|z|)`[ns$Plot4 == "C"] <- 1
ns$sig <- map_vec(ns$`Pr(>|z|)`, f_star)
levels(ns$Plot4) <- list(C  = "C", `W (5 yrs)` = "y05", `W (27 yrs)` = "y27")
names(ns)[1:2] <- c("Growth form", "Plot treatment")

Plot:

p <- ns |> 
  ggplot(aes(x = `Growth form`, y = rate, color = `Plot treatment`)) + 
  ylab("Number of species per plot") + 
  geom_pointrange(aes(ymin = Lower, ymax = Upper), position = position_dodge(width = 0.3)) +
  geom_text(aes(label = sig, x = `Growth form`, y = Upper + 0.17, group = `Plot treatment`), angle = 90, position = position_dodge(width = 0.34), vjust = 0.5, show.legend = FALSE) +
  theme_bw()
p

Weekly, all species:

df0 <- expand_grid(PlotU = factor(plot_summaries_2020$PlotU), Week = 2:9) |> 
  mutate(Plot4 = factor(if_else(
    PlotU %in% c(1:16, 33:48, 65:80, 97:112), "C", 
    if_else(PlotU %in% c(17:24, 49:56), "y05", 
            if_else(PlotU %in% c(25:32, 57:64, 89:96, 121:128), "y02",
                    if_else(PlotU %in% c(81:88, 113:120), "y27", NA))))),
    Location = if_else(PlotU %in% 1:64, "Imnavait", "Toolik"),
    Site = if_else(PlotU %in% c(1:32, 65:96), "Dry", "Moist"))
df1 <- weekly_summaries_2020 |> 
  right_join(df0) |> 
  mutate(n_sp = replace(n_sp, is.na(n_sp), 0))

f <- function(w){
  m <- glmer(n_sp ~ Plot4 + (1 | Site), family = poisson, data = filter(df1, Week == w))
  emms <- as_tibble(emmeans(m, ~ Plot4, type = "response"))
  s <- f_qp(m) # quasi-poisson adjustment
  return(bind_cols(emms[,1:3], s[,3:4]))
}

# p-values
all <- map_df(2:9, f) |> 
  mutate(Week = rep(2:9, each = 3), .before = Plot4) |> 
  mutate(Group = "All", .before = Plot4)
knitr::kable(all, digits = 4)

Week	Group	Plot4	rate	SE	z value	Pr(>\|z\|)
2	All	C	1.4545	0.3636	1.5903	0.1118
2	All	y05	1.0000	0.4082	-0.8305	0.4062
2	All	y27	1.0000	0.4082	-0.8305	0.4062
3	All	C	2.8182	0.5062	6.6681	0.0000
3	All	y05	3.0000	0.7071	0.2439	0.8073
3	All	y27	4.0000	0.8165	1.4888	0.1365
4	All	C	2.5456	0.4899	4.9768	0.0000
4	All	y05	3.4954	0.7750	1.1244	0.2608
4	All	y27	3.8283	0.8123	1.4843	0.1377
5	All	C	2.3147	0.6698	3.1553	0.0016
5	All	y05	3.3495	1.0248	1.3748	0.1692
5	All	y27	3.0305	0.9515	0.9746	0.3298
6	All	C	1.7023	0.4985	2.0015	0.0453
6	All	y05	2.5820	0.8026	1.3564	0.1750
6	All	y27	1.6138	0.5904	-0.1511	0.8799
7	All	C	0.7273	0.2571	-0.9373	0.3486
7	All	y05	1.8333	0.5528	2.0707	0.0384
7	All	y27	0.6667	0.3333	-0.1479	0.8825
8	All	C	0.5455	0.2227	-1.5853	0.1129
8	All	y05	1.3333	0.4714	1.7671	0.0772
8	All	y27	0.3333	0.2357	-0.6440	0.5196
9	All	C	0.2727	0.1575	-2.1754	0.0296
9	All	y05	0.8333	0.3727	1.4785	0.1393
9	All	y27	0.3333	0.2357	0.2125	0.8317

Weekly, by growth forms:

df0 <- expand_grid(PlotU = factor(plot_summaries_2020$PlotU), Week = 2:9, Group = c("F", "ES", "DS", "GS")) |> 
  mutate(Plot4 = factor(if_else(
    PlotU %in% c(1:16, 33:48, 65:80, 97:112), "C", 
    if_else(PlotU %in% c(17:24, 49:56), "y05", 
            if_else(PlotU %in% c(25:32, 57:64, 89:96, 121:128), "y02",
                    if_else(PlotU %in% c(81:88, 113:120), "y27", NA))))),
    Location = if_else(PlotU %in% 1:64, "Imnavait", "Toolik"),
    Site = if_else(PlotU %in% c(1:32, 65:96), "Dry", "Moist"))
df1 <- weekly_summaries_GF_2020 |> 
  right_join(df0) |> 
  mutate(n_sp = replace(n_sp, is.na(n_sp), 0))

# auxiliary function to compute EMMs, SEs, and p-values
f <- function(w, gf){
  m <- glmer(n_sp ~ Plot4 + (1 | Site), family = poisson, data = filter(df1, Week == w, Group == gf))
  emms <- as_tibble(emmeans(m, ~ Plot4, type = "response"))
  s <- f_qp(m) # quasi-poisson adjustment
  return(bind_cols(emms[,1:3], s[,3:4]))
}

# Forbs
comps <- map_df(2:7, \(x)f(x, "F")) # week 8 not enough obs
wk <- rep(2:7, each = 3)
comps_F <- comps |> 
  mutate(Group = "F", .before = Plot4) |> 
  mutate(Week = wk, .before = Group)
aux <- tibble(Week = rep(8:9, each = 3), Group = "F", Plot4 = rep(c("C", "y05", "y27"), each = 2), rate = numeric(6), SE = numeric(6), `z value` = rep(NA, 6), `Pr(>|z|)` = rep(NA, 6))
comps_F <- bind_rows(comps_F, aux)
comps_F$`Pr(>|z|)`[comps_F$Plot4 == "y05" & comps_F$Week %in% 6:7] <- 1 # method did not converge for weeks 5 and 6

# Evergreen shrubs
comps <- map_df(2:9, \(x)f(x, "ES")) 
wk <- rep(2:9, each = 3)
comps_ES <- comps |> 
  mutate(Group = "ES", .before = Plot4) |> 
  mutate(Week = wk, .before = Group)

# Deciduous shrubs
comps <- map_df(2:9, \(x)f(x, "DS"))
wk <- rep(2:9, each = 3)
comps_DS <- comps |> 
  mutate(Group = "DS", .before = Plot4) |> 
  mutate(Week = wk, .before = Group)
comps_DS$`Pr(>|z|)`[comps_DS$Week %in% 6:7 & comps_DS$Plot4 %in% c("y05", "y27")] <- 1 # method did not converge

# Graminoid
comps <- map_df(2:8, \(x)f(x, "GS"))
wk <- rep(2:8, each = 3)
comps_G <- comps |> 
  mutate(Group = "G", .before = Plot4) |> 
  mutate(Week = wk, .before = Group)
aux <- tibble(Week = 9, Group = "G", Plot4 = c("C", "y05", "y27"), rate = numeric(3), SE = numeric(3), `z value` = rep(NA, 3), `Pr(>|z|)` = rep(NA, 3))
comps_G <- bind_rows(comps_G, aux)

nsp <- all |> 
  bind_rows(comps_F) |> 
  bind_rows(comps_ES) |> 
  bind_rows(comps_DS) |>
  bind_rows(comps_G)
nsp$`z value`[nsp$Plot4 == "C"] <- 0
nsp$`Pr(>|z|)`[nsp$Plot4 == "C"] <- 1
knitr::kable(nsp)

Week	Group	Plot4	rate	SE	z value	Pr(>\|z\|)
2	All	C	1.4545455	0.3636364	0.0000000	1.0000000
2	All	y05	1.0000000	0.4082483	-0.8305299	0.4062392
2	All	y27	1.0000000	0.4082483	-0.8305299	0.4062392
3	All	C	2.8181818	0.5061604	0.0000000	1.0000000
3	All	y05	3.0000000	0.7071068	0.2438714	0.8073304
3	All	y27	4.0000000	0.8164966	1.4888287	0.1365325
4	All	C	2.5456180	0.4898746	0.0000000	1.0000000
4	All	y05	3.4954013	0.7749976	1.1243930	0.2608463
4	All	y27	3.8283142	0.8122781	1.4842553	0.1377412
5	All	C	2.3146672	0.6697912	0.0000000	1.0000000
5	All	y05	3.3495134	1.0247906	1.3747721	0.1692021
5	All	y27	3.0305336	0.9514988	0.9745518	0.3297826
6	All	C	1.7022845	0.4984824	0.0000000	1.0000000
6	All	y05	2.5820101	0.8026105	1.3564336	0.1749612
6	All	y27	1.6137518	0.5903706	-0.1510752	0.8799164
7	All	C	0.7272727	0.2571297	0.0000000	1.0000000
7	All	y05	1.8333333	0.5527708	2.0706562	0.0383909
7	All	y27	0.6666667	0.3333334	-0.1478614	0.8824521
8	All	C	0.5454545	0.2226809	0.0000000	1.0000000
8	All	y05	1.3333333	0.4714045	1.7670875	0.0772136
8	All	y27	0.3333333	0.2357022	-0.6439967	0.5195776
9	All	C	0.2727273	0.1574592	0.0000000	1.0000000
9	All	y05	0.8333333	0.3726780	1.4784656	0.1392832
9	All	y27	0.3333333	0.2357023	0.2124941	0.8317216
2	F	C	0.1818182	0.1285649	0.0000000	1.0000000
2	F	y05	0.0000000	0.0000000	-0.0000019	0.9999985
2	F	y27	0.0000000	0.0000000	-0.0000019	0.9999985
3	F	C	0.3636364	0.1818182	0.0000000	1.0000000
3	F	y05	0.6666667	0.3333333	0.9986143	0.3179816
3	F	y27	0.1666667	0.1666667	-0.8129067	0.4162716
4	F	C	0.2727273	0.1574592	0.0000000	1.0000000
4	F	y05	1.0000000	0.4082483	1.9425448	0.0520712
4	F	y27	0.1666667	0.1666667	-0.4508876	0.6520706
5	F	C	0.3463506	0.3352968	0.0000000	1.0000000
5	F	y05	0.6906307	0.6651269	2.2325365	0.0255795
5	F	y27	0.0000000	0.0000045	-0.0346357	0.9723702
6	F	C	0.2660096	0.0019801	0.0000000	1.0000000
6	F	y05	0.6891326	0.0072543	258.0200743	1.0000000
6	F	y27	0.0000000	0.0000016	-0.0486983	0.9611597
7	F	C	0.0909091	0.0909091	0.0000000	1.0000000
7	F	y05	0.1666667	0.1666667	0.4823764	1.0000000
7	F	y27	0.0000000	0.0000000	-0.0000015	0.9999988
8	F	C	0.0000000	0.0000000	0.0000000	1.0000000
8	F	C	0.0000000	0.0000000	0.0000000	1.0000000
8	F	y05	0.0000000	0.0000000	NA	NA
9	F	y05	0.0000000	0.0000000	NA	NA
9	F	y27	0.0000000	0.0000000	NA	NA
9	F	y27	0.0000000	0.0000000	NA	NA
2	ES	C	0.4545454	0.2032789	0.0000000	1.0000000
2	ES	y05	0.0000000	0.0000000	-0.0000042	0.9999967
2	ES	y27	0.0000000	0.0000000	-0.0000020	0.9999984
3	ES	C	1.4545455	0.3636364	0.0000000	1.0000000
3	ES	y05	1.0000000	0.4082483	-0.7841845	0.4329319
3	ES	y27	2.3333333	0.6236096	1.2938354	0.1957223
4	ES	C	1.6363636	0.3856946	0.0000000	1.0000000
4	ES	y05	1.1666667	0.4409585	-0.8809501	0.3783448
4	ES	y27	2.5000000	0.6454973	1.4060518	0.1597088
5	ES	C	1.2727273	0.3401507	0.0000000	1.0000000
5	ES	y05	1.5000000	0.5000000	0.5211860	0.6022372
5	ES	y27	2.1666667	0.6009252	1.8720350	0.0612018
6	ES	C	1.0000000	0.3015114	0.0000000	1.0000000
6	ES	y05	1.0000000	0.4082483	0.0000000	1.0000000
6	ES	y27	1.1666667	0.4409586	0.4363663	0.6625710
7	ES	C	0.3636364	0.1818182	0.0000000	1.0000000
7	ES	y05	1.0000000	0.4082483	1.8946102	0.0581441
7	ES	y27	0.3333333	0.2357023	-0.1214648	0.9033229
8	ES	C	0.2727273	0.1574591	0.0000000	1.0000000
8	ES	y05	0.8333333	0.3726780	1.6462421	0.0997139
8	ES	y27	0.1666667	0.1666666	-0.4590619	0.6461897
9	ES	C	0.0909091	0.0909090	0.0000000	1.0000000
9	ES	y05	0.3333333	0.2357023	1.0608609	0.2887531
9	ES	y27	0.1666667	0.1666667	0.4286030	0.6682122
2	DS	C	0.6363636	0.2405228	0.0000000	1.0000000
2	DS	y05	0.5000000	0.2886751	-0.4371546	0.6619992
2	DS	y27	1.0000000	0.4082483	1.0162335	0.3095182
3	DS	C	0.6363636	0.2405229	0.0000000	1.0000000
3	DS	y05	0.6666667	0.3333334	0.0903926	0.9279753
3	DS	y27	1.5000000	0.5000000	2.0721896	0.0382478
4	DS	C	0.3504875	0.1960243	0.0000000	1.0000000
4	DS	y05	0.6293552	0.3549984	0.8823806	0.3775710
4	DS	y27	1.1013625	0.5077639	1.9466495	0.0515768
5	DS	C	0.2431016	0.1746863	0.0000000	1.0000000
5	DS	y05	0.2861036	0.2373714	0.1767770	0.8596836
5	DS	y27	0.8583226	0.5195062	1.7669850	0.0772307
6	DS	C	0.2234131	0.0014749	0.0000000	1.0000000
6	DS	y05	0.2605796	0.0024321	23.4365493	1.0000000
6	DS	y27	0.3908477	0.0036483	85.1705876	1.0000000
7	DS	C	0.1701507	0.0017842	0.0000000	1.0000000
7	DS	y05	0.3052209	0.0045994	48.6568671	1.0000000
7	DS	y27	0.3052301	0.0045995	48.6594278	1.0000000
8	DS	C	0.1818182	0.1285648	0.0000000	1.0000000
8	DS	y05	0.1666667	0.1666666	-0.0565388	0.9549126
8	DS	y27	0.1666667	0.1666666	-0.0565388	0.9549126
9	DS	C	0.1818182	0.1285648	0.0000000	1.0000000
9	DS	y05	0.1666667	0.1666666	-0.0565388	0.9549126
9	DS	y27	0.1666667	0.1666666	-0.0565388	0.9549126
2	G	C	0.1814175	0.1288718	0.0000000	1.0000000
2	G	y05	0.4984584	0.2895464	1.3928406	0.1636680
2	G	y27	0.0000000	0.0000001	-0.0000150	0.9999880
3	G	C	0.3636364	0.1818182	0.0000000	1.0000000
3	G	y05	0.6666667	0.3333333	1.2454906	0.2129515
3	G	y27	0.0000000	0.0000000	-0.0000019	0.9999985
4	G	C	0.2727273	0.1574592	0.0000000	1.0000000
4	G	y05	0.6666667	0.3333333	1.6131228	0.1067178
4	G	y27	0.0000000	0.0000000	-0.0000066	0.9999947
5	G	C	0.2727273	0.1574592	0.0000000	1.0000000
5	G	y05	0.5000000	0.2886751	0.9756544	0.3292358
5	G	y27	0.0000000	0.0000000	-0.0000017	0.9999986
6	G	C	0.0909091	0.0909091	0.0000000	1.0000000
6	G	y05	0.3333333	0.2357023	1.2358646	0.2165089
6	G	y27	0.0000000	0.0000000	-0.0000017	0.9999986
7	G	C	0.0909091	0.0909091	0.0000000	1.0000000
7	G	y05	0.3333333	0.2357023	1.2358646	0.2165089
7	G	y27	0.0000000	0.0000000	-0.0000017	0.9999986
8	G	C	0.0909091	0.0909091	0.0000000	1.0000000
8	G	y05	0.3333333	0.2357023	1.2358646	0.2165089
8	G	y27	0.0000000	0.0000000	-0.0000017	0.9999986
9	G	C	0.0000000	0.0000000	0.0000000	1.0000000
9	G	y05	0.0000000	0.0000000	NA	NA
9	G	y27	0.0000000	0.0000000	NA	NA

Weekly plot by growth forms:

df <- nsp |> 
  mutate(Plot4 = factor(Plot4), Group = factor(Group))

df$sig <- map_vec(df$`Pr(>|z|)`, f_star)

# Lower and upper bounds for Mean +- SE
df$Lower <- df$rate - df$SE
df$Upper <- df$rate + df$SE

# Plot groups
df$paired <- rep(NA, nrow(df))
for (i in 1:nrow(df)) {
  df$paired[i] <- paste0(df$Plot4[i], "_", df$Group[i])
}

# Position of stars
df$star_pos <- rep(NA, nrow(df))
df$star_pos[df$Group == "All"] <- 3.5
df$star_pos[df$Group == "F"] <- 1.5
df$star_pos[df$Group == "ES"] <- 1.6
df$star_pos[df$Group == "DS"] <- 2.1
df$star_pos[df$Group == "G"] <- 0.5

levels(df$Plot4) <- list(C  = "C", `W (5 yrs)` = "y05", `W (27 yrs)` = "y27")
levels(df$Group) <- list(All = "All", F  = "F", ES = "ES", DS = "DS", G = "G")
names(df)[2] <- "Growth form"
names(df)[3] <- "Plot treatment"

p <- df %>%
  ggplot(aes(x = Week, y = rate, group = `Plot treatment`, color = `Plot treatment`)) + 
  ylab(as.expression("Number of species per plot")) + 
  geom_point(mapping = aes(x = Week, y = star_pos + 0.25), color = "white", alpha = 0) +
  geom_point(position = position_dodge(width = 0.1)) +
  geom_line(aes(group = paired, color = `Plot treatment`)) +
  geom_errorbar(aes(ymin = Lower, ymax = Upper), width = 0, 
                position = position_dodge(width = 0.1)) +
  geom_text(aes(label = sig, y = star_pos), angle = 90, position = position_dodge(width = 0.15), vjust = 0.6, show.legend = FALSE) +
  facet_grid(rows = vars(`Growth form`), scales = "free_y") +
  theme_bw() 
p

Polen limitation

Data:

fruit_set <- read.csv("fruit_set3.csv")

# Subsets by species
subCT <- filter(fruit_set, Species == "CT")
subDO <- filter(fruit_set, Species == "DO")
subRT <- filter(fruit_set, Species == "RT")
subBO <- filter(fruit_set, Species == "BO")

CT

Sample size by bag treatment X plot treatment:

table(subCT$Bag_Treatment, subCT$Plot_Treatment2, subCT$Location, subCT$Site)

## , ,  = Imnavait,  = Wet
## 
##       
##        1y 4y  C
##   Bag  11  1 11
##   Open  6  0  6
## 
## , ,  = Toolik,  = Wet
## 
##       
##        1y 4y  C
##   Bag   3  0  1
##   Open  2  0  0

Mann-Whitney tests:

subCT.E <- filter(subCT, Plot_Treatment == "W")
w <- wilcox.test(Fruit_set ~ Bag_Treatment, data = subCT.E)
qnorm(w$p.value/2)

## [1] -1.031691

w$p.value

## [1] 0.3022167

w <- wilcox.test(Fruit_set ~ Plot_Treatment, data = subCT)
qnorm(w$p.value/2)

## [1] -2.719638

w$p.value

## [1] 0.006535345

DO

Sample size by bag treatment X plot treatment X location X community type :

table(subDO$Bag_Treatment, subDO$Plot_Treatment2, subDO$Location, subDO$Site)

## , ,  = Toolik,  = Dry
## 
##       
##        26y  C
##   Bag   12  7
##   Open  14  7

Testing for an interaction:

g <- glm(fruits ~ Bag_Treatment * Plot_Treatment, family = "quasipoisson", data = subDO)
s <- summary(g)
s$coefficients[4,]

##   Estimate Std. Error    t value   Pr(>|t|) 
##  0.8343021  0.5799589  1.4385539  0.1589147

Main effects:

g <- glm(fruits ~ Bag_Treatment + Plot_Treatment, family = "quasipoisson", data = subDO)
summary(g)

## 
## Call:
## glm(formula = fruits ~ Bag_Treatment + Plot_Treatment, family = "quasipoisson", 
##     data = subDO)
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         3.1815     0.2505  12.698 4.68e-15 ***
## Bag_TreatmentOpen  -0.2265     0.2232  -1.015   0.3168    
## Plot_TreatmentW     0.5747     0.2631   2.184   0.0354 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 16.05206)
## 
##     Null deviance: 936.25  on 39  degrees of freedom
## Residual deviance: 838.20  on 37  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5

RT

Sample size by bag treatment X plot treatment:

table(subRT$Bag_Treatment, subRT$Plot_Treatment2, subRT$Location, subRT$Site)

## , ,  = Imnavait,  = Wet
## 
##       
##        1y 26y 4y  C
##   Bag   5   0  6  8
##   Open  6   0  2  6
## 
## , ,  = Toolik,  = Wet
## 
##       
##        1y 26y 4y  C
##   Bag   3  13  0  6
##   Open  3   8  0  2

Testing for an interaction:

g <- glmer(fruits ~ offset(log(flowers)) + Bag_Treatment * Plot_Treatment + (1 | Location), family = "poisson", data = subRT)
f_qp(g)[4,]

## # A tibble: 1 × 4
##   Estimate `Std. Error` `z value` `Pr(>|z|)`
##      <dbl>        <dbl>     <dbl>      <dbl>
## 1    0.192        0.324     0.592      0.554

Main effects:

g <- glmer(fruits ~ offset(log(flowers)) + Bag_Treatment + Plot_Treatment + (1 | Location), family = "poisson", data = subRT)
f_qp(g)

## # A tibble: 3 × 4
##   Estimate `Std. Error` `z value`   `Pr(>|z|)`
##      <dbl>        <dbl>     <dbl>        <dbl>
## 1   -1.13         0.205     -5.52 0.0000000335
## 2    0.150        0.118      1.27 0.205       
## 3    0.724        0.157      4.60 0.00000422

BO

Sample size by bag treatment X plot treatment:

table(subBO$Bag_Treatment, subBO$Plot_Treatment2, subBO$Location, subBO$Site)

## , ,  = Imnavait,  = Dry
## 
##       
##        1y 4y  C
##   Bag   2  1 13
##   Open  1  2  9
## 
## , ,  = Imnavait,  = Wet
## 
##       
##        1y 4y  C
##   Bag   1  2  9
##   Open  2  2 10

Testing for an interaction:

g <- glmer(fruits ~ offset(log(flowers)) + Bag_Treatment * Plot_Treatment + (1 | Site), family = "poisson", data = subBO[-45,])
f_qp(g)[4,]

## # A tibble: 1 × 4
##   Estimate `Std. Error` `z value` `Pr(>|z|)`
##      <dbl>        <dbl>     <dbl>      <dbl>
## 1     1.06        0.888      1.19      0.235

Main effects:

g <- glmer(fruits ~ offset(log(flowers)) + Bag_Treatment + Plot_Treatment + (1 | Site), family = "poisson", data = subBO[-45,])
f_qp(g)

## # A tibble: 3 × 4
##   Estimate `Std. Error` `z value`   `Pr(>|z|)`
##      <dbl>        <dbl>     <dbl>        <dbl>
## 1   -1.68         0.295     -5.70 0.0000000119
## 2    0.694        0.257      2.71 0.00683     
## 3   -0.695        0.353     -1.97 0.0491

Effects of experimental warming in the Low Arctic - floral density and number of species per plot.

Phenology

Floral density (2019 only)

Number of species

2019

2020

Polen limitation

CT

DO

RT

BO