# MAGNI 2000
# The model describes Model#4 of the original publication.
# This model is used for simulation purposes of 7 subjects included in the dataset.
#
# Limitations:
# - The current implementation does not include distributions on the model parameters as required by the Bayesian context
# - Categorical covariates (HSTATUS, FEMALE) have been encoded as continuous covariates
#
# Compatibility:
# NONMEM
# (Monolix cannot handle covariate-dependent population parameters)
#
# Dataset:
# Constructed from Table 3 of the paper, with DVs simulated via original Matlab code provided by the authors in stochastic framework
#
# Evaluation method for this model:
# use the automatically generated NONMEM code for simulation;
# check the consistency of predictions with the DVs provided by the authors as predictions in stochastic framework
#
# By Elisa Borella, Letizia Carrara, Silvia Lavezzi, Paolo Magni, Giulia Muzio, Lorenzo Pasotti, Elena Tosca
# BMS lab - Univ Pavia
# Copyright 2015
#


magni2000_data = dataObj{

DECLARED_VARIABLES{Y1::observation Y2::observation Y3::observation Y4::observation}

   DATA_INPUT_VARIABLES{
      ID:{use is id}
      TIME:{use is idv} #ALWAYS ZERO
      WEIGHT : { use is covariate }
      HEIGHT : { use is covariate }
      #HSTATUS : { use is catCov withCategories {normal when 1, obese when 2, diabetic when 3} } 
      HSTATUS : { use is covariate }
      AGE : { use is covariate }
      #FEMALE : { use is catCov withCategories {yes when 1, no when 0}}
      FEMALE : { use is covariate}
      DV : { use is dv, define={1 in DVID as Y1, 2 in DVID as Y2, 3 in DVID as Y3, 4 in DVID as Y4} }
      DVID : { use is dvid }
   }
   
   SOURCE{
       srcfile: {file = "Simulated_magni2000_subjects.csv", 
            	  inputFormat is nonmemFormat#, 
            	  #ignore="@"
            	}
   }
}


magni2000_pri = priorObj{
	
	PRIOR_PARAMETERS{

	  	MU_mtsn = 5
		MU_mtso = 5
		MU_mtsd = 5
		MU_mFn = 1
		MU_mFo = 1
		MU_mFd = 1
		MU_atl = 30
		MU_btl = 1
		MU_aVm = 1
		MU_bVm = 1
		MU_aVf = 1
		MU_bVf = 1
		
		VAR_mtsn = 5^2
		VAR_mtso = 5^2
		VAR_mtsd = 5^2
		VAR_mFn = 1
		VAR_mFo = 1
		VAR_mFd = 1
		VAR_atl = 30^2
		VAR_btl = 1
		VAR_aVm = 1
		VAR_bVm = 1
		VAR_aVf = 1
		VAR_bVf = 1
		
		R = inverse([[0.5, 0, 0, 0;
		0, 0.1, 0, 0;
		0, 0, 3, 0;
		0, 0, 0, 0.4]])
		rho = 10
	}
	
	PRIOR_VARIABLE_DEFINITION{
					
		mtsn ~ Normal(mean=MU_mtsn, var=VAR_mtsn)
		mtso ~ Normal(mean=MU_mtso, var=VAR_mtso)
		mtsd ~ Normal(mean=MU_mtsd, var=VAR_mtsd)
		mFn ~ Normal(mean=MU_mFn, var=VAR_mFn)
		mFo ~ Normal(mean=MU_mFo, var=VAR_mFo)
		mFd ~ Normal(mean=MU_mFd, var=VAR_mFd)
		atl ~ Normal(mean=MU_atl, var=VAR_atl)
		btl ~ Normal(mean=MU_btl, var=VAR_btl)
		aVm ~ Normal(mean=MU_aVm, var=VAR_aVm)
		bVm ~ Normal(mean=MU_bVm, var=VAR_bVm)
		aVf ~ Normal(mean=MU_aVf, var=VAR_aVf)
		bVf ~ Normal(mean=MU_bVf, var=VAR_bVf)
		
		invOMEGA_PAR ~ Wishart1(scaleMatrix=R, degreesOfFreedom=rho)
		MATRIX_1 = inverse(invOMEGA_PAR)
		
		ADD_SMALL=1E-20
	}
}



magni2000_par = parObj{

   STRUCTURAL{
   
	mtsn : {value=4.99}
	mtso : {value=4.50}
	mtsd : {value=4.70}
	mFn : {value=0.767}
	mFo : {value=0.781}
	mFd : {value=0.778}
	atl : {value=26.70}
	btl : {value=0.209}
	aVm : {value=0.417}
	bVm : {value=0.89}
	aVf : {value=2.04}
	bVf : {value=1.77}
   	}
   
   VARIABILITY{

	MATRIX_1 : { matrixValue = triangle([1.3,
										0.0062, 0.0022,
										3.24, 0.071, 32.94, 
										0.59, -0.0055, 1.89, 0.71], 4, true)}
	ADD_SMALL : {value=1E-20, fix=true}
	}

}

magni2000_mdl = mdlObj{
 	
 	IDV {T}
 	
 	COVARIATES{
		#HSTATUS withCategories {normal, obese, diabetic}
		#FEMALE withCategories {yes, no}
		HSTATUS
		FEMALE
		AGE
		WEIGHT
		HEIGHT
	}
 	
 	VARIABILITY_LEVELS{
   	DV : { level=1, type is observation }
   	ID : { level=2, type is parameter }
   }

	
	STRUCTURAL_PARAMETERS{
	mtsn
	mtso
	mtsd
	mFn
	mFo
	mFd
	atl
	btl
	aVm
	bVm
	aVf
	bVf
	}
	
	VARIABILITY_PARAMETERS{
	MATRIX_1::matrix
	ADD_SMALL
	}
		
   	
   	GROUP_VARIABLES {
   	
	   	BSA = (WEIGHT^0.425)*(HEIGHT^0.725)*0.007184
   	
   		#GROUP_ts = if (HSTATUS==HSTATUS.normal) then mtsn
   		#	 	   else if (HSTATUS==HSTATUS.obese) then mtso
   		#	 	   else mtsd;
   			 
   		#GROUP_F = if (HSTATUS==HSTATUS.normal) then mFn
   		#		  else if (HSTATUS==HSTATUS.obese) then mFo
   		#		  else mFd;
   			
   		GROUP_tl = atl+btl*AGE
   		
   		#GROUP_V = if(FEMALE==FEMALE.no) then aVm+bVm*BSA
   		#		  else aVf+bVf*BSA
   				  
 		GROUP_ts = if (HSTATUS==1) then mtsn
   			 	   elseif (HSTATUS==2) then mtso
   			 	   elseif (HSTATUS==3) then mtsd
   			 
   		GROUP_F = if (HSTATUS==1) then mFn
   				  elseif (HSTATUS==2) then mFo
   				  elseif (HSTATUS==3) then mFd
   			
   		GROUP_V = if(FEMALE==0) then aVm+bVm*BSA
   				  else aVf+bVf*BSA
   	}
   	
   	INDIVIDUAL_VARIABLES{    	
    	ts_IND = GROUP_ts + eta_4[1]
		F_IND = GROUP_F + eta_4[2]
		tl_IND = GROUP_tl + eta_4[3]
		V_IND = GROUP_V + eta_4[4]
   	}
   	
   RANDOM_VARIABLE_DEFINITION(level=ID){
      eta_4 ~ MultivariateNormal1(mean=[0,0,0,0], covarianceMatrix=MATRIX_1)
   }# end RANDOM_VARIABLE_DEFINITION
	
   RANDOM_VARIABLE_DEFINITION(level=DV){
	  EPS  ~ Normal(mean=0, var=1 )
   }# end RANDOM_VARIABLE_DEFINITION
   
  
	MODEL_PREDICTION {
		ts_PRED = ts_IND
		F_PRED = F_IND
		tl_PRED = tl_IND
		V_PRED = V_IND
	}	
   
   
   OBSERVATION {	
		Y1 : { type is additiveError, additive = ADD_SMALL, eps = EPS, prediction = ts_PRED }
		Y2 : { type is additiveError, additive = ADD_SMALL, eps = EPS, prediction = F_PRED }
		Y3 : { type is additiveError, additive = ADD_SMALL, eps = EPS, prediction = tl_PRED }
		Y4 : { type is additiveError, additive = ADD_SMALL, eps = EPS, prediction = V_PRED }
 	}
}
   
magni2000_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
        }
	}
}
	
magni2000_mog =mogObj {
	OBJECTS{
	magni2000_data : { type is dataObj }
	magni2000_mdl : { type is mdlObj }
	magni2000_pri : { type is priorObj }
	magni2000_task : { type is taskObj }
	}
	
}
