drin_hmc: Developmental Rate using HMC

In mkondakis/DRIN: Inferring Developmental rate thermal thresholds statistically in a Bayesian framework

Description

Estimation of Ecological model parameters using HMC sampling

Usage

drin_hmc(
  c_data,
  sparam = c(4, 11000, 10000, 1),
  mtype = "bieri",
  lik = "gauss",
  prior = "NULL",
  init_list = "NULL",
  predplot = FALSE,
  ...
)

Arguments

sparam	vector indicating the number of chains, iterrations, burn_in and thinning in Rstan execution. Default is c(4,11000,10000,1);
mtype	is the ecological model type between "bieri" (default), "briere", "analytis" and "lactin";
lik	is the likelihood choice between "gauss" (default), "igamma" (inverse gamma) in rstan form;
prior	is a list with parameter priors. If not declared default weekly informative;
init_list	is the initial values of the parameters for the shake of the HMC algorithm. It should be a list or a function of a list
C_data	data list including: y the response values, x the predictor values, status of censoring 0 if anthropod is not evolved ; NP is the number predictions to be made.

Value

List that contains: (est) s4 type Stanfit object, (summary) Stanfit summary of basic parameters, (diagnostics) Stanfit overal diagnostics report

See Also

stan

Examples

## Not run: 
Do not run
DataBoot<-xls::read.xlsx(paste(bugname,".xlsx",sep=""), sheetIndex=1, header=FALSE, endRow = 1000)
DataBoot<-na.omit(DataBoot)
if (inum==1) {colnames(DataBoot)<- c("Temp","gr","status")} # if censored zero data exist
else {colnames(DataBoot)<- c("Temp","gr")} # if no zeros exist
y<-c(DataBoot$gr) #Observed days until development ocurrence
temp<-c(DataBoot$Temp) #create predictor for observed values (i.e. temperatures)
N <-dim(DataBoot)[1] #number of data
if (inum==1) status<-c(DataBoot$status) else status<-rep(1,N)
y<-1/y # create rensponse variable of observed developmental rates
if (inum==1) {data <- list(N = N, y = y, x=temp, status= status, t=10.0)} else {data <- list(N = N, y = y, x=temp, status= status, t=10.0)}  
data$NP<-round(2*(max(data$x)-min(data$x)),0) #number of predictions
data$xpred<-seq(min(data$x),max(data$x),length.out = data$NP) #generate NP predictor values for model predictions
modbieri_stan<-drin_hmc(data,c(1,11000,10000,1)) #run the developmental 
modbieri_stan<-drin_hmc(data,c(1,11000,10000,1),"bieri","gauss")

## End(Not run)

mkondakis/DRIN documentation built on Dec. 21, 2021, 7:06 p.m.