Purpose

Analyze selection on and plastic responses of floral and vegetative traits of Ipomopsis aggregata to manipulations of snowmelt and summer precipitation. Determine whether plasticity is adaptive (improves survival and reproduction). Part of a larger project to understand the potential for evolutionary rescue in alpine plants under climate change.

Treatments

library(tidyverse)
library(lubridate)

library(lme4)
library(lmerTest)
library(emmeans)
library(broom)
library(vegan)

library(RColorBrewer)
library(viridis)
library(GGally)

library(knitr)
knitr::opts_chunk$set(comment="", cache=T, warning = F, message = F, fig.path = "images/")

signif.num <- function(x) {
    as.character(symnum(x, corr = FALSE, na = FALSE, legend = FALSE,
           cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1), 
           symbols = c("***", "**", "*", ".", " ")))
}

load("data/maxfield_data.rda")

library(ggnewscale)
subplot_offset <- 1/6
size_meter <- 6
ggplot(treatments_map) + coord_fixed(clip="off")+
  geom_point(aes(x=plot_x, y=plot_y, fill=snow), size=7*size_meter, shape=22, stroke=1)+
  scale_fill_manual("Snowmelt", values=snow_pal, guide=guide_legend(override.aes = list(size = 2*size_meter)))+
  new_scale_fill()+
  geom_point(aes(x=plot_x+subplot_x*subplot_offset, y=plot_y+subplot_y*subplot_offset, fill=water4), size=2*size_meter, shape=22, stroke=1) + 
    geom_text(aes(x=plot_x+subplot_x*subplot_offset, y=plot_y+subplot_y*subplot_offset, label=subplot)) + 
  scale_fill_manual("Precipitation", values=water4_pal) +
  geom_text(aes(x=plot_x, y=plot_y, label=plot))+
  theme(legend.position="left", axis.title=element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), panel.background = element_blank(), panel.grid = element_blank(), legend.box.background = element_blank())

Snow melt vs. treatments

ggplot(hobo, aes(x=intensity, y=temp, color=hour(time))) +  geom_point(alpha=0.5, size=0.5) + geom_density_2d(color="black")+ 
  scale_x_log10(labels=function(x) format(x, scientific=F, big.mark=","), limits=c(10,1e6), breaks=10^(0:6)) +
  scale_color_gradientn(colours = rainbow(5))+
  geom_vline(xintercept = melt_threshold_sun) + 
  geom_hline(yintercept=melt_threshold_warm)+ geom_hline(yintercept=-melt_threshold_warm) +
  labs(x="Light intensity (lux)", y="Temperature (C)", color="Hour") +
  theme_dark() + theme(panel.grid.minor.x=element_blank())

ggplot(hobo %>% filter(hour(time)%in%5:17 & 
                  (year!="2019"&yday(time) %in% 110:130 | 
                   year=="2019"&yday(time) %in% 145:160)), 
       aes(x=intensity, y=temp, color=hour(time))) +   geom_point(alpha=0.7) +geom_density_2d(color="black")+
  scale_x_log10(labels=function(x) format(x, scientific=F, big.mark=","), limits=c(10,1e6), breaks=10^(0:6)) +
  scale_color_gradientn(colours = rainbow(5))+
  geom_vline(xintercept = melt_threshold_sun) + 
  geom_hline(yintercept=melt_threshold_warm)+ geom_hline(yintercept=-melt_threshold_warm) +
  labs(title="Temperature and light thresholds for snow melt", subtitle = "During the day during melt weeks",
       x="Light intensity (lux)", y="Temperature (C)", color="Hour") +
  theme_dark() + theme(panel.grid.minor.x=element_blank())

ggplot(hobo %>% filter(month(time)%in%4:6), aes(x=time, y=temp, color=plot)) + 
  facet_wrap(vars(year), scales="free_x", ncol=1) + 
  geom_hline(yintercept=melt_threshold_warm) + geom_hline(yintercept=-melt_threshold_warm) +geom_line() +
  geom_segment(data=meltdates, aes(x=warm_time, y=-10, xend=sun_time, yend=-15, color=plot), size=1) +
    geom_line(data=daily_temp %>% filter(year!="2017", month(date)%in%4:6), aes(x=as.POSIXct(date), y=av_temp_2m_C), inherit.aes = F, color=brewer.pal(3, "Oranges")[3]) +
  labs(x="Date",y="Temperature (C)", color="Plot")

ggplot(hobo %>% filter(month(time)%in%4:6), aes(x=time, y=intensity, color=plot)) + 
  facet_wrap(vars(year), scales="free_x", ncol=1) + 
  geom_line() + geom_hline(yintercept=melt_threshold_sun) +
  geom_segment(data=meltdates, aes(x=sun_time, y=-2e4, xend=melt_time, yend=-4e4, color=plot), size=1) +
  labs(x="Date",y="Light intensity (lux)", color="Plot")

meltdates %>% select(c("year","plot", "snow", ends_with("date"))) %>% pivot_longer(ends_with("date")) %>% arrange(name) %>% 
  ggplot(aes(y=name, x=value, color=plot, group=plot, linetype=snow)) + facet_wrap(vars(year), ncol=1) +
  geom_path(size=1) + #geom_text(aes(label=plot), hjust=1.5) +
  scale_linetype_manual(values=c(4,1)) + scale_x_continuous(minor_breaks=110:160)+
  labs(x="Day of year",y="Criteron", color="Plot", linetype="Snowmelt")

meltdates %>% ggplot(aes(x=melt_offset, fill=snow)) + geom_dotplot(method="histodot", binwidth=1, color="white") +
  facet_wrap(vars(year), ncol=1)+ scale_fill_manual(values=snow_pal) +
  scale_y_continuous(expand=expansion(0,0))+ scale_x_continuous(breaks=-10:4)+
  theme(panel.grid.minor = element_blank()) +
  labs(x="Offset from mean normal snow melt date (days)", y="Number of plots", fill="Snowmelt")

left_join(sm.subplotmonth, meltdates) %>% 
  mutate(mo = paste0("Month ",mo,", ",month.abb[mo])) %>% 
  ggplot(aes(x=sun_date, y=VWC, shape=year,color=water)) + 
  facet_wrap(vars(mo), nrow=1)+
  geom_point() + geom_smooth(se=F, method="lm") +
  scale_color_manual(values=water_pal) +
  labs(x="Snowmelt date", y="Mean soil moisture in subplot", color="Water", shape="Year")

Weather data

ggplot(daily_precip, aes(y=precip_mm, x=julian)) + geom_col(fill=brewer.pal(3, "Greens")[3]) + facet_wrap(vars(year), ncol=1) + scale_y_continuous(limits=c(0,125)) + labs(x="Day of year", y="Precipiation (mm)") + theme_minimal()

groundcover %>% mutate(first_snow_day=first_snow_day-365) %>% 
  pivot_longer(ends_with("day"), names_to="threshold", values_to="day") %>% 
  mutate(threshold = fct_reorder(factor(threshold),day, na.rm=T, .desc=T)) %>% 
  ggplot(aes(x=year, y=day, color=threshold)) +  
  annotate("rect", xmin=2017.5, xmax=2020.5, ymin=-Inf, ymax=Inf, fill="lightblue")+
  geom_line(aes(group=year), size=2)+ geom_smooth( method="lm", se=F) +geom_point() +
  scale_x_continuous(minor_breaks = 1976:2020, position = "top") + scale_y_continuous(n.breaks=10)+
  scale_color_grey(labels=function(x) str_to_sentence(str_replace_all(x, "_", " ")), start=0, end=0.85) +
  labs(x="Year", y="Day of year", color="Threshold") + theme_minimal()

cat(range(groundcover$year), sep=" - ")

1976 - 2020

cat(paste("SD of first bare ground =", round(sd(groundcover$first_0_cm_day),1),"days"))

SD of first bare ground = 13.5 days

cat(paste("First bare ground trends", round(coef(lm(first_0_cm_day ~ year, data=groundcover))[[2]], 2),"days per year"))

First bare ground trends -0.21 days per year

Soil moisture vs. treatments

ggplot(sm, aes(x=yday(date), shape=snow, linetype=snow, color=water, y = VWC)) + 
  facet_wrap(vars(year), ncol=1) +
  geom_jitter(height=0, width=0.3) + geom_smooth(se=F, span=0.5) +
  geom_col(data=daily_precip %>% filter(year!="2017"), aes(x=julian, y=precip_mm), inherit.aes = F, fill=brewer.pal(3, "Greens")[3]) + 
  geom_line(data=daily_temp %>% filter(year!="2017"), aes(x=julian, y=av_temp_2m_C), inherit.aes = F, color=brewer.pal(3, "Oranges")[3]) + 
  geom_vline(data=meltdates, aes(xintercept=sun_date, linetype=snow), color=brewer.pal(3, "Purples")[3])+
  geom_line(data = snowcloth, aes(x=yday(date), y=mean_snow_depth_cm/10, group=year), inherit.aes=F)+
  geom_point(data = snowcloth, aes(x=yday(date), y=mean_snow_depth_cm/10), inherit.aes=F)+
  geom_line(data=groundcover %>% pivot_longer(ends_with("day"), names_to="threshold", values_to="day") %>% filter(threshold!="first_snow_day", year>2017) %>% mutate(snow_depth_cm = as.integer(str_split_fixed(threshold,"_",n=4)[,2])),
            aes(x=day, y=snow_depth_cm/10), inherit.aes=F, color="grey40")+
  scale_x_continuous(limits=c(85, 230), breaks=scales::breaks_extended(10), expand=expansion(0,0)) + scale_y_continuous(limits=c(-3, 20)) +
  labs(y="Soil VWC (%) / Daily precipitation (mm) / Mean temp. (C) / Snowpack (dm)", 
       x="Day of year", color="Precipitation", shape="Snow", linetype="Snow") +
  scale_color_manual(values=water_pal) + scale_linetype_manual(values=c(2,1)) + theme_minimal() + 
  theme(legend.key.size = unit(2, "lines"))

Census

The annual census at the beginning of the summer contains information on mortality, flowering status, and plant size (longest leaf and number of leaves for each rosette).

Load census data

Tally the number of each status combination for each year.

nrf_count <- bind_rows(
  allowed_nrf %>% filter(year==2019) %>% left_join(cen19 %>% mutate(nrf=paste(notes, rosettes, flowering)) %>% group_by(nrf) %>% tally()),
  allowed_nrf %>% filter(year==2020) %>% left_join(cen20 %>% mutate(nrf=paste(notes, rosettes, flowering)) %>% group_by(nrf) %>% tally()))

kable(nrf_count %>% select(-nrf) %>% arrange(year, status, desc(n)))

year	notes	rosettes	flowering	status	n
2019	nf	blank	blank	dead_nf	215
2019	dead	zero	n	dead_nf	143
2019	blank	zero	n	dead_nf	18
2019	nf	blank	n	dead_nf	7
2019	blank	blank	n	dead_nf	6
2019	dead	blank	n	dead_nf	3
2019	dead	one_or_more	n	dead_nf	3
2019	dead	zero	blank	dead_nf	1
2019	dead	zero	blank	dead_nf	1
2019	nf	zero	n	dead_nf	1
2019	blank	blank	y	flowering	110
2019	blank	one_or_more	y	flowering	51
2019	blank	indist	y	flowering	3
2019	nf	one_or_more	blank	recheck	3
2019	blank	zero	y	recheck	2
2019	blank	blank	blank	recheck	1
2019	chewed	zero	y	recheck	1
2019	nf	blank	y	recheck	1
2019	tagnf	one_or_more	blank	recheck	1
2019	tagnf	blank	blank	tagnf	46
2019	tagnf	blank	n	tagnf	3
2019	tagnf	zero	blank	tagnf	1
2019	blank	one_or_more	n	vegetative	433
2019	blank	indist	n	vegetative	75
2019	chewed	zero	n	vegetative	23
2019	blank	one_or_more	blank	vegetative	19
2019	chewed	one_or_more	n	vegetative	7
2019	chewed	indist	n	vegetative	2
2020	nf	zero	blank	dead_nf	490
2020	dead	zero	blank	dead_nf	120
2020	blank	one_or_more	y	flowering	151
2020	blank	blank	y	flowering	50
2020	chewed	one_or_more	y	flowering	17
2020	chewed	indist	y	flowering	1
2020	chewed	blank	blank	recheck	3
2020	chewed	zero	blank	recheck	1
2020	chewed	blank	y	recheck	1
2020	chewed	zero	y	recheck	1
2020	tagnf	blank	blank	tagnf	76
2020	blank	one_or_more	n	vegetative	233
2020	chewed	one_or_more	n	vegetative	6
2020	chewed	zero	n	vegetative	2

Duplicates

When there is more than one record for a plant tag in the census, which duplicate should be kept? Currently resolved by the following rules:

sort by status (keep flowering > vegetative > recheck > dead_nf > tagnf)
sort by date within status (if statuses are the same, keep the observation with the newest date)
keep the first row and discard the other duplicates

This is automatic, but it may be better to resolve duplicates manually in some cases. The duplicates are written out to spreadsheets

#write out the duplicates
#cen19 %>% group_by(id) %>% filter(n()>1) %>% arrange(id, status, desc(date)) %>% 
#  write_sheet(ss=filter(datasheets, name=="Maxfield Results"), sheet="2019duplicates")
#cen20 %>% group_by(id) %>% filter(n()>1) %>% arrange(id, status, desc(date)) %>% 
#  write_sheet(ss=filter(datasheets, name=="Maxfield Results"), sheet="2020duplicates")

#distinct() keeps the first row after sorting, discards the other rows
cen19 <- cen19 %>% arrange(id, status, desc(date)) %>% distinct(id, .keep_all = TRUE)
cen20 <- cen20 %>% arrange(id, status, desc(date)) %>% distinct(id, .keep_all = TRUE)

Status vs. treatment

ggplot(drop_na(cen19, water, snow) %>% filter(status %in% status_ok), aes(x=water, fill=status)) + facet_wrap(vars(snow), labeller=labeller(.cols=label_both)) + geom_bar(position="fill") + scale_fill_manual(values=status_pal) + labs(title="2019 status", x="Precipitation", y="Proportion", fill="Status")

ggplot(drop_na(cen20, water, snow) %>% filter(status %in% status_ok), aes(x=water, fill=status)) + facet_wrap(vars(snow), labeller=labeller(.cols=label_both)) + geom_bar(position="fill") + scale_fill_manual(values=status_pal) + labs(title="2020 status", x="Precipitation", y="Proportion", fill="Status")

Join 2019 and 2020 datasets

Check these tags missing in each year to make sure tags were not mis-read. When tags were not found in the other year, the part of the row for that year is “NA”.

Recorded in 2019 but not in 2020:

sort(setdiff(cen19[cen19$notes != "tagnf",]$id, cen20$id))

 [1] "-642"   "-643"   "1B-244" "1B-245" "1B-246" "1B-247" "1B-248" "1B-252"
 [9] "1B-662" "1B-679" "1B-680" "1B-681" "1B-682" "1B-683" "1D-31"  "1D-32" 
[17] "3A-659" "3B-602" "3B-603" "3B-604" "3B-605" "3B-611" "3B-612" "3B-615"
[25] "3B-617" "3B-618" "3B-622" "3B-623" "3D-852" "4A-100" "4A-66"  "4A-83" 
[33] "4A-90"  "4B-162" "4B-172" "4B-192" "4C-115" "4C-29"  "4C-32"  "4D-39" 
[41] "4D-61"  "5A-172" "5A-177" "5B-14"  "5B-287" "5C-205" "5C-213" "5D-231"
[49] "5D-236" "5D-284" "6B-32"  "6C-103" "6D-80"  "6D-81"

Recorded in 2020 but not in 2019:

sort(setdiff(cen20[cen20$notes != "tagnf",]$id, cen19$id))

 [1] "1A-244" "1B-199" "1B-66"  "2A-100" "2A-589" "2A-684" "2A-689" "2A-99" 
 [9] "2B-171" "2C-256" "2D-225" "3C-120" "3C-185" "3C-213" "3C-86"  "3C-94" 
[17] "3D-255" "4A-124" "4A-168" "4B-594" "4B-596" "4B-680" "4C-641" "4C-642"
[25] "4C-643" "4D-66"  "5B-30"  "5B-32"  "5B-88"  "5C-"    "5C-683" "6A-171"
[33] "6A-685" "6B-120" "6B-243" "6B-266" "6C-271" "6C-447" "6D-18"  "6D-86"

Status transitions

What did plants of a given status in 2019 look like now in 2020?

with(cen, table(`2019`=paste(status_19), `2020`=paste(status)))

            2020
2019         dead_nf flowering  NA recheck tagnf vegetative
  dead_nf        294        21  27       1    18         51
  flowering      127        15   9       1     7         16
  NA              21        15   0       0     1          4
  recheck          5         0   0       0     3          1
  tagnf           13         7  15       0    12          5
  vegetative     164       167  24       4    38        175

ggplot(filter(cen, status_19 %in% status_ok, status %in% status_ok), aes(x=status_19, fill=status)) + geom_bar(position="fill") + scale_fill_manual(values=status_pal) + labs(fill="2020 Status", y="Proportion", x="2019 Status")

Survival and flowering vs. plant size

ggplot(filter(cen, status %in% status_ok), aes(x=r1_longest_19, fill=status)) + facet_wrap(vars(status), ncol=1) +
  geom_density() + scale_fill_manual(values=status_pal) + labs(x="2019 Rosette 1 Longest Leaf (mm)", fill="2020 status")

ggplot(filter(cen, status %in% status_ok), aes(x=r1_leaves_19, fill=status)) + facet_wrap(vars(status), ncol=1) +
  geom_density() + scale_fill_manual(values=status_pal) + labs(x="2019 Rosette 1 # Leaves", fill="2020 status")

Plant growth vs. treatments

Only includes plants that remained vegetative and did not die or flower.

ggplot(drop_na(cen, water, snow), aes(x=r1_longest_19, y=r1_longest, color=water, shape=snow)) + geom_point() + geom_abline(intercept=0, slope=1) + scale_color_manual(values=water_pal) + theme_minimal()+ labs(x="2019 Rosette 1 Longest Leaf (mm)", y="2020 Rosette 1 Longest Leaf (mm)",  color="Precipitation", shape="Snowmelt")

ggplot(drop_na(cen, water, snow), aes(x=r1_leaves_19, y=r1_leaves, color=water, shape=snow)) + geom_point() + geom_abline(intercept=0, slope=1 )+ scale_color_manual(values=water_pal) + theme_minimal()+ labs(x="2019 Rosette 1 # Leaves", y="2020 Rosette 1 # Leaves",  color="Precipitation", shape="Snowmelt")

Plant size vs. treatments

ggplot(drop_na(cen19, water, snow), aes(x=water, color=snow, y=r1_longest)) + geom_boxplot() + scale_color_manual(values=snow_pal) + labs(x="Precipitation", y="2019 Rosette 1 Longest Leaf (mm)",  color="Snowmelt")

ggplot(drop_na(cen19, water, snow), aes(x=water, color=snow, y=r1_leaves)) + geom_boxplot() + scale_color_manual(values=snow_pal)+ labs(x="Precipitation", y="2019 Rosette 1 # Leaves",  color="Snowmelt")

ggplot(drop_na(cen20, water, snow), aes(x=water, color=snow, y=r1_longest)) + geom_boxplot() + scale_color_manual(values=snow_pal) + labs(x="Precipitation", y="2020 Rosette 1 Longest Leaf (mm)",  color="Snowmelt")

ggplot(drop_na(cen20, water, snow), aes(x=water, color=snow, y=r1_leaves)) + geom_boxplot() + scale_color_manual(values=snow_pal) + labs(x="Precipitation", y="2020 Rosette 1 # Leaves",  color="Snowmelt")

Leaf traits

ggplot(lt, aes(x=wet_weight_g, y=dry_weight_g)) + geom_point() + geom_abline(slope=1, intercept = 0) + geom_abline(slope=1-median(lt$water_content), intercept = 0)

ggplot(lt, aes(x=leaf_area_cm2, y=dry_weight_g)) + geom_point() + geom_smooth(method="lm")

ggplot(lt, aes(x=leaf_area_cm2, y=wet_weight_g)) + geom_point() + geom_smooth(method="lm")

ggplot(lt, aes(x=water, color=snow, y = sla)) + geom_boxplot() + scale_color_manual(values=snow_pal) + labs(x="Precipitation", y="Specific Leaf Area (cm2/g)", color="Snowmelt") + facet_wrap(vars(year, round), nrow=1) + scale_y_continuous(limits=c(0,300))

ggplot(lt, aes(x=sla, y=water_content, color=water, shape=snow, linetype=snow)) + geom_point() +  scale_color_manual(values=water_pal) + theme_minimal() + 
  theme(legend.key.size = unit(2, "lines"))

ggplot(lt, aes(x=water, color=snow, y = water_content)) + geom_boxplot() + scale_color_manual(values=snow_pal) + labs(x="Precipitation", y="Water content", color="Snowmelt") + facet_wrap(vars(year, round), nrow=1)

ggplot(lt, aes(x=water, color=snow, y = trichome_density)) + geom_boxplot() + scale_color_manual(values=snow_pal) + labs(x="Precipitation", y="Trichome density (/cm2)", color="Snowmelt") + facet_wrap(vars(year, round), nrow=1)

Floral traits

Floral morphology and nectar

Floral trait correlations

traits <- c("corolla_length", "corolla_width", "style_length", "sepal_width", "nectar_24_h_ul", "nectar_conc","nectar_sugar_24_h_mg") #exclude anthers
traitnames <- setNames(c("Corolla Length (mm)", "Corolla Width (mm)", "Style Length (mm)", "Sepal Width (mm)", "Nectar Volume (uL/day)", "Nectar Concentration (% sucrose by mass)", "Nectar Sucrose (mg/day)"), traits)

ggpairs(mt %>% select(corolla_length,   corolla_width,  style_length,   anther_max, anther_min, sepal_width))

ggplot(mt, aes(x=nectar_24_h_ul, y = nectar_conc)) + geom_point() + labs(x="Nectar volume (uL / day)", y="Nectar concentration (% sucrose by mass)") + geom_smooth(se=F)

Floral plasticity: water treatments and snowmelt treatments

This is a split-plot analysis, based on the original experimental design. The models use:

mt.subplot, the 24 * 3 subplot means within each year
snow, the snowmelt treatment (except for plot 2 in 2019)
water, the water treatments with mock rainouts and controls combined
plot as a random effect (subplots vary within this larger unit)

The models are specified as:

lmer(trait ~ year * snow * water + (1|plot), data = mt.subplot)

The models for nectar volume and total sugar are singular (variance at plot level near 0)
Removing the interactions with year does not make those two models not singular
The model fits are mostly singular when a random effect of plotid (subplot) is added for repeated measures among years

options(contrasts = c("contr.sum", "contr.poly"))
mt.mod <- mt.mod.test <- mod.emm <- mt.mod.plot <- list()
for(tr in traits) {
  mt.mod[[tr]] <- lmer(mt.subplot[[tr]] ~ year * snow * water + (1|plot), data=mt.subplot)
  mt.mod.test[[tr]] <- anova(mt.mod[[tr]])
  mod.emm[[tr]] <- summary(emmeans(ref_grid(mt.mod[[tr]]), ~ year + snow * water))
  
  mt.mod.plot[[tr]] <- 
    ggplot(mod.emm[[tr]], aes(x=water, color=snow, y = emmean)) + facet_wrap(vars(year), nrow=1) +
        labs(x="Precipitation", y=traitnames[[tr]], color="Snowmelt", shape="Year") +
        geom_point(data=mt.subplot, aes_string(x="water", color="snow", y = tr), inherit.aes = F, position=position_dodge(0.7)) + 
        geom_point(position=position_dodge(0.7), size=4, shape=0) + 
        geom_linerange(aes(ymin=emmean-SE, ymax=emmean+SE), position=position_dodge2(0.6)) +
        guides(x=guide_axis(angle = 90)) + scale_color_manual(values=snow_pal)
  print(mt.mod.plot[[tr]])
  cat(paste(tr,"\n"))
  print(mt.mod.test[[tr]])
}

corolla_length 
Type III Analysis of Variance Table with Satterthwaite's method
                Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
year            79.898  39.949     2 50.123 22.2232 1.231e-07 ***
snow             0.313   0.313     1  6.733  0.1738   0.68969    
water            9.865   4.933     2 49.183  2.7439   0.07417 .  
year:snow        0.218   0.109     2 51.925  0.0605   0.94136    
year:water       4.777   1.194     4 49.183  0.6644   0.61972    
snow:water       0.518   0.259     2 49.183  0.1440   0.86629    
year:snow:water  3.254   0.813     4 49.183  0.4525   0.77005    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

corolla_width 
Type III Analysis of Variance Table with Satterthwaite's method
                 Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
year            2.99723 1.49862     2 49.546 56.9010 1.462e-13 ***
snow            0.00342 0.00342     1  7.538  0.1300 0.7283296    
water           0.50804 0.25402     2 47.896  9.6448 0.0003021 ***
year:snow       0.01566 0.00783     2 50.143  0.2974 0.7440619    
year:water      0.14078 0.03520     4 47.896  1.3363 0.2701610    
snow:water      0.02453 0.01226     2 47.896  0.4656 0.6305408    
year:snow:water 0.19365 0.04841     4 47.896  1.8382 0.1369561    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

style_length 
Type III Analysis of Variance Table with Satterthwaite's method
                Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
year            69.663  34.832     2 50.026 20.5058 3.133e-07 ***
snow             0.153   0.153     1  6.477  0.0898    0.7738    
water            7.542   3.771     2 49.084  2.2200    0.1194    
year:snow        0.049   0.024     2 51.928  0.0143    0.9858    
year:water       5.425   1.356     4 49.084  0.7985    0.5320    
snow:water       1.564   0.782     2 49.084  0.4604    0.6338    
year:snow:water 10.740   2.685     4 49.084  1.5807    0.1943    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sepal_width 
Type III Analysis of Variance Table with Satterthwaite's method
                  Sum Sq  Mean Sq NumDF  DenDF F value   Pr(>F)   
year            0.000199 0.000199     1 34.077  0.0241 0.877513   
snow            0.000888 0.000888     1 11.500  0.1074 0.748998   
water           0.128507 0.064253     2 31.783  7.7706 0.001789 **
year:snow       0.016125 0.016125     1 34.183  1.9501 0.171583   
year:water      0.095841 0.047920     2 31.783  5.7953 0.007148 **
snow:water      0.003098 0.001549     2 31.783  0.1873 0.830090   
year:snow:water 0.034221 0.017110     2 31.783  2.0693 0.142956   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

nectar_24_h_ul 
Type III Analysis of Variance Table with Satterthwaite's method
                Sum Sq Mean Sq NumDF DenDF F value   Pr(>F)   
year            4.8137 2.40683     2    54  6.7666 0.002386 **
snow            1.5321 1.53213     1    54  4.3075 0.042722 * 
water           2.3693 1.18467     2    54  3.3306 0.043255 * 
year:snow       1.1263 0.56316     2    54  1.5833 0.214682   
year:water      0.4961 0.12403     4    54  0.3487 0.843800   
snow:water      0.2006 0.10030     2    54  0.2820 0.755398   
year:snow:water 1.1010 0.27524     4    54  0.7738 0.547040   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

nectar_conc 
Type III Analysis of Variance Table with Satterthwaite's method
                Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
year            626.29 313.143     2 51.355 32.9352 6.255e-10 ***
snow              8.03   8.028     1 12.023  0.8444  0.376216    
water           123.73  61.863     2 49.819  6.5065  0.003087 ** 
year:snow        11.87   5.934     2 51.246  0.6241  0.539760    
year:water       14.38   3.594     4 49.819  0.3780  0.823275    
snow:water       29.12  14.562     2 49.819  1.5315  0.226211    
year:snow:water  62.30  15.574     4 49.819  1.6380  0.179412    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

nectar_sugar_24_h_mg 
Type III Analysis of Variance Table with Satterthwaite's method
                Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
year            4839.6 2419.78     2    54  8.0412 0.0008774 ***
snow             499.7  499.67     1    54  1.6605 0.2030319    
water            349.0  174.49     2    54  0.5799 0.5634262    
year:snow        918.3  459.13     2    54  1.5257 0.2266851    
year:water      1041.9  260.47     4    54  0.8656 0.4905963    
snow:water       307.7  153.83     2    54  0.5112 0.6026421    
year:snow:water 1223.4  305.86     4    54  1.0164 0.4071820    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

mt.mod.tests <- map_dfr(mt.mod.test, tidy, .id="trait")
mt.mod.tests %>% select(c("trait", "term", "p.value")) %>% mutate(p.value = signif.num(p.value)) %>% 
  pivot_wider(names_from=trait, values_from=p.value) %>% 
  bind_rows(lapply(mt.mod, function(x) ifelse(isSingular(x),"Y","N"))) %>% mutate(term = replace_na(term,"Singular?")) %>% kable

term	corolla_length	corolla_width	style_length	sepal_width	nectar_24_h_ul	nectar_conc	nectar_sugar_24_h_mg
year	***	***	***		**	***	***
snow					*
water	.	***		**	*	**
year:snow
year:water				**
snow:water
year:snow:water
Singular?	N	N	N	N	Y	N	Y

Floral plasticity: water treatments and snowmelt offset

The differences in snowmelt dates vary among years and plots, so this analysis uses a melt date offset. The models use:

mt.subplot, the 24 * 3 subplot means within each year
water, the water treatments with mock rainouts and controls combined
melt_offset, the number of days between the snowmelt date (when the light entering the HOBO logger exceeds 10,000 lux) and the mean snowmelt date for the normal plots in that year (according to snow)
plotid (subplot) as a random effect (plants vary within this larger unit)

The models are specified as:

lmer(trait ~ year * water * melt_offset + (1|plotid), data = mt.subplot)

options(contrasts = c("contr.sum", "contr.poly"))
mt.mod.raw <- mt.mod.plantyr <- mt.mod.subplot <- mt.mod.subplot.test <- mt.mod.plot <- list()
for(tr in traits) {
  #mt.mod.raw[[tr]]     <- lmer(mt[[tr]]         ~ year * water * melt_offset + (1|plotid) + (1|plantid), data=mt)
  #mt.mod.plantyr[[tr]] <- lmer(mt.plantyr[[tr]] ~ year * water * melt_offset + (1|plotid), data=mt.plantyr)
  mt.mod.subplot[[tr]] <- lmer(mt.subplot[[tr]] ~ year * water * melt_offset + (1|plot), data=mt.subplot)
  mt.subplot[[paste0("fit.",tr)]] <- NA
  mt.subplot[[paste0("fit.",tr)]][!is.na(mt.subplot[[tr]])]<- predict(mt.mod.subplot[[tr]], re.form=NA)
  mt.mod.subplot.test[[tr]] <- anova(mt.mod.subplot[[tr]])
  
  mt.mod.plot[[tr]] <- 
    ggplot(mt.subplot, aes_string(x="melt_offset", color="water", y = tr)) + geom_point() +
    geom_line(data=mt.subplot %>% drop_na(paste0("fit.",tr)), aes_string(y=paste0("fit.",tr)), size=2)+
     scale_color_manual(values=water_pal)+ facet_wrap(vars(year))+
        labs(x="Offset from mean normal snow melt date (days)", y=traitnames[[tr]], color="Precipitation")
  print(mt.mod.plot[[tr]])
}

mt.mod.subplot.tests <- map_dfr(mt.mod.subplot.test, tidy, .id="trait")
mt.mod.subplot.tests %>% select(c("trait", "term", "p.value")) %>% mutate(p.value = signif.num(p.value)) %>% 
  pivot_wider(names_from=trait, values_from=p.value) %>% 
  bind_rows(lapply(mt.mod.subplot, function(x) ifelse(isSingular(x),"Y","N"))) %>% mutate(term = replace_na(term,"Singular?")) %>% kable

term	corolla_length	corolla_width	style_length	sepal_width	nectar_24_h_ul	nectar_conc	nectar_sugar_24_h_mg
year	***	***	***		**	***	*
water		***		*		**
melt_offset					*
year:water		.		.
year:melt_offset
water:melt_offset
year:water:melt_offset
Singular?	N	N	Y	N	Y	N	Y

Floral plasticity: soil moisture

Both snowmelt timing and summer precipitation affect soil moisture and its timing. This analysis regresses the floral traits against the mean soil moisture in a subplot. Slopes may differ among years. The models use:

mt.plantyr, the trait data averaged by plant within year
VWC, the summer mean of the approximately weekly mean volumetric water content in a subplot. This may need to change to a subset of the summer because it includes soil moisture taken after traits were measured. See section on within-season variation for timing of measurements in each year.
plotid (subplot) as a random effect (plants vary within this larger unit)

The models are specified as:

lmer(trait ~ year * VWC + (1|plotid), data = mt.plantyr)

The model for total sugar is singular (variance at subplot level near 0)

options(contrasts = c("contr.sum", "contr.poly"))
mt.sm.mod <- mt.sm.mod.test <- mt.sm.mod.raw <- coef.VWC <- emm.VWC <- mt.sm.mod.plot <-list()
for(tr in traits) {
  #mt.sm.mod.raw[[tr]] <- lmer(mt[[tr]]    ~ year * VWC + (1|plot) + (1|plantid), data=mt)
  mt.sm.mod[[tr]] <-lmer(mt.plantyr[[tr]] ~ year * VWC + (1|plotid), data=mt.plantyr)
  mt.plantyr[[paste0("fit.sm.",tr)]] <- NA
  mt.plantyr[[paste0("fit.sm.",tr)]][!is.na(mt.plantyr[[tr]])]<- predict(mt.sm.mod[[tr]], re.form=NA)
  mt.sm.mod.test[[tr]] <- anova(mt.sm.mod[[tr]])
  
  #Z-score the data first to get standardized slopes
  mt.sm.mod[[tr]] <- lmer(scale(mt.plantyr[[tr]]) ~ year * VWC + (1|plotid), data=mt.plantyr)
  coef.VWC[[tr]] <- summary(mt.sm.mod[[tr]])$coefficients["VWC",]
  emm.VWC[[tr]] <- emtrends(mt.sm.mod[[tr]],  ~ year, var = "VWC")

   mt.sm.mod.plot[[tr]] <- 
    ggplot(mt.plantyr, aes_string(x="VWC", color="year", y = tr)) + geom_point() +
    geom_line(data=mt.plantyr %>% drop_na(paste0("fit.sm.",tr)), aes_string(y=paste0("fit.sm.",tr)), size=2)+
     scale_color_manual(values=year_pal)+
        labs(x="Summer mean of soil moisture in subplot (VWC)", y=traitnames[[tr]], color="Year")
  print(mt.sm.mod.plot[[tr]])
}

mt.sm.mod.tests <- map_dfr(mt.sm.mod.test, tidy, .id="trait")
mt.sm.mod.tests %>% select(c("trait", "term", "p.value")) %>% mutate(p.value = signif.num(p.value)) %>% 
  pivot_wider(names_from=trait, values_from=p.value) %>% 
    bind_rows(lapply(mt.sm.mod, function(x) ifelse(isSingular(x),"Y","N"))) %>% mutate(term = replace_na(term,"Singular?")) %>% kable

term	corolla_length	corolla_width	style_length	sepal_width	nectar_24_h_ul	nectar_conc	nectar_sugar_24_h_mg
year	**	**		**		***	**
VWC	*	***	**	*	***	*	***
year:VWC		*		**			*
Singular?	N	N	N	N	N	N	Y

coef.VWC.df <- as.data.frame(t(as.data.frame(coef.VWC))) %>% rownames_to_column("trait")
emm.VWC.ls <- lapply(emm.VWC, as.data.frame)
emm.VWC.ls$sepal_width[3,] <- NA
emm.VWC.df  <- bind_rows(emm.VWC.ls, .id="trait")

ggplot(emm.VWC.df %>% drop_na(), aes(x=trait, y=VWC.trend, color=year)) +  
  geom_hline(yintercept=0) + geom_pointrange(aes(ymin=VWC.trend - SE, ymax = VWC.trend + SE), position=position_dodge(width=0.3), size=0.4) + 
  labs(y="Slope of z-scored trait vs. mean annual soil moisture", x="Floral Trait", color="Year") +
  theme_classic() + scale_x_discrete(labels = traitnames) + coord_flip() + scale_color_manual(values=year_pal)

Multivariate response

mt.morph <- mt[,c("corolla_length", "corolla_width","style_length","anther_max","anther_min","year","VWC")] %>% drop_na(1:5) #drop sepal_width
mt.pca <- prcomp(mt.morph[,-(6:7)], scale.=TRUE)

mt.morph.pca <- mt.morph %>%
  bind_cols(as.data.frame(mt.pca$x)) %>% 
  tibble::rowid_to_column("index") %>% 
  mutate(x_0=0,y_0=0,
         x_1=-corolla_width/2,y_1=corolla_length,
         x_2= corolla_width/2,y_2=corolla_length) %>% 
  pivot_longer(c(starts_with("x"),starts_with("y")), names_to=c(".value", "point"), names_sep="_")

scl <- 100
zoom <- 5.5
mt.pca.vars <- as.data.frame(mt.pca$rotation) %>% rownames_to_column("variable")
mt.pca.labs <- paste0("PC", 1:2, " (", round(100*mt.pca$sdev[1:2]^2/sum(mt.pca$sdev^2), 0), "% explained var.)") 
ggplot(mt.morph.pca) + theme_classic() +
  geom_polygon(aes(x= x/scl+PC1, y= y/scl+PC2, group=index, color=VWC, fill=VWC)) +
  #geom_point(aes(x=PC1, y=PC2, size=sepal_width/2), color="darkgreen") + 
  geom_point(aes(x=PC1, y=PC2 + anther_min/scl), color="yellow", size=0.3) +
  geom_point(aes(x=PC1, y=PC2 + anther_max/scl), color="yellow", size=0.3) +
  geom_point(aes(x=PC1, y=PC2 + style_length/scl), color="grey60", size=0.5, shape=4) +
  geom_segment(data=mt.pca.vars,aes(x=0, y=0, xend=PC1*zoom, yend=PC2*zoom), 
               arrow = arrow(length = unit(1/2, "picas")), color="black") +
  geom_text(data =mt.pca.vars, aes(label = gsub("_"," ",variable), x = PC1*zoom*1.4, y = PC2*zoom*1.4,
                                   angle =  (180/pi) * atan(PC2/PC1), hjust = 0), color="black")+
  scale_size_identity() + coord_fixed() + 
  scale_fill_viridis(direction=-1, option="plasma") + scale_color_viridis(direction=-1, option="plasma") + 
  labs(x=mt.pca.labs[1], y=mt.pca.labs[2])

ggsave("images/mt_morph_pca_noSW.pdf", height=15, width=20, device=cairo_pdf)

ggplot(mt, aes(x=corolla_width, y=corolla_length, color=VWC)) + scale_color_viridis(direction=-1, option="plasma") + geom_point() + theme_minimal()

ggplot(mt, aes(y=style_length / corolla_length, x=VWC, color=VWC)) + scale_color_viridis(direction=-1, option="plasma") + geom_point() + theme_minimal() + geom_smooth()

options(contrasts = c("contr.treatment", "contr.poly"))
mt.cap <- capscale(scale(mt.morph[,-(6:7)]) ~ VWC + year, data=mt.morph)
mt.cap

Call: capscale(formula = scale(mt.morph[, -(6:7)]) ~ VWC + year, data =
mt.morph)

              Inertia Proportion Rank
Total          5.0000     1.0000     
Constrained    1.0367     0.2073    3
Unconstrained  3.9633     0.7927    5
Inertia is mean squared Euclidean distance 
Species scores projected from 'scale' 'mt.morph[, -(6:7)]' 

Eigenvalues for constrained axes:
  CAP1   CAP2   CAP3 
0.9909 0.0447 0.0011 

Eigenvalues for unconstrained axes:
  MDS1   MDS2   MDS3   MDS4   MDS5 
2.6291 0.6967 0.4033 0.1385 0.0957

anova(mt.cap, by="term") %>% tidy %>% kable

term	df	Variance	statistic	p.value
VWC	1	0.4492375	107.79564	0.001
year	2	0.5874772	70.48329	0.001
Residual	951	3.9632853	NA	NA

library(ggvegan)
autoplot(mt.cap) + guides(color=F, shape=F)

mt.morph.wm <- mt[,c("corolla_length", "corolla_width","style_length","anther_max","anther_min","year","water","melt_offset")] %>% drop_na(1:6) #and nectar (19,20) too
mt.cap.wm <- capscale(scale(mt.morph.wm[,-(6:8)]) ~ year * melt_offset * water, data=mt.morph.wm)
mt.cap.wm

Call: capscale(formula = scale(mt.morph.wm[, -(6:8)]) ~ year *
melt_offset * water, data = mt.morph.wm)

              Inertia Proportion Rank
Total          5.0000     1.0000     
Constrained    1.1459     0.2292    5
Unconstrained  3.8541     0.7708    5
Inertia is mean squared Euclidean distance 
Species scores projected from 'scale' 'mt.morph.wm[, -(6:8)]' 

Eigenvalues for constrained axes:
  CAP1   CAP2   CAP3   CAP4   CAP5 
1.0587 0.0517 0.0284 0.0048 0.0023 

Eigenvalues for unconstrained axes:
  MDS1   MDS2   MDS3   MDS4   MDS5 
2.5579 0.6711 0.3985 0.1337 0.0930

anova(mt.cap.wm, by="term") %>% tidy %>% kable

term	df	Variance	statistic	p.value
year	2	0.9178296	111.5715011	0.001
melt_offset	1	0.0344690	8.3801187	0.003
water	2	0.1267323	15.4055996	0.001
year:melt_offset	2	0.0080217	0.9751170	0.397
year:water	4	0.0289412	1.7590459	0.083
melt_offset:water	2	0.0062569	0.7605915	0.540
year:melt_offset:water	4	0.0236900	1.4398801	0.167
Residual	937	3.8540593	NA	NA

#plot CAP without interactions
autoplot(capscale(scale(mt.morph.wm[,-(6:8)]) ~ year + melt_offset + water, data=mt.morph.wm)) + guides(color=F, shape=F)

Floral traits within season

No obvious effect of time within a season for corolla length/width or sepal width. The timing of trait measurements was different each year:

ggplot(sm, aes(x=yday(date))) + 
  facet_wrap(vars(year), ncol=1) +
  #geom_vline(xintercept=175) + geom_vline(xintercept=183) +
  geom_point(data=mt, aes(y= corolla_length), color="purple")+
    geom_point(data=mt, aes(y= nectar_24_h_ul), color="orange")+
  geom_point(aes(y = VWC), color="black") +
  geom_vline(data=meltdates, aes(xintercept=sun_date, linetype=snow), color="grey40")+
  labs(title="Timing of measurements", subtitle="VWC (black) / Morphology (purple) / Nectar (orange)", x="Day of year", shape="Snow", linetype="Snow")+
  scale_linetype_manual(values=c(2,1)) + theme_minimal() + 
  theme(legend.key.size = unit(2, "lines"))

ggplot(data=mt, aes(x=yday(date), y= corolla_length, color=water)) + 
  facet_wrap(vars(year), ncol=1) +
  geom_point()+   geom_smooth(se=F)+ 
  labs(y="Corolla Length (mm)", x="Day of year")+
  scale_color_manual(values=water_pal) + theme_minimal()

ggplot(data=mt, aes(x=yday(date), y= nectar_24_h_ul, color=water)) + 
  facet_wrap(vars(year), ncol=1) +
  geom_point()+   geom_smooth(se=F)+ 
  labs(y="Nectar amount (48 hr, mm)", x="Day of year")+
  scale_color_manual(values=water_pal) + theme_minimal()

ggplot(data=nt, aes(x=yday(date), y= nectar_conc, color=water)) + 
  facet_wrap(vars(year), ncol=1) +
  geom_point()+   geom_smooth(se=F)+ 
  labs(y="Nectar concentration", x="Day of year")+
  scale_color_manual(values=water_pal) + theme_minimal()

Phenology

ggplot(ph20, aes(x=julian, y=height_cm, shape=snow, color=water))+ geom_point() + geom_path(aes(group=plantid))+
  labs(y="Plant height", 
       x="Day of year", color="Precipitation", shape="Snow") +
  scale_color_manual(values=water_pal) + theme_minimal() + 
  theme(legend.key.size = unit(2, "lines"))

#ggplot(ph20, aes(x=julian, y=open+buds, shape=snow, linetype=snow, color=water))+ geom_point() + geom_path(aes(group=plantid)) +
#  labs(y="Total flowers and buds", 
#       x="Day of year", color="Precipitation", shape="Snow", linetype="Snow") +
#  scale_color_manual(values=water_pal) + scale_linetype_manual(values=c(2,1)) + theme_minimal() + 
#  theme(legend.key.size = unit(2, "lines"))

# ggplot(ph20, aes(x=julian, y=eggs, shape=snow, linetype=snow, color=water))+ geom_point() + geom_path(aes(group=plantid)) +
#   labs(y="Total eggs", 
#        x="Day of year", color="Precipitation", shape="Snow", linetype="Snow") +
#   scale_color_manual(values=water_pal) + scale_linetype_manual(values=c(2,1)) + theme_minimal() + 
#   theme(legend.key.size = unit(2, "lines"))

ph <-bind_rows(list("2020"= ph20,
                    "2018"= ph18 %>% mutate(plant=as.character(plant))), .id="year_df") %>% mutate(julian=as.integer(as.character(julian)))

ph %>% group_by(year_df, snow, water, julian) %>% summarize_at(c("open","buds"), mean) %>% 
  ggplot(aes(x=julian, y=open, shape=snow, linetype=snow, color=water))+ geom_point() + 
  facet_wrap(vars(year_df), nrow=2)+
  geom_path(aes(group=paste(snow, water))) +
  labs(y="Mean total flowers", 
       x="Day of year", color="Precipitation", shape="Snow", linetype="Snow") +
  scale_color_manual(values=water_pal) + scale_linetype_manual(values=c(2,1)) + theme_minimal() + 
  theme(legend.key.size = unit(2, "lines"))

ph %>% group_by(year_df, snow, water, julian) %>% summarize_at("flowering", mean, na.rm=T) %>% 
  ggplot(aes(x=julian, y=flowering, shape=snow, linetype=snow, color=water))+ geom_point() + 
    facet_wrap(vars(year_df), nrow=2)+
  geom_path(aes(group=paste(snow, water))) +
  labs(y="Proportion flowering", 
       x="Day of year", color="Precipitation", shape="Snow", linetype="Snow") +
  scale_y_continuous(labels = scales::percent)+
  scale_color_manual(values=water_pal) + scale_linetype_manual(values=c(2,1)) + theme_minimal() + 
  theme(legend.key.size = unit(2, "lines"))

ph %>% group_by(year_df, snow, water, julian) %>% summarize_at("eggs", mean, na.rm=T) %>% 
  ggplot(aes(x=julian, y=eggs, shape=snow, linetype=snow, color=water))+ geom_point() + 
    facet_wrap(vars(year_df), nrow=2)+
  geom_path(aes(group=paste(snow, water))) +
  labs(y="Mean total eggs", 
       x="Day of year", color="Precipitation", shape="Snow", linetype="Snow") +
  scale_color_manual(values=water_pal) + scale_linetype_manual(values=c(2,1)) + theme_minimal() + 
  theme(legend.key.size = unit(2, "lines"))

ph %>% group_by(year_df, snow, water, julian) %>% summarize_at("has_egg", mean, na.rm=T) %>% 
  ggplot(aes(x=julian, y=has_egg, shape=snow, linetype=snow, color=water))+ geom_point() + 
    facet_wrap(vars(year_df), nrow=2)+
  geom_path(aes(group=paste(snow, water))) +
  labs(y="Proportion of plants with eggs", 
       x="Day of year", color="Precipitation", shape="Snow", linetype="Snow") +
  scale_y_continuous(labels = scales::percent)+
  scale_color_manual(values=water_pal) + scale_linetype_manual(values=c(2,1)) + theme_minimal() + 
  theme(legend.key.size = unit(2, "lines"))

Floral display size

2020

ph20.plantid <- ph20 %>% left_join(sm.subplotyear) %>% group_by(plotid, plantid, snow, water) %>% summarize_at(c("open","height_cm", "VWC"), mean, na.rm=T)

ggplot(ph20.plantid, aes(x=height_cm, y=open, color=water)) + geom_point() + geom_smooth(method="lm", se=F) + scale_color_manual(values=water_pal) +
  labs(x="Mean inflorescence height (cm)", y="Mean number of open flowers", color="Precipitation")

ggplot(ph20.plantid, aes(x=water, color=snow, y = height_cm)) +
    labs(x="Precipitation", y="Mean inflorescence height (cm)", color="Snowmelt") +
    geom_boxplot() + scale_color_manual(values=snow_pal)

ggplot(ph20.plantid, aes(x=VWC, y=height_cm)) + geom_point() + geom_smooth() +
    labs(x="Mean VWC", y="Mean inflorescence height (cm)")

ggplot(ph20.plantid, aes(x=water, color=snow, y = open)) +
    labs(x="Precipitation", y="Mean number of open flowers", color="Snowmelt") +
    geom_boxplot() + scale_color_manual(values=snow_pal)

ggplot(ph20.plantid, aes(x=VWC, y=open)) + geom_point() + geom_smooth()+
    labs(x="Mean VWC", y="Mean number of open flowers")

Seeds

2018

ggplot(sds18, aes(color=water, fruits, seeds_total)) + geom_point() + geom_smooth(se=F) + scale_color_manual(values=water_pal)

ggplot(sds18, aes(color=water, water, seeds_total)) + geom_boxplot() + scale_color_manual(values=water_pal)
ggplot(sds18, aes(color=water, water, fruits)) + geom_boxplot() + scale_color_manual(values=water_pal)
ggplot(sds18, aes(color=water, water, seeds_total/fruits)) + geom_boxplot() + scale_color_manual(values=water_pal)
ggplot(sds18, aes(color=water, water, fruits_fly_no_seeds)) + geom_boxplot() + scale_color_manual(values=water_pal)
ggplot(sds18, aes(color=water, water, fruits_fly_no_seeds/fruits)) + geom_boxplot() + scale_color_manual(values=water_pal)

Maxfield Meadow Ipomopsis Experiment

John Powers

2020-11-11

Purpose

Treatments

Snow melt vs. treatments

Weather data

Soil moisture vs. treatments

Census

Load census data

Duplicates

Status vs. treatment

Join 2019 and 2020 datasets

Status transitions

Survival and flowering vs. plant size

Plant growth vs. treatments

Plant size vs. treatments

Leaf traits

Floral traits

Floral morphology and nectar

Floral trait correlations

Floral plasticity: water treatments and snowmelt treatments

Floral plasticity: water treatments and snowmelt offset

Floral plasticity: soil moisture

Multivariate response

Floral traits within season

Phenology

Floral display size

2020

Seeds

2018