Back to Article

###### 2-2 cage - locust physiology anaylses

Download Notebook

This script analyzes the locust data (protein and carbohydrates) for the field cage experiments

Please see manuscript for specific methodology.

This script does the following:

1. analyze locust weight, surivial, and # adult in field cages
2. analyze locust nutrient rebalancing through time

Manuscript figure preparation happens in the `03_field_cages_intake_target_mauscript_figure_code.ipynb` file

    In [2]:

    load_packages <- function(packages) {
      # Check for uninstalled packages
      uninstalled <- packages[!packages %in% installed.packages()[,"Package"]]

      # Install uninstalled packages
      if(length(uninstalled)) install.packages(uninstalled, dependencies = TRUE)

      # Load all packages
      for (pkg in packages) {
        if (!require(pkg, character.only = TRUE, quietly = TRUE)) {
          message(paste("Failed to load package:", pkg))
        }
      }
    }

    # List of packages to check, install, and load
    packages <- c("mgcv", "MuMIn", "multcomp", "gratia",
                  "ggpubr", "patchwork", "broom", "knitr", "janitor",
                  "here","ggpubr","MetBrewer","GGally","tidyverse")

    load_packages(packages)

    i_am('README.md')

    # Functions

    ## Standard error
    std <- function(x) sd(x)/sqrt(length(x))

    # Setting R options for jupyterlab
    options(repr.plot.width = 10, repr.plot.height = 10,repr.matrix.max.rows=10)

    here() starts at /home/datascience/herbivore_nutrient_interactions

# Section 1 - analyze locust weight, surivial, and # adult in field cages

## read in data

    In [3]:

    adult_survival_prop <- read_csv(here("data/raw/field_cage_experiment/LocustCageData_v02.xlsx - by cage.csv"),show_col_types = FALSE)

    final_mass <- read_csv(here("data/raw/field_cage_experiment/LocustCageData.csv"),show_col_types = FALSE) |>
      mutate(Treatment = fct_recode(Treatment, "None" = "0N", "Medium" = "MedH", "High" = "HighN"))

    New names:
    • `` -> `...24`
    • `` -> `...25`
    • `` -> `...26`
    • `` -> `...27`
    • `` -> `...28`

## Now create two dataframes

- survival and adult proporation at end of field cage experiment
- final mass of all surviving locusts at end of experiment

    In [6]:

    adult_survival_prop2 <- adult_survival_prop |>
        dplyr::select(`N level`,cageNum,AdultProportion,SurviveProportion) |>
        clean_names() |>
        mutate(n_level = factor(n_level),
               treatment = fct_recode(n_level, "None" = '0', "Medium" = '87.5', "High" = '175')) |>
        select(!n_level)

    In [7]:

    final_mass2 <- final_mass |>
        clean_names() |>
        filter(!complete.cases(oedaleus_australis)) |>
        select(treatment,mass_g,sex)

## save to disk for figure construction

    In [8]:

    write.csv(adult_survival_prop2, here("data/raw/field_cage_experiment/adult_locust_proportion_survival_by_treatment.csv"),row.names=FALSE)

    write.csv(final_mass2, here("data/raw/field_cage_experiment/adult_locust_final_mass_by_treatment.csv"),row.names=FALSE)

# Section 2 - analyze locust nutrient rebalancing through time

## read in data

We need to read in the data and conduct some significant management

    In [9]:

    # Read and clean the names
    dat <- suppressMessages(read_csv(here("data/raw/field_cage_experiment/Field_Cage_IT.csv"),
                    show_col_types = FALSE, skip = 1) |>
      clean_names() |>
      as_tibble())

    dat |> head()

A tibble: 6 × 53
locust_id_num | field_plot | field_cage_num | date_collected_from_field | color | sex | locust_initialmass_g | locusts_final_mass_g | notes | molt_date | ⋯ | high_prot_final_g_44 | x45 | diet_change5_start_time | diet_change5_end_time | high_carb_dish_num_48 | high_carb_initial_g_49 | high_carb_final_g_50 | high_prot_dish_num_51 | high_prot_initial_g_52 | high_prot_final_g_53
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
<dbl> | <chr> | <dbl> | <chr> | <chr> | <chr> | <dbl> | <dbl> | <chr> | <chr> | ⋯ | <dbl> | <lgl> | <chr> | <chr> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl>
101 | C1 | 5 | Nov 23, 2015 7:00 PM | B | F | 0.353 | 0.330 | NA | Nov 30, 2015 | ⋯ | 2.888 | NA | Nov 29, 2015 8:00 PM | Dec 2, 2015 | 123 | 2.954 | 2.941 | 90 | 2.924 | 2.9250
102 | C1 | 5 | Nov 23, 2015 7:00 PM | B | M | 0.167 | 0.165 | NA | Dec 1, 2015 | ⋯ | 2.887 | NA | NA | NA | 146 | 2.868 | 2.868 | 85 | 2.856 | 2.8540
103 | C1 | 6 | Nov 23, 2015 7:00 PM | B | F | 0.300 | NA | Dead2nd day; prob no eat | NA | ⋯ | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA
104 | C1 | 7 | Nov 23, 2015 7:00 PM | G | F | 0.262 | NA | Died 1 Dec | NA | ⋯ | 2.932 | NA | NA | NA | 106 | 2.841 | 2.839 | 101 | 2.852 | 2.8532
105 | H1 | 1 | Nov 23, 2015 7:00 PM | B | F | 0.156 | NA | Dead2nd day; prob no eat | NA | ⋯ | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA
106 | H1 | 2 | Nov 23, 2015 7:00 PM | B | F | 0.211 | 0.328 | NA | NA | ⋯ | 2.905 | NA | NA | NA | 116 | 2.894 | 2.858 | 58 | 2.906 | 2.8680

## data management step 1

- dropping NAs
- identifying treatments based on cage number
- creating time series for locust IT (step 1-5)

    In [11]:

    dat2 <- dat |>
      drop_na(locusts_final_mass_g) |>
      mutate(
        sex = if_else(sex == 'B', 'M', sex),
        treatment = case_when(
          field_cage_num %in% c(1:4, 21:24, 29:32) ~ "High",
          field_cage_num %in% c(9:12, 17:20, 25:28) ~ "Med",
          field_cage_num %in% c(5:8, 13:16, 33:36) ~ "none"),
        total_c_step1 = (`high_carb_initial_g_14` - `high_carb_final_g_15`) * 0.42,
        total_p_step1 = (`high_prot_initial_g_17` - `high_prot_final_g_18`) * 0.42,
        total_c_step2 = (`high_carb_initial_g_22` - `high_carb_final_g_23`) * 0.42,
        total_p_step2 = (`high_prot_initial_g_25` - `high_prot_final_g_26`) * 0.42,
        total_c_step3 = (`high_carb_initial_g_31` - `high_carb_final_g_32`) * 0.42,
        total_p_step3 = (`high_prot_initial_g_34` - `high_prot_final_g_35`) * 0.42,
        total_c_step4 = (`high_carb_initial_g_40` - `high_carb_final_g_41`) * 0.42,
        total_p_step4 = (`high_prot_initial_g_43` - `high_prot_final_g_44`) * 0.42,
        total_c_step5 = (`high_carb_initial_g_49` - `high_carb_final_g_50`) * 0.42,
        total_p_step5 = (`high_prot_initial_g_52` - `high_prot_final_g_53`) * 0.42,
      )

    In [12]:

    dat2

A tibble: 30 × 64
locust_id_num | field_plot | field_cage_num | date_collected_from_field | color | sex | locust_initialmass_g | locusts_final_mass_g | notes | molt_date | ⋯ | total_c_step1 | total_p_step1 | total_c_step2 | total_p_step2 | total_c_step3 | total_p_step3 | total_c_step4 | total_p_step4 | total_c_step5 | total_p_step5
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
<dbl> | <chr> | <dbl> | <chr> | <chr> | <chr> | <dbl> | <dbl> | <chr> | <chr> | ⋯ | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl>
101 | C1 | 5 | Nov 23, 2015 7:00 PM | B | F | 0.353 | 0.330 | NA | Nov 30, 2015 | ⋯ | 0.00000 | 0.00000 | 0.01848 | 0.00168 | 0.00840 | 0.00252 | 0.00042 | 0.00126 | 0.00546 | -0.00042
102 | C1 | 5 | Nov 23, 2015 7:00 PM | B | M | 0.167 | 0.165 | NA | Dec 1, 2015 | ⋯ | 0.00000 | 0.00042 | 0.00756 | 0.00000 | 0.00966 | 0.00252 | 0.00000 | 0.00420 | 0.00000 | 0.00084
106 | H1 | 2 | Nov 23, 2015 7:00 PM | B | F | 0.211 | 0.328 | NA | NA | ⋯ | 0.00798 | 0.00000 | 0.01050 | 0.00000 | 0.03738 | 0.00882 | 0.02478 | 0.01512 | 0.01512 | 0.01596
107 | H1 | 2 | Nov 23, 2015 7:00 PM | G | F | 0.216 | 0.313 | NA | NA | ⋯ | 0.01764 | 0.00042 | 0.01512 | 0.00168 | 0.01974 | 0.01218 | 0.02226 | 0.01680 | 0.00084 | 0.00000
108 | H1 | 3 | Nov 23, 2015 7:00 PM | B | F | 0.192 | 0.252 | NA | NA | ⋯ | 0.01344 | 0.00000 | 0.01008 | 0.00168 | 0.01176 | 0.00336 | 0.02016 | 0.00336 | 0.01764 | 0.00168
⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋱ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮
133 | C3 | 34 | Nov 23, 2015 7:00 PM | M | M | 0.121 | 0.138 | NA | NA | ⋯ | 0.00252 | 0.00000 | 0.00420 | 0 | 0.00714 | 0.00168 | 0.00672 | 0.00126 | 0.00336 | 0.00420
134 | C3 | 34 | Nov 23, 2015 7:00 PM | B | F | 0.293 | 0.213 | Wing buds separating | Nov 30, 2015 | ⋯ | 0.00546 | 0.00000 | 0.00336 | 0 | 0.00252 | 0.00042 | 0.00126 | 0.00084 | -0.83664 | 0.00084
136 | C3 | 36 | Nov 23, 2015 7:00 PM | B | F | 0.231 | 0.225 | NA | NA | ⋯ | -0.00126 | 0.00000 | 0.00000 | 0 | 0.01764 | 0.00168 | 0.00588 | 0.00126 | 0.00042 | 0.00042
137 | H3 | 31 | Nov 23, 2015 7:00 PM | G | F | 0.215 | 0.336 | NA | NA | ⋯ | 0.00294 | 0.00042 | 0.00588 | 0 | 0.00714 | 0.00210 | 0.01302 | 0.01050 | 0.02142 | 0.01680
138 | H3 | 32 | Nov 23, 2015 7:00 PM | B | M | 0.134 | 0.135 | NA | NA | ⋯ | 0.00126 | 0.00000 | 0.00462 | 0 | 0.00630 | 0.00168 | 0.00462 | 0.00252 | 0.00924 | 0.00336

## data management step 2

- filter consumption step for negative numbers
- creating a row ID
- collasping the time steps into a longer format
- recoding the time step into a simple number

    In [13]:

    dat3 <- dat2|>
      filter(
        total_c_step1 >= 0, total_p_step1 >= 0,
        total_c_step2 >= 0, total_p_step2 >= 0,
        total_c_step3 >= 0, total_p_step3 >= 0,
        total_c_step4 >= 0, total_p_step4 >= 0,
        total_c_step5 >= 0, total_p_step5 >= 0
      ) |>
      mutate(id = row_number()) |>
      pivot_longer(cols = c(starts_with("total_c_"), starts_with("total_p_")),
                   names_to = c(".value", "interval"),
                   names_pattern = "(.*)_(.*)") |>
      select(field_plot,field_cage_num, sex,locust_initialmass_g,locusts_final_mass_g, treatment,id,interval,total_c,total_p) |>
      mutate(interval = recode(interval,
                               step1 = "1",
                               step2 = "2",
                               step3 = "3",
                               step4 = "4",
                               step5 = "5"),
            id = factor(id),
            treatment = factor(treatment),
            interval = factor(interval))  |>
        group_by(id) |>
        arrange(id,interval) %>%
      mutate(
        cumulative_total_c = cumsum(total_c),
        cumulative_total_p = cumsum(total_p),
        day_interval = case_when(interval == 1 ~'Day 1',
                                 interval == 2 ~'Day 2',
                                 interval == 3 ~'Day 3-4',
                                 interval == 4 ~'Day 5-6',
                                 interval == 5 ~'Day 7-9')
      ) |>
      ungroup()

## raw data visualization

    In [14]:

    treatment_raw_plot <- dat3 |>
        ggplot(aes(x=total_p,y=total_c,color=treatment)) +
            geom_point() +
            coord_equal(ratio=1) +
            xlab("Protein consumed (g)") +
            ylab("Carbohydrates consumed (g)") +
            theme_pubr() +
            geom_abline(slope = 1,linetype=1) +
            geom_abline(slope = 2,linetype=2) +
            facet_wrap(~day_interval) +
            scale_color_manual(values=c("#4daf4a","#984ea3")) +
            xlim(0,0.04) +
            ylim(0,0.04)

    sex_raw_plot <- dat3 |>
        ggplot(aes(x=total_p,y=total_c,color=sex)) +
            geom_point() +
            coord_equal(ratio=1) +
            xlab("Protein consumed (g)") +
            ylab("Carbohydrates consumed (g)") +
            theme_pubr() +
            geom_abline(slope = 1,linetype=1) +
            geom_abline(slope = 2,linetype=2) +
            facet_wrap(~day_interval) +
            scale_color_manual(values=c("#dd5129","#0f7ba2")) +
            xlim(0,0.04) +
            ylim(0,0.04)

### treatment raw plot

    In [15]:

    treatment_raw_plot

### sex raw plot

    In [16]:

    sex_raw_plot

### save plots to disk

    In [17]:

    # Define the output directory
    output_dir <- here('output/field_cage_experiment/raw_visualization/')

    # List of plots and their corresponding file names
    plots <- list(
      "intake_target_treatment_raw_graph" = treatment_raw_plot,
      "intake_target_sex_raw_graph" = sex_raw_plot
    )

    # Loop through each plot and save it
    for (name in names(plots)) {
      ggsave(
        filename = file.path(output_dir, paste0(name, ".png")),
        plot = plots[[name]],
        width = 15, height = 10, dpi = 300
      )
    }

## Modeling

### Model Comparison Overview: Locust Nutrient Consumption

This analysis focuses on the nutrient cimulative consumption of carbohydrates (cumulative_total_c) and proteins (cumulative_total_p) by locusts over time (interval), with key predictors such as sex, treatment, and other covariates. The following six Generalized Additive Models (GAMs) are used to assess the impact of these factors:

#### Model 1: Full Model with Random Effects for Repeated Measures

    mod1 <- gam(list(
      cumulative_total_c ~ sex + interval + treatment + s(id, bs="re"),
      cumulative_total_p ~ sex + interval + treatment + s(id, bs="re")),
      family=mvn(d=2), select=TRUE, data=dat3
    )

- Fixed Effects: Includes `sex`, `interval`, and `treatment`.
- Random Effects:
  - `s(id, bs="re")`: Random effect to account for repeated measures on individual locusts (`id`), but note the sample size may not be sufficient for this structure.
- Purpose: This model captures the effects of sex, time, and treatment while accounting for individual locust variability.

#### Model 2: Interaction Between Interval and Treatment, No Smooth Terms

    mod2 <- gam(list(
      cumulative_total_c ~ sex + (interval) * treatment,
      cumulative_total_p ~ sex + (interval) * treatment),
      family=mvn(d=2), select=TRUE, data=dat3
    )

- Fixed Effects: Includes `sex`, `interval`, `treatment`, and their interaction.
- Purpose: Simplifies the model by removing the smooth term for initial mass, focusing solely on the interaction between time and treatment.

#### Model 3: Fixed Effects for Sex Only

    mod3 <- gam(list(
      cumulative_total_c ~ sex,
      cumulative_otal_p ~ sex),
      family=mvn(d=2), select=TRUE, data=dat3
    )

- Fixed Effects: Only includes `sex`.
- Purpose: This model isolates the effect of sex on nutrient consumption without considering time, treatment, or mass.

#### Model 4: Fixed Effects for Treatment Only

    mod4 <- gam(list(
      cumulative_total_c ~ treatment,
      cumulative_total_p ~ treatment),
      family=mvn(d=2), select=TRUE, data=dat3
    )

- Fixed Effects: Only includes `treatment`.
- Purpose: This model tests the direct effect of treatment on nutrient consumption without accounting for sex, time, or mass.

#### Model 5: Intercept-Only Model (Null Model)

    mod5 <- gam(list(
      cumulative_total_c ~ 1,
      cumulative_total_p ~ 1),
      family=mvn(d=2), select=TRUE, data=dat3
    )

- Fixed Effects: Only includes an intercept term.
- Purpose: Serves as a null model, estimating the overall mean of carbohydrate and protein consumption without considering any predictors.

### Summary

By comparing these models, we can assess the influence of sex, treatment, interval, and initial mass on locust nutrient consumption. The models range from full random effects models (Model 1) to simplified intercept-only models (Model 6), allowing for a thorough exploration of the factors driving nutrient consumption.

    In [18]:

    mod1 <- gam(list(
      cumulative_total_c ~  sex  +  interval + treatment + s(id, bs="re"),
      cumulative_total_p ~  sex + interval + treatment  + s(id, bs="re")),
      family=mvn(d=2),select=TRUE, data=dat3
    )

    mod2 <- gam(list(
      cumulative_total_c ~  sex + interval * treatment,
      cumulative_total_p ~  sex + interval * treatment),
      family=mvn(d=2),select=TRUE, data=dat3
    )

    mod3 <- gam(list(
      cumulative_total_c ~  sex ,
      cumulative_total_p ~  sex),
      family=mvn(d=2),select=TRUE, data=dat3
    )

    mod4 <- gam(list(
      cumulative_total_c ~   treatment,
      cumulative_total_p ~   treatment),
      family=mvn(d=2),select=TRUE, data=dat3
    )

    mod5 <- gam(list(
      cumulative_total_c ~   1,
      cumulative_total_p ~   1),
      family=mvn(d=2),select=TRUE, data=dat3
    )

### Model Diagnostics

    In [19]:

    gam.check(mod1)

    Method: REML   Optimizer: outer newton
    full convergence after 6 iterations.
    Gradient range [-0.000203556,0.0005737438]
    (score -954.8838 & scale 1).
    Hessian positive definite, eigenvalue range [1.80807,7.847553].
    Model rank =  63 / 63

    Basis dimension (k) checking results. Low p-value (k-index<1) may
    indicate that k is too low, especially if edf is close to k'.

              k'  edf k-index p-value
    s(id)   23.0 17.6      NA      NA
    s.1(id) 23.0 10.6      NA      NA

These diagnostic plots indicate these models could be improved. I think for descriptive uses these models are okay.

Assuming a gaussian distribution is probably limiting here.

### Model selection criteria

    In [20]:

    # Model selection using AIC, BIC, and AICc
    bic <- BIC(mod1, mod2, mod3, mod4, mod5) |>  # Bayesian Information Criterion
        mutate(deltaBIC = BIC - min(BIC)) |>
        rownames_to_column('model') |>
        dplyr::select(model,BIC,deltaBIC)

    aic <- AIC(mod1, mod2, mod3, mod4, mod5)  |>  # Akaike Information Criterion
        mutate(deltaAIC = AIC - min(AIC)) |>
        rownames_to_column('model') |>
        dplyr::select(model,AIC,deltaAIC)

    aicc <- AICc(mod1, mod2, mod3, mod4, mod5)  |>  # AIC adjusted for small sample size
        mutate(deltaAICc = AICc - min(AICc)) |>
        rownames_to_column('model') |>
        dplyr::select(model,AICc,deltaAICc)

    model_selection_tables <- bic |>
        left_join(aic,by='model') |>
        left_join(aicc,by='model') |>
        mutate(experiment = 'locust_intake_target_through_time')

    In [21]:

    model_selection_tables |> arrange(deltaAIC)

A data.frame: 5 × 8
model | BIC | deltaBIC | AIC | deltaAIC | AICc | deltaAICc | experiment
--- | --- | --- | --- | --- | --- | --- | ---
<chr> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <chr>
mod1 | -1959.594 | 0.00000 | -2086.993 | 0.0000 | -2021.876 | 0.00000 | locust_intake_target_through_time
mod2 | -1907.268 | 52.32555 | -1975.891 | 111.1021 | -1961.285 | 60.59134 | locust_intake_target_through_time
mod3 | -1900.026 | 59.56801 | -1919.240 | 167.7533 | -1918.193 | 103.68257 | locust_intake_target_through_time
mod4 | -1895.604 | 63.98977 | -1914.818 | 172.1751 | -1913.772 | 108.10433 | locust_intake_target_through_time
mod5 | -1898.807 | 60.78667 | -1912.532 | 174.4619 | -1911.981 | 109.89482 | locust_intake_target_through_time

All model selection criteria agreed that `mod1` was the best model

### Final Model Summary

    In [22]:

    summary(mod1)

    Family: Multivariate normal
    Link function:

    Formula:
    cumulative_total_c ~ sex + interval + treatment + s(id, bs = "re")
    cumulative_total_p ~ sex + interval + treatment + s(id, bs = "re")

    Parametric coefficients:
                      Estimate Std. Error z value Pr(>|z|)
    (Intercept)      0.0125868  0.0036504   3.448 0.000565 ***
    sexM            -0.0108931  0.0041592  -2.619 0.008818 **
    interval2        0.0066835  0.0025232   2.649 0.008077 **
    interval3        0.0160330  0.0025232   6.354 2.09e-10 ***
    interval4        0.0259670  0.0025232  10.291  < 2e-16 ***
    interval5        0.0350061  0.0025232  13.874  < 2e-16 ***
    treatmentnone   -0.0061015  0.0040947  -1.490 0.136202
    (Intercept).1    0.0022792  0.0014613   1.560 0.118829
    sexM.1          -0.0036943  0.0014086  -2.623 0.008726 **
    interval2.1      0.0005296  0.0014992   0.353 0.723909
    interval3.1      0.0034148  0.0014992   2.278 0.022739 *
    interval4.1      0.0071035  0.0014992   4.738 2.16e-06 ***
    interval5.1      0.0125817  0.0014992   8.392  < 2e-16 ***
    treatmentnone.1 -0.0009902  0.0013868  -0.714 0.475224
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

    Approximate significance of smooth terms:
              edf Ref.df Chi.sq p-value
    s(id)   17.59     20  484.7  <2e-16 ***
    s.1(id) 10.61     20  110.7   0.381
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

    Deviance explained = 73.4%
    -REML = -954.88  Scale est. = 1         n = 115

### Key Findings

- Sex Effect: Males consumed significantly less protein and carbohydrates than females.
- Time (Interval) Effect:
  - Protein consumption increased significantly in later intervals, especially interval 5 (`p < 0.001`).
  - Carbohydrate consumption also showed increases in later intervals, notably in intervals 3, 4, and 5.
- Fertilization Treatment: Locusts in plots without nitrogen fertilization consumed less protein (`p = 0.0285`), but the effect was not significant for carbohydrates.
- Random Effect: There was no significant individual effect (`s(id)`), indicating individual locusts did not contribute heavily to variation.
- Deviance Explained: The model explained 20.5% of the variation in nutrient consumption.

### Model outputs

I am going to get the estimated marginal means for `inter`,`date`, and `treatment:date` and output as a dataframe

I will write out the transformed raw data to disk

I will also write out to disk the model object

    In [23]:

    intake_target_exp_preds <- expand.grid(
        sex = unique(dat3$sex),
        interval = unique(dat3$interval),
        treatment = unique(dat3$treatment),
        id= unique(dat3$id)
    )

    # Predict carbohydrate consumption
    intake_target_exp_preds$pred_carb <- predict(mod1, newdata = intake_target_exp_preds, type = "response")[, 1]

    # Predict protein consumption
    intake_target_exp_preds$pred_protein <- predict(mod1, newdata = intake_target_exp_preds, type = "response")[, 2]

    intake_target_exp_preds

A data.frame: 460 × 6
sex | interval | treatment | id | pred_carb | pred_protein
--- | --- | --- | --- | --- | ---
<fct> | <fct> | <fct> | <fct> | <dbl> | <dbl>
M | 1 | none | 1 | 0.005251196 | 0.002089128
F | 1 | none | 1 | -0.005641895 | -0.001605175
M | 2 | none | 1 | 0.011934674 | 0.002618693
F | 2 | none | 1 | 0.001041583 | -0.001075610
M | 3 | none | 1 | 0.021284240 | 0.005503910
⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮
F | 3 | High | 23 | 0.01278549 | 0.001897175
M | 4 | High | 23 | 0.03361250 | 0.009280174
F | 4 | High | 23 | 0.02271940 | 0.005585871
M | 5 | High | 23 | 0.04265163 | 0.014758435
F | 5 | High | 23 | 0.03175854 | 0.011064132

    In [24]:

    # Now averaging arcoss species, plot, and cage
    intake_target_exp_preds2 <- intake_target_exp_preds |>
        group_by(treatment,interval) |>
        summarize(pred_carb = mean(pred_carb),
               pred_protein = mean(pred_protein))

    intake_target_exp_preds2

    `summarise()` has grouped output by 'treatment'. You can override using the
    `.groups` argument.

A grouped_df: 10 × 4
treatment | interval | pred_carb | pred_protein
--- | --- | --- | ---
<fct> | <fct> | <dbl> | <dbl>
High | 1 | 0.007140269 | 4.320264e-04
High | 2 | 0.013823747 | 9.615916e-04
High | 3 | 0.023173312 | 3.846809e-03
High | 4 | 0.033107225 | 7.535505e-03
High | 5 | 0.042146356 | 1.301377e-02
none | 1 | 0.001038814 | -5.581555e-04
none | 2 | 0.007722292 | -2.859025e-05
none | 3 | 0.017071858 | 2.856627e-03
none | 4 | 0.027005771 | 6.545323e-03
none | 5 | 0.036044901 | 1.202358e-02

    In [26]:

    write.csv(intake_target_exp_preds2,here('data/processed/field_cage_experiment/intake_target_through_time_predictions.csv'),row.names=FALSE)
    write.csv(dat3,here('data/processed/field_cage_experiment/intake_target_through_time_predictions_raw_data.csv'),row.names=FALSE)

    write.csv(model_selection_tables, here("output/field_cage_experiment/intake_target_through_time_model_selection_criteria_results.csv"),row.names=FALSE)

    # Combine the lists into a single list
    models <- list('mod1'= mod1, 'mod2' = mod2, 'mod3'= mod3, 'mod4'= mod4, 'mod5'= mod5)

    # Save the combined list to an RDS file
    saveRDS(models, here("output/field_cage_experiment/model_objects/intake_target_through_time_models.rds"))

# Section 3 - Final locust mass analysis

This model looks to see if the final mass of locusts differed between treatments.

    In [27]:

    locust_dat <- read_csv(here('data/raw/field_cage_experiment/LocustCageData.csv'),show_col_types = FALSE) |>
        clean_names() |>
        filter(!complete.cases(oedaleus_australis)) |>
        select(treatment,cage_num,mass_g,carb_mg_mg,protein_mg_mg,sex) |>
        drop_na() |>
        mutate(across(c(sex,cage_num),~factor(.)),
              treatment = factor(case_when(treatment == 'HighN' ~ 'high',
                                    treatment == 'MedH' ~ 'medium',
                                    TRUE ~ 'none'), levels= c('none','medium','high')))

    New names:
    • `` -> `...24`
    • `` -> `...25`
    • `` -> `...26`
    • `` -> `...27`
    • `` -> `...28`

    In [28]:

    locust_dat |> glimpse()

    unique(locust_dat$sex)
    unique(locust_dat$treatment)
    unique(locust_dat$cage_num)

    Rows: 292
    Columns: 6
    $ treatment     <fct> high, high, high, high, high, high, high, high, high, hi…
    $ cage_num      <fct> 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,…
    $ mass_g        <dbl> 0.242, 0.314, 0.288, 0.257, 0.181, 0.115, 0.132, 0.336, …
    $ carb_mg_mg    <dbl> 0.1210000, 0.1210000, 0.1210000, 0.1210000, 0.1210000, 0…
    $ protein_mg_mg <dbl> 0.1729659, 0.1729659, 0.1729659, 0.1729659, 0.1729659, 0…
    $ sex           <fct> F, F, F, F, M, M, M, F, F, M, M, M, F, F, F, F, F, F, F,…

1. F
2. M

Levels:

1. 'F'
2. 'M'

1. high
2. none
3. medium

Levels:

1. 'none'
2. 'medium'
3. 'high'

1. 1
2. 2
3. 3
4. 4
5. 5
6. 6
7. 7
8. 8
9. 9
10. 10
11. 11
12. 12
13. 13
14. 14
15. 15
16. 16
17. 17
18. 18
19. 19
20. 20
21. 21
22. 22
23. 23
24. 24
25. 25
26. 26
27. 27
28. 28
29. 29
30. 30
31. 31
32. 32
33. 33
34. 34
35. 35
36. 36

Levels:

1. '1'
2. '2'
3. '3'
4. '4'
5. '5'
6. '6'
7. '7'
8. '8'
9. '9'
10. '10'
11. '11'
12. '12'
13. '13'
14. '14'
15. '15'
16. '16'
17. '17'
18. '18'
19. '19'
20. '20'
21. '21'
22. '22'
23. '23'
24. '24'
25. '25'
26. '26'
27. '27'
28. '28'
29. '29'
30. '30'
31. '31'
32. '32'
33. '33'
34. '34'
35. '35'
36. '36'

## Building model

    In [29]:

    mod <- gam(mass_g ~
                   sex +
                   treatment +
                   s(carb_mg_mg,protein_mg_mg) +
                   s(cage_num,bs='re'),
            family=scat(),
            select=TRUE,
            data=locust_dat)

## Model appraisal

    In [30]:

    gratia::appraise(mod)

## Model summary

    In [31]:

    summary(mod)

    Family: Scaled t(3.062,0.036)
    Link function: identity

    Formula:
    mass_g ~ sex + treatment + s(carb_mg_mg, protein_mg_mg) + s(cage_num,
        bs = "re")

    Parametric coefficients:
                     Estimate Std. Error z value Pr(>|z|)
    (Intercept)      0.326188   0.007210  45.243   <2e-16 ***
    sexM            -0.147694   0.005507 -26.820   <2e-16 ***
    treatmentmedium  0.014935   0.009518   1.569    0.117
    treatmenthigh   -0.011055   0.010095  -1.095    0.273
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

    Approximate significance of smooth terms:
                                      edf Ref.df Chi.sq  p-value
    s(carb_mg_mg,protein_mg_mg)  0.001665     29  0.002    0.416
    s(cage_num)                 17.399299     33 42.160 1.78e-05 ***
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

    R-sq.(adj) =  0.617   Deviance explained = 61.6%
    -REML = -427.13  Scale est. = 1         n = 292

## Key Findings

- Sex Effect: Males are significantly lighter than females (`p < 0.001`).
- Fertilization Treatment:
  - No significant effect of medium or high fertilization treatments on locust mass (`p > 0.1`).
- Nutrient Content: The interaction between carbohydrate and protein content did not significantly affect locust mass (`p = 0.416`).
- Cage Effect: Significant variability in mass was attributed to differences between cages (`p = 1.78e-05`).
- Deviance Explained: The model explains 61.6% of the variance in locust mass.

## Model outputs

I dont need the estimated marginal means since there was no difference. I will just save the model to disk?

    In [32]:

    saveRDS(mod, here("output/field_cage_experiment/model_objects/final_locust_weight_by_field_cage_treatment_models.rds"))