Further model implementation

[ScymnusRIP]

Is there a way to add one or more custom models to the available models or to the "fitting many models" routine? I'm imaging two ways:

• provinding all necessary data in order to let you improve the package version with X models
• letting users to add custom models to the routine, something like what available in termPerf package at following link Thanks in advance

[padpadpadpad]

Hi ScymnusRIP

In thermPerf a method of adding custom models is needed because the fitting process is done within the thermPerf package.

rTPC is basically a set of helper functions to fit the most common rTPC models, but the fitting process is done using nls.multstart and other packages. Consequently if you just wrote your own function, it would like work with the "fitting many models" routine, apart from you would have to set your own start values and lower and upper limits. I think calc_params() should also work out the box as it just uses the model object to calculate the extra traits.

That being said if you did want to add a function in I would happily help you add it. Things we would need are:

• a function for the model (named for the year of first use in the literature)
• know its sensible upper and lower limits (based either on the data or on parameter values from the literature) and these added to get_lower_lims() and get_upper_lims()
• know how to get sensible start values (based either on the data or on parameter values from the literature) and add these to get_start_vals()
• add the model name to get_model_names()

Happy to help where needed.

[Shuangen-YU]

Hello, your rTPC is really helpful. During the process of studying thermal performance curves, I came across another highly recognized model. Lately, I've been eager to use it. If possible, I suggest considering its integration into rTPC. This model has a Gaussian curve before Topt and a parabolic curve after Topt. The link to the article is https://doi.org/10.1073/pnas.0709472105. Thanks in advace

[padpadpadpad]

Hi @Shuangen-YU have you got this model as a function in R or have you fitted it in R before?

[Shuangen-YU]

Yes, I have tried to calculate the parameters of this model using nonlinear fitting, but due to the limited amount of data I used for my research, I haven't fully implemented it at the moment.

In my codes, I assume that Topt is known. However, in most cases, we are unsure about where Topt is. This is also where I am currently puzzled. Perhaps, with a sufficient number of data points, is there a way to calculate Topt similar to the other models used in rTPC?

Partial codes:

gaussian (before Topt)

gaussian <- function(x, a, b, c) { a exp(-((x - b)/(2c))^2) }

parabola (after Topt)

parabola <- function(x, a, b, c) { a * (x - b)^2 + c }

input data

x <- c(9, 13, 17, 21, 25, 28) y <- c(0.007, 0.01, 0.025, 0.028, 0.009, 0)

here, I assume that Topt is 20.09305

However, In most cases, we are unsure about what Topt is.

gauss_data <- data.frame(x = x[x <= 20.09305 ], y = y[x <= 20.09305 ]) gauss_model <- nls(y ~ gaussian(x, a, 20.09305 , c), data = gauss_data, start = list(a = 1, c = 1), trace = FALSE)

parabola_data <- data.frame(x = x[x >= 20.09305 ], y = y[x >= 20.09305 ]) parabola_model <- nls(y ~ parabola(x, a, 20.09305 , gauss_coef[1]), data = parabola_data, start = list(a = 1), trace = FALSE)

segmented_function <- function(x) { ifelse(x <= 20.09305 , predict(gauss_model, newdata = data.frame(x = x)), predict(parabola_model, newdata = data.frame(x = x))) }

plot_data <- data.frame(x = x, y = y) aaa <- ggplot(plot_data, aes(x, y)) + geom_point() + stat_function(fun = segmented_function, color = "blue") + geom_text(x = min(x), y = max(y), label = paste("Gaussian R^2 =", round(gauss_r_squared, 3)), color = "blue", hjust = 0, vjust = 1) + geom_text(x = min(x), y = max(y), label = paste("Parabola R^2 =", round(parabola_r_squared, 3)), color = "red", hjust = 0, vjust = 3) + geom_vline(xintercept = 20.09305, linetype = "dashed", color = "red") + labs(x = "Temperature", y = "Growth rate kg/month") + theme_bw()+ geom_vline(xintercept = 20.09305, linetype = "dashed", color = "red")

[padpadpadpad]

I have had a first go here:

# attempted deutsch function
deutsch_2008 <- function(temp, a, topt, ctmax){
  est <- {ifelse(temp < topt, exp(-((temp - topt)/(2*a))^2), (1 - ((temp - topt)/(topt - ctmax))^2)  )}
  return(est)
}

# test it out
temp = 1:40
a = 5; topt = 30; ctmax = 40
rate = deutsch_2008(temp, a, topt, ctmax)

# plot it
plot(rate ~ temp)

Seems to work fine. I am unsure what the benefits are of this curve over others though? For suppose I could do with some data on what values a ($\sigma_{p}$) in the Deutsch paper can take for setting start values and sensible limits.

Also the implementation in the Deutsch paper has no way for the fit to move up and down (I dont think), so the rmax is always 1. Can easily add an rmax term in to help with that though.

[Shuangen-YU]

Your test works fine. The difference between this curve and others (e.g., rezende_2019 in rTPC, another segmented model) is that Deutsh's model used gaussian curve before Topt. This seems to be accepted by many famous physiologists. This formulation of performance conforms to both observed and theoretical considerations. Of course, other research also has a high level of recognition.

I think the sigma_p is a shape parameter determing the steepness of the rising portion of the curve. In the study, they defined the asymmetry, α, of a performance curve as the ratio of its width on the warm and cold sides of Topt. α = (CTmax-Topt)/(4*sigma_p). They used a range of target values of α from 2 to 8, and obtained a family of plausible-fit curves. At last they used α =3 for tupical insects.

Well, I am agree with you on the point that we can add an rmax. In the Deutsch paper, they used the relative fitness (not the ture rate, I think).

[padpadpadpad]

@Shuangen-YU I have added deutsch_2008() with a parameter for rmax to allow the height of the curve to change. I found it difficult to work out alpha was doing in the equation, so for the start values i set them to the max - min temperatures / 5. Seems to work but can be tweaked.

Let me know how you get on.