R code for looking at how removing data points impacts a Sharpe-Schoolfield model fit.

dirty_TPC_censor_code.R

# have a look at whether having multiple 0 values impacts the results of model fitting to TPCs

# load libraries
librarian::shelf(rTPC, nls.multstart, tidyverse, broom)

# simulate data using Sharpe Schoolfield model
d_simulate <- tibble(temp = seq(from = 10, to = 55, length.out = 10),
                     rate = sharpeschoolhigh_1981(temp = temp, r_tref = 0.5, e = 0.6, eh = 5, th = 33, tref = 17.5)) %>%
  # add noise
  mutate(noise = rnorm(n(), 0, 0.1),
         rate2 = rate + noise,
         # if rate is < 0, set to 0
         rate2 = ifelse(rate2 < 0, 0, rate2))

# plot rates
ggplot(d_simulate, aes(temp, rate2)) +
  geom_point(size = 3) +
  theme_bw(base_size = 16)

# fit sharpe schoolfield to whole dataset using nls_multstart
fit <- nls_multstart(rate2 ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref = 17.5),
                     data = d_simulate,
                     iter = 100,
                     start_lower = get_start_vals(d_simulate$temp, d_simulate$rate2, model_name = 'sharpeschoolhigh_1981') - 5,
                     start_upper = get_start_vals(d_simulate$temp, d_simulate$rate2,model_name = 'sharpeschoolhigh_1981') + 5, 
                     supp_errors = 'Y')

# plot results
pred1 <- augment(fit)

# fit model to only temperatures <=45ºC
d_simulate2 <- d_simulate %>%
  filter(temp <= 45)

# fit sharpe schoolfield to whole dataset using nls_multstart
fit2 <- nls_multstart(rate2 ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref = 17.5),
                     data = d_simulate2,
                     iter = 100,
                     start_lower = get_start_vals(d_simulate$temp, d_simulate$rate2, model_name = 'sharpeschoolhigh_1981') - 5,
                     start_upper = get_start_vals(d_simulate$temp, d_simulate$rate2,model_name = 'sharpeschoolhigh_1981') + 5, 
                     supp_errors = 'Y')

# plot results
pred2 <- augment(fit2)

# fit model to only temperatures <=40ºC
d_simulate3 <- d_simulate %>%
  filter(temp <= 40)

# fit sharpe schoolfield to whole dataset using nls_multstart
fit3 <- nls_multstart(rate2 ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref = 17.5),
                      data = d_simulate3,
                      iter = 100,
                      start_lower = get_start_vals(d_simulate$temp, d_simulate$rate2, model_name = 'sharpeschoolhigh_1981') - 5,
                      start_upper = get_start_vals(d_simulate$temp, d_simulate$rate2,model_name = 'sharpeschoolhigh_1981') + 5, 
                      supp_errors = 'Y')

# plot results
pred3 <- augment(fit3)

# plot the results
ggplot() +
  geom_point(data = d_simulate, aes(temp, rate2), size = 3) +
  geom_line(data = pred1, aes(temp, .fitted), color = 'red') +
  geom_line(data = pred2, aes(temp, .fitted), color = 'blue') +
  geom_line(data = pred3, aes(temp, .fitted), color = 'green') +
  theme_bw(base_size = 16)