# RUSSU 2009 Emax
# C = (e^(ALPHA)*D / (e^(BETA)+D)) (1 + EPS)
# ALPHA = THETA1+ETA_1
# BETA = THETA2+ETA_2

# This model is used for estimation.
#
# Limitations:
# - The current implementation does not include prior distributions and Bayesian estimation
# - DOSE has been replaced by TIME for compatibility reasons for NONMEM
#
# Compatibility:
# NONMEM and Monolix
#
# Dataset:
# simulated dataset provided by the authors
#
# Evaluation method for this model:
# run NONMEM estimation via TEL
# run MONOLIX estimation via TEL
# compare results between the two tools and with the figures in the original publication
# (although results were obtained via Bayesian estimation and thus they will not be identical)
#
# By Elisa Borella, Letizia Carrara, Silvia Lavezzi, Paolo Magni, Giulia Muzio, Lorenzo Pasotti, Elena Tosca
# BMS lab - Univ Pavia
# Copyright 2015
#
# THIS IS A BAYESIAN VERSION OF THE MODEL

russu_Emax_data = dataObj{

DECLARED_VARIABLES{Y::observation}

   DATA_INPUT_VARIABLES{
      ID:{use is id}
      TIME:{use is idv}
      DV:{use is dv, variable=Y} # CMAX
   }
   
   SOURCE {
   	srcfile : {
       file="Simulated_russu_2009_emax_data.csv",  
       inputFormat is nonmemFormat}
       }
   }

russu_pri = priorObj{
	
	PRIOR_PARAMETERS{
	
	  # values have to be updated
	  MU_THETA1 = 1
      MU_THETA2 = 1
      VAR_THETA1 = 1
      VAR_THETA2 = 1

	  # in the original publication, b are provided as rates
      a_1 = 1
      b_1 = 1
      a_2 = 1
      b_2 = 1
      a_cv = 1
      b_cv = 1
	}
	
	PRIOR_VARIABLE_DEFINITION{
	  THETA1 ~ Normal(mean=MU_THETA1, var=VAR_THETA1)
	  THETA2 ~ Normal(mean=MU_THETA2, var=VAR_THETA2)
	  invOMEGA_1 ~ Gamma1(shape=a_1, scale=b_1)
	  invOMEGA_2 ~ Gamma1(shape=a_2, scale=b_2)
	  invCV2 ~ Gamma1(shape=a_cv, scale=b_cv)
	  
	  OMEGA_1 = 1/invOMEGA_1
	  OMEGA_2 = 1/invOMEGA_2
	  CV = sqrt(1/invCV2)
	}
}

russu_Emax_mdl = mdlObj{

 	IDV { T }
 	
 	VARIABILITY_LEVELS{
   	DV : { level=1, type is observation }
   	ID : { level=2, type is parameter }
   }

	
	STRUCTURAL_PARAMETERS{
      THETA1 
      THETA2 
	  CV  
	}
	
	VARIABILITY_PARAMETERS{
	  OMEGA_1 
	  OMEGA_2 
	}
		
   	
   	INDIVIDUAL_VARIABLES {
	   	ALPHA : { type is linear, pop = THETA1, ranEff = ETA_1 }
		BETA : { type is linear, pop = THETA2, ranEff = ETA_2 }
   	}
   	
   	 RANDOM_VARIABLE_DEFINITION(level=ID){
      ETA_1 ~ Normal(mean=0, var=OMEGA_1 )
      ETA_2 ~ Normal(mean=0, var=OMEGA_2 )
   }# end RANDOM_VARIABLE_DEFINITION
	
   RANDOM_VARIABLE_DEFINITION(level=DV){
    EPS  ~ Normal(mean=0, var=1 )
   }# end RANDOM_VARIABLE_DEFINITION
   
  
	MODEL_PREDICTION {
	C = exp(ALPHA)*T/(exp(BETA)+T)
	}	
   
   
   OBSERVATION {
		Y : { type is proportionalError, proportional = CV, eps = EPS, prediction = C }	
 		   }
}
   
russu_Emax_task = taskObj{ 
	ESTIMATE {
      set algo is mcmc
        TARGET_SETTINGS(target="BUGS"){
            set 
            nchains = 1,
            burnin = 1,
            niter = 1,
        #    parameters = "LAMBDA0_POP,LAMBDA1_POP,K1_POP,K2_POP,W0_POP,LAMBDA0,K2,CV",
         #    thinupdater=1,
          #  thinsamples=5, # uno dei due
#           odesolver="LSODA", #oppure LSODA
        #    dic="true",
        	winbugsgui="true" # oppure no
        }
        }
}
	
russu_Emax_mog =mogObj {
	OBJECTS{
	russu_Emax_data : { type is dataObj }
	russu_pri : { type is priorObj }
	russu_Emax_mdl : { type is mdlObj }
	russu_Emax_task : { type is taskObj }
	}

}
