---
title: "bioassay_CutFlower_consumption"
author: "Sonja Glasser"
date: '2024-11-20'
output: html_document
---

I will be combining the cutflower consumption data with the evaporation controls for each pollen type and dry time. 
I will use this data to generate curves to calculate the amount of pollen consumed. 

Sometimes the pollen was dried for two days and sometimes the pollen was dried for three days. 
The pollen dried for 2 days were either dried for 45 or 48 hours. For these pollen diets I will use the evaporation controls dried for each pollen type for 48 hours. 
The pollen dried for three days were either dried for 64, 70, 71, or 72 hours, for this pollen I will use the evaporation controls dried for 68 hours pollen for all pollen types. 

```{r setup, include=FALSE}
require(here)#telling R where to look for a file
require(lme4) # generalized mixed models 
library(performance) # check model performance
library(tidyverse) # language & ggplot
library(car) # Anova
library(emmeans) #extract the estimated marginal means from the model
library(effects)#allEffects()
library(MuMIn)#dredge () model selection
library(AICcmodavg) #compare models 
library(GGally) # pairwise plots 
library(lattice) # another graphing package
library("glmmTMB")
library("bbmle")
library(multcomp)
library(broom) #tidy() function
```


```{r}
counts <- read.csv(here("data", "counts_rdataset.csv"))
consume <- read.csv(here("data", "consumption_bioassay_submit.csv"))
evap_curve <- read.csv(here("data", "evaporation_controls_submit.csv"))
```


```{r}
consume<-consume%>%
  rename(Bee_ID = "Bee.ID..")
cols<- c("Bee_ID","type","tot_time_dry", "evap_cntrl_time")
consume[cols] <- lapply(consume[cols], factor)
consume<-consume[,c(1,2,5,14,16,17,19,20)]

summary(consume)

cols2<- c("ID","type","evap_cntrl_time")
evap_curve[cols2] <- lapply(evap_curve[cols2], factor)
evap_curve<-evap_curve[,1:13]
evap_curve<-evap_curve%>% 
  filter(evap_cntrl_time!="NA")
summary(evap_curve)
```


First let's visualize the evaporation controls 

```{r}
library(viridis)
evap_controls_48 <- evap_curve %>%
  filter(evap_cntrl_time == "48")

ggplot(aes(x=pol_wet, y=pol_dry, col=type), data=evap_controls_48) + geom_point() + labs(title = "Regression for control pollen 48 hours dried") + theme_bw()


evap_controls_68 <- evap_curve %>%
  filter(evap_cntrl_time == "68")

ggplot(aes(x=pol_wet, y=pol_dry, col=type), data=evap_controls_68) + geom_point() + labs(title = "Regression for control pollen 68 hours dried") + theme_bw()

#difference between 48 and 68 hours... I will just apply the 68 hours dry time to the 48 hour pollen type A treatments. 
evap_B <- evap_curve %>%
  filter(type == "B")
ggplot(aes(x=pol_wet, y=pol_dry, col=evap_cntrl_time), data=evap_B) + geom_point() + labs(title = "Regression for pollen B 48 & 68 hours dried") + theme_bw()


```

First, I will do this for just one pollen type for dried for 2 days (Treatment = = "B"). 

```{r}
# evap = weights of the pollen without bees (evaporation controls)
# first filter out the diet we want for the control weights
evap_B_48 <- evap_curve %>%
  filter(type == "B", evap_cntrl_time == "48") 
# make a linear regression with the old & new control weights
evap_B_lg_48 <- lm(pol_dry ~ pol_wet + 0, data =evap_B_48)
summary(evap_B_lg_48)
# plot to see the relationship
plot(pol_dry ~ pol_wet, data =evap_B_48)
abline(evap_B_lg_48)
##Now with the intercept is NOT constrained to zero because of poor fit with the data
##this likely indicates a non-linear regression near the origin (e.g. very fast evaporation at very small pollen ball initial weights)
evap_B_lg_48 <- lm(pol_dry ~ pol_wet, data =evap_B_48)
summary(evap_B_lg_48)
# plot to see the relationship
plot(pol_dry ~ pol_wet, data = evap_B_48)
abline(evap_B_lg_48)
```


With the regression we can predict, given an initial weight, what the final "should" be. The amount that the bees consumed is the predicted final weight (lost to evaporation) subtracted by the measured final weight (evaporation + consumption).

Here I will try just with one species of pollen (B) at 48 hours dry time. 
```{r}
# cons = weights of the pollen from consumption data
# now take the values from the pollen fed to the bees
cons_48_B <- consume %>%
  filter(evap_cntrl_time == "48" & type == "B") 
# now predict what those values would be based on our evaporation model:
cons_48_B $predicted_final_weight <- predict(evap_B_lg_48, cons_48_B )

cons_48_B <- cons_48_B %>%
  mutate(pollen_consumed = predicted_final_weight - pol_dry)
plot(cons_48_B$pol_wet,cons_48_B$pollen_consumed) # I am not sure what this graph is telling me... but I have the
hist(cons_48_B$pollen_consumed)

```

Here is a function to calculate across pollen types. 

```{r}
evap_lg <- evap_curve %>%
  group_by(type, evap_cntrl_time) %>%
  do(pollen_lm = tidy(lm(pol_dry ~ pol_wet, data = .))) %>% 
  unnest(pollen_lm) %>%
  group_by(type) %>%
  pivot_wider(names_from=term, values_from=c("estimate":"p.value"))

evap_cons_merged<- merge(consume, evap_lg, by=c("type","evap_cntrl_time"))  %>%
  rename(intercept_estimate = "estimate_(Intercept)") %>% 
  mutate(predicted_final_pollen = estimate_pol_wet * pol_wet   + intercept_estimate) %>%
  mutate(pollen_consumed = predicted_final_pollen - pol_dry)

#test to make sure this worked. 
sun_48<-consume%>%
  filter(type=="SUN", evap_cntrl_time=="48")
sun_48_merge<-evap_cons_merged%>%
  filter(type=="SUN", evap_cntrl_time=="48")

testing_SUN <- merge(sun_48, sun_48_merge, by="Bee_ID")
  
plot(testing_SUN$pol_consumed.x,testing_SUN$pol_consumed.y) 
#looks good

histogram(evap_cons_merged$pollen_consumed)

```

Now I need to merge consumption data and the count data. 

```{r}
evap_cons_merged_filt <- evap_cons_merged %>%
  dplyr::select("Bee_ID", "pollen_consumed")

count_cons<- merge(evap_cons_merged_filt, counts, by="Bee_ID")
```


I will now make a new r dataset and use that for the data analysis. 

```{r}
write.csv(count_cons,file=here::here("data","count_cons_rdataset.csv"), row.names=FALSE)
```