Reputation: 13
I am working with biological data, like growth responses to temperature changes or nutrient concentrations and stuff like that. My goal is to build an automated routine to load data, fit different simple, non-linear models onto the same data set and write the output (i.e. model coefficients) into an object, so that i can afterwards judge, which model fits best. For this is have written to following code:
#make up some example data
data<-data.frame(x=c(0,1,2,3,4,5,6,7,8,9,10),
y=c(0,2,4,5.5, 6.2, 6.4, 6.5, 6.6, 6.6,6.65,6.65))
#plot data
ggplot() +
geom_point(data = data, aes(x = x, y = y))+
scale_colour_discrete(type = getOption("ggplot2.discrete.colour"))
##########
# For this example, i want to fit the Michaelis-Menten formula: y=((Vmax*x)/(Km+x)
#First, i try to fit the data with the bare function:
MiMefit<-nls(y~(Vmax*x)/(Km+x),
data=data,
start = c(Vmax=5, Km=2))
#This seems to work, and I can predict & visualize modelled values
MiMefit
data$fit<-predict(MiMefit)
ggplot() +
geom_point(data = data, aes(x = x, y = y))+
geom_line(data = data, aes(x = x, y = fit), color="red")+
scale_colour_discrete(type = getOption("ggplot2.discrete.colour"))
##########
#Since i ´do not want to write every model by hand in every fitting procedure, i deposit the model #in an object and define the Michaelis-Menten equation as "standalone" model
MiMe <- function(x, Vmax, Km) {
return((Vmax*x)/(Km+x))
}
#test for functionality is positive
MiMe(x=4, Vmax = 100, Km = 9)
#Now i can do the fitting by calling the model name instead of writing the function:
#MiMefit2<-nls(y~MiMe(x = data$x, Vmax, Km),
data=data,
start = c(Vmax=5, Km=2))
MiMefit2
##################
#Since i want to automate this, i want to call the model by a *generic* name, that will later
#(in a for loop) "encode" other models with other equations.
#I make up an object that links to the MiMe formula, called "model"
model<-MiMe
#seems to have worked:
str(model)
#Now i can do the fitting with using the generic name:
MiMefit3<-nls(y~model(x = data$x, Vmax, Km),
data=data,
start = c(Vmax=5, Km=2))
MiMefit3
###################
#Still, i need to specify the parameters of the model (Vmax, Km) manually. Other models have other #parameter names, and diff. number of parameters. Thus i´d like to deposit the names of the #parameters in a vector:
params<-c("Vmax","Km")
str(params)
My problem is now, that i can not manage to edit the command in such way that i can parse the
# necessary parameter names to the function:
MiMefit4<-nls(y~I(model(x = data$x, params)),
data=data,
start = c(Vmax=5, Km=2))
You can see that parsing the parameters to the function using the vector fails
I checked all articles i could find on nls, but could not find other examples of my problem. It also seems to that it has to do with the syntax required by functions, but my searches in there were not successful either.
I have tried for days now to rephrase the command e.g.
MiMefit4<-nls(y~I(model(x = data$x, params[1], params[2])),
data=data,
start = c(Vmax=5, Km=2))
or
MiMefit5<-nls(y~I(model(as.list(noquote(paste("x = data$x","," ,noquote(params[1]),",",noquote(params[2])))))),
data=data,
start = c(Vmax=5, Km=2))
MiMefit5
without success. I have tried the noquote command: noquote(params[1]); I tried diverse as.xxxx functions, trying to make it work with as.formula, as.function, as.list and so on. I think that a list should somehow be expected, but that did not work.
Also, ideally the vector with the params is called once,inserting all contained parameters into the list of model arguments, since in the future, as stated, i´d like to specify models with different numbers of parameters.
As always, thanks for your time and effort trying to help with my problem!! Sebastian
EDIT:
Thanks for your reply. Your approach works, but does not really fit into my larger scheme of automization. I am using the package rTPC and the function nls_multstart, which is a derivative of nls as far as i can see. The package contains predesigned models already and deposits them as well as a vector/list:
#these are the packages in case you wanna try
#remotes::install_github("padpadpadpad/rTPC")
library(rTPC)
#install.packages("nls.multstart")
library(nls.multstart)
#the package rTPC contains many models:
get_model_names()
#I can use object 'model' to link to these:
model<-get_model_names()[5]
model
str(model)
get(model)
#I use the below command to estimate start values and limits:
start_vals<-get_start_vals(x = data$x, y=data$y, model_name = model)
lower_lims<-get_lower_lims(x = data$x, y=data$y, model_name = model)
upper_lims<-get_upper_lims(x = data$x, y=data$y, model_name = model)
#Since the package automatically makes named numeric vectors, it #deposits also the names of the parameters:
params<-(c("x", names(lower_lims)))
# Here, i add the x value!
params
#the classical,i.e., manual approach from the author of the package:
fit<-nls_multstart(formula =y~beta_2012(x,a,b,c,d,e)
,data = d
,start_lower = lower_lims
,start_upper = upper_lims
,iter = 500
,supp_errors = "Y")
fit
#to automate it and make it suitable for a for loop, i want to generalize the fit expression.
#i found that with get(), i can make it work halfway:
fit2<-nls_multstart(formula =y~get(model)(x,a,b,c,d,e)
,data = d
,start_lower = lower_lims
,start_upper = upper_lims
,iter = 500
,supp_errors = "Y")
fit2
Now i am still struggling to implement the latter part, i.e. the supply of arguments. I have tried several formulations of the formula call:
formula =y~get(model)(x,a,b,c,d,e) WORKS
formula =y~get(model)(params)
formula =y~get(model)(paste(params))
formula =y~get(model)(paste0(params))
formula =y~do.call(model,params)
formula =y~do.call(what = get(model),args = params)
formula =y~get(model)paste(as.character(params), collapse=", ")
formula =y~get(model)paste(as.character(params), sep="' '", collapse=", ")
All these without any success. Maybe you have a better idea about the nature of that argument and how it needs to be parsed and how you can make that from
Upvotes: 1
Views: 199
Reputation: 18581
To make things easier, I think its best to create one function that returns a formula and feed a list of parameters into that function using do.call
.
# make up some example data
data<-data.frame(x=c(0,1,2,3,4,5,6,7,8,9,10),
y=c(0,2,4,5.5, 6.2, 6.4, 6.5, 6.6, 6.6,6.65,6.65))
model <- function(lhs, x, Vmax, Km) {
rhs <- paste0("(", Vmax, "*", x, ")/(", Km, "+", x, ")")
reformulate(rhs, lhs)
}
params <- list(lhs = "y",
x = "x",
"Vmax",
"Km")
str(params)
#> List of 4
#> $ lhs: chr "y"
#> $ x : chr "x"
#> $ : chr "Vmax"
#> $ : chr "Km"
MiMefit4<-nls(do.call("model", params),
data = data,
start = c(Vmax=5, Km=2))
MiMefit4
#> Nonlinear regression model
#> model: y ~ (Vmax * x)/(Km + x)
#> data: data
#> Vmax Km
#> 8.503 2.051
#> residual sum-of-squares: 1.629
#>
#> Number of iterations to convergence: 7
#> Achieved convergence tolerance: 2.723e-06
Created on 2023-04-26 with reprex v2.0.2
Upvotes: 1