I am having a problem using the tryCatch() function in R in a function I created.
What I want to do is this:
gamlss modelpredict function to extract model predictions over a new range of valuesMy main problem is that my model is somewhat unstable and once in a while predictions are kind of wild, which in turn generates an error when I try to analyze it with gamlss. My objective is to write a tryCatch statement within my simulation function and to basically simply run the simulation/prediction code a second time in the event that an error occurs. (I know this is not optimal, I could also write it in a recursive statement using repeat for example and run it until I don't get an error but I get few enough errors that the probability of getting two in a row is quite low, and I'm having enough troube with this task as it is.)
So I simplified my code as much as I could and created a dummy dataframe for which the modelling still works.
I wrote in the code where I believe the error is (with the predict function which does not find the mod_sim object). It is likely there since the cat just above this line prints while the one just below doesn't print.
I think there are some things about how tryCatch works that I don't understand well enough and I'm having a hard time to understand which objects are kept in which parts of functions and when they can be called or not...
Here is the code I have so far. The error occurs at l.84 (identified in the script). The data and code can be found here.
library(tidyverse)
library(gamlss)
library(gamlss.dist)
#Load data
load('DHT.RData')
#Run original model
mod_pred<-gamlss(harvest_total ~ ct,
data = DHT,
family = DPO)
#Function to compute predictions based on model
compute_CI_trad_gamlss<-function(n.sims=200, mod){#,
#DF for simulations
df_sims<-as.data.frame(DHT)
#Dateframe with new data to predict over
new.data.ct<<-expand.grid(ct=seq(from=5, to=32, length.out=50))
#matrix to store predictions
preds.sim.trad.ct <<- matrix(NA, nrow=nrow(new.data.ct), ncol=n.sims)
#Number of obs to simulate
n<-nrow(df_sims)
#Simulation loop (simulate, analyze, predict, write result)
for(i in 1:n.sims){
#Put in tryCatch to deal with potential error on first run
tryCatch({
#Create matrix to store results of simulation
y<-matrix(NA,n,1)
#in DF for simulations, create empty row to be filled by simulated data
df_sims$sim_harvest<-NA
#Loop to simulate observations
for(t in 1:n){
#Simulate data based on model parameters
y[t]<-rDPO(n=1, mu=mod$mu.fv[t], sigma = mod$sigma.fv[t])
}#enf of simulation loop
#Here I want the result of the simulation loop to be pasted in the df_sims dataset
df_sims$sim_harvest<-y
#Analysis of simulated data
mod_sim<-gamlss(sim_harvest ~ ct,
data = df_sims,
family = DPO)
#Refit the model if convergence not attained
if(mod_sim$converged==T){
#If converged do nothing
} else {
#If not converged refit model
mod_sim<-refit(mod_sim)
}
cat('we make it to here!\n')
#Store results in object
ct <<-as.vector(predict(mod_sim, newdata = new.data.ct, type='response'))
cat('but not to here :( \n')
#If we made it down here, register err as '0' to be used in the if statement in the 'finally' code
err<<-0
},
#If error register the error and write it!
error = function(e) {
#If error occured, show it
cat('error at',i,'\n')
#Register err as 1 to be used in the if statement in the finally code below
err<<-1
},
finally = {
if(err==0){
#if no error, do nothing and keep going outside of tryCatch
}#End if err==0
else if (err==1){
#If error, re-simulate data and do the analysis again
y<-matrix(NA,n,1)
df_sims$sim_harvest<-NA
#Loop to simulate observations
for(t in 1:n){
#Simuler les données basées sur les résultats du modèle
y[t]<-rDPO(n=1, mu=mod$mu.fv[t], sigma = mod$sigma.fv[t])
}#enf of simulation loop
#Here I want the result of the simulation loop to be pasted in the df_sims dataset
df_sims$sim_harvest<-y
#Analysis of simulated data
mod_sim<-gamlss(sim_harvest ~ ct,
data = df_sims,
family = DPO)
cat('we also make it here \n')
#Store results in object
ct <<-as.vector(predict(mod_sim, newdata = new.data.ct, type='response'))
cat('but not here... \n')
}#End if err==1,
}#End finally
)#End tryCatch
#Write predictions for this iteration to the DF and start over
preds.sim.trad.ct[,i] <<-ct
#Show iteration number
cat(i,'\n')
}
#Do some more stuff here
#Return results
return(preds = list(ct= list(predictions=preds.sim.trad.ct)))
}
#Run simulation and store object
result<-compute_CI_trad_gamlss(n.sims=20, mod=mod_pred)
Anyway I hope someone can help!
Thanks a lot!
So after a bit of trial and error I managed to make it work. I believe the problem lies in the mod_sim object that is not saved to the global environment. predict (or predict.gamlss here) is probably not looking in the function environment for the mod_sim object although I don't understand why it wouldn't. Anyway using <<- (i.e. assigning the object in the global environment from the function) for every object created in the function seemed to do the trick. If anyone has an explanation on why this happens though I'd be glad to understand what I'm doing wrong!