#Fit the rumen model from script 3 to observed duodenal flow data #Execute Load Observed Data.R before running this script #Load the model function from execise_1 R script 3 #Load the Obj.f and pred.f functions at the bottom of this script require(plyr) require(data.table) require(FME) #Create a vector of Initial Model Inputs and Parameters parameters <- c(fDMI = 20, #kg DM/h cCp = 0.196*0.644, #kg/kg DM, 0.53+.114 = .644 from publication cCho = 0.50, KdCp = 0.04 * 24, #per day KdCho = 0.03 * 24, Kp = 0.05 * 24, #per hour KCpAbs = 0.75, KChoAbs = 0.45, KRdp = 0.25, VmMiG = 0.25) #test the parameter definitions using the model function Gut(t,parameters) #A list of model parameters to be derived from the data (can be more than 1) fitParmNames <- c("Kp") #fitParmNames <- c("KdCho") #fitParmNames <- c("Kp","KdCho") #fitParmNames <- c("Kp","KdCho","KRdp") #fitParmNames <- c("Kp","KdCho","KRdp","VmMiG") #if we assume Kd for CP from in situ is correct, Kp can be derived to achieve the correct FCpSI. #If ruminal concentrations of Cp were provided, the Kd and Kp could both be solved. #If Kp is defined by FCpSi and applies to CHO as well, then KdCho can be derived from FChoSI. # #Create a vector of initial parameter estimates for optimization parms_init <- parameters[fitParmNames] #upper and lower bounds that can be used. Must have the same names as parameters lb <- parameters[fitParmNames] #pick up the names #lb[] <- rep(0,2) #set all to 0 #or specify each lb['Kp'] <- 0.01 * 24 #lb['KdCho'] <- 0.01 * 24 #lb['KRdp'] <- 0.0 #lb['VmMiG'] <- 0.0 ub <- parameters[fitParmNames] #pick up the names #ub[] <- rep(0.1,2) #set all to 0.1 #or specify each ub['Kp'] <- 0.08 * 24 #ub['KdCho'] <- 0.15 * 24 #List of model variables to be compared to observed. Var names and units must match between model and obs obsVars <- c("FCpSi") #obsVars <- c("FChoSi") #obsVars <- c("FCpSi","FChoSi") #obsVars <- c("FCpSi","FCpMiSi","FChoSi") scale_res <- TRUE #scale residuals to the mean observed value of each #List the model output time points to be compared to observed data tout <- c(4) #if a single steady state time is required, could specify as the length of t. #Create a function for running model predictions that can be called with teh latest parameters estimates pred.f = function(parameters) { #predict from observed dietary inputs and P using apply function pred <- apply(o[,c('fDMI', 'Dt_CP','Dt_CHO','KdRUP')], 1, function(x) { #extract diet inputs to x and call fn parameters['fDMI'] <- x['fDMI'] #update inputs from each row of data parameters['cCp'] <- x['Dt_CP']/100 parameters['cCho'] <- x['Dt_CHO']/100 parameters['KdCp'] <- x['KdRUP'] out <- Gut(t,parameters) #call the model and collect output pred <- out[out[,"time"] == tout, obsVars] #retrieve only the specified times and obsVars return(pred) #return the predicted values #return(list(pred=pred,pars=parameters)) #use this statement only to check that diet vals are being changed }) return(pred) } pred.f(parameters) #The model objective function; generates a vector of residual errors to be minimized #Uses the pred.f function to generate predicted values Obj.f <- function(P) { #update the parameters vector with values passed from the optimizer (P) parameters[names(P)] <- P #extract the observed values obs <- as.matrix(o[, obsVars]) colnames(obs) <- obsVars obsmean <- sapply(as.data.frame(obs), FUN = mean) obst <- melt(as.matrix(obs), id.vars=obsVars, var2.name="name", value.name="obs") obst <- rename(obst, c("Var2"="name")) #simulate predicted from observed dietary inputs and P using apply function pred <- as.matrix(pred.f(parameters)) #above creates a column vector for 1 var but a wide matrix for 2 or more if(dim(pred)[1] == length(obsVars)) {pred <- t(pred)} #convert wide matrix to long colnames(pred) <- obsVars predt <- melt(as.matrix(pred), id.vars=obsVars, var2.name="name", value.name="mod") predt <- rename(predt, c("Var2"="name")) residuals <- merge(obst,predt, sort=FALSE) residuals$res.unweighted <- residuals$mod - residuals$obs #Calc unscaled resid errors var_indx <- as.vector(match(residuals$name,obsVars)) residuals$scale <- as.vector(obsmean[var_indx]) if(scale_res==FALSE) residuals$scale <- rep(1,length(residuals$name)) residuals$res <- residuals$res.unweighted * residuals$scale #res <- unlist(res, use.names = FALSE) res <- residuals[,c("res")] return(res) #return residuals in res } Obj.f(parms_init) #test the objective function using parms_init sum(Obj.f(parms_init)) #Fit the model to the data; be patient. Must load the functions below first. #It is not a speedy process in R. wait for the > to return. scale_res <- TRUE #scale residuals to the mean observed value of each #scale_res <- FALSE #scale residuals to the mean observed value of each m1 <- modFit(f=Obj.f,p=parms_init, lower=lb) #m1 <- modFit(f=Obj.f, p=parms_init, lower=lb, upper=ub) summary(m1) m1$par/24*100 #express the derived rates as % per hour for literature reference attributes(m1) #Transfer the final parameter estimates to the parameters vector and run the model #this section needs further work. pred for 1 variable is a matrix. For 2+, it needs to be transposed. parameters[fitParmNames] <- m1$par pred <- as.vector(pred.f(parameters)) #pred <- as.data.frame(t(pred.f(parameters))) #had problems with this statement when fitting to 1 var obs <- o[, obsVars] studies <- o[, "PubID"] res <- pred-obs #plot(pred,res) #plot(pred,obs) resprplt(obs,pred,"RUP") #plot obs and res vs pred resprgrpplt(obs,pred, studies,"RUP") #plot obs and res vs pred with lines w/in PubID #resprplt(obs$FCpSi,pred$FCpSi,"RUP") #resprplt(obs$FChoSi,pred$FChoSi,"CHO") #resprplt(obs$FCpMiSi,pred$FCpMiSi,"MiP") RMSE(obs,pred) #use for a single predicted variable RMSE(obs$FCpSi,pred$FCpSi) RMSE(obs$FChoSi,pred$FChoSi)