# Prior_Example3321 : Parameters V, k and Vpop, kpop are correlated, compact notation 
# Dosing regimen: single IV bolus administration (at time 0)
# Structural model : 1 compartment model (V, K)
# 1 cmpt (V) + linear elimination (K)
# C(t) = Dose / V*exp(-K*t)
# Variability model: 
#   - IIV on V, K.
##    ETA_K ~ Normal(mean, stdev); ETA_V ~ Normal(mean, stdev)
#   - Additive residual error model
#########
## PRIOR SPECIFICATION:
## Population mean parameters are correlated
## lPOP_P ~ MultivariateNormal1(mean=lMU_POP_P, covarianceMatrix=SIGMA_POP_P)
##
## Population between subject variabilities are correlated:
## OMEGA_P ~ InverseWishart1(scaleMatrix=R, degreesOfFreedom=rho) # dof = degrees of freedom
##
## Residual error
## TAU ~ Gamma1(a_POP_T, b_POP_T)

example331_dat = dataObj {
   DECLARED_VARIABLES{ DOSE::dosingTarget Y::observation}
   DATA_INPUT_VARIABLES{
      ID:{use is id}
      AMT:{use is amt, variable=DOSE}
      TIME:{use is idv}
      DV:{use is dv, variable=Y}
      MDV : {use is mdv}
   }# end DATA_INPUT_VARIABLES

   SOURCE{
   mysrcfile:{file="TestData.csv",
      inputFormat is nonmemFormat}
   }# end SOURCE
     
}

example3321_pri = priorObj{

  PRIOR_PARAMETERS{
		# prior on "THETA"
		MU_POP_K = 0.1
		MU_POP_V = 8
		
        SIGMA_POP_P = [[1,0.1;0.1,1]] 
 		
		# prior on "OMEGA"
        a_OMEGA_T = 0.3
		b_OMEGA_T = 0.3
        
        R = [[1,0.1;0.1,1]]
        rho = 2 
        
        # prior on "SIGMA"
        a_POP_T = 0.1
		b_POP_T = 0.1
   } # end PRIOR

   PRIOR_VARIABLE_DEFINITION{

		# prior on "THETA"
	    lMU_POP_K = ln(MU_POP_K)
	    lMU_POP_V = ln(MU_POP_V)
	    
	    lMU_POP_P = [lMU_POP_K, lMU_POP_V]    	       
	 	lPOP_P ~ MultivariateNormal1(mean=lMU_POP_P, covarianceMatrix=SIGMA_POP_P)
	 	
	    POP_K = exp(lPOP_P[1])
	    POP_V = exp(lPOP_P[2]) 
	    
	    # prior on "OMEGA"
	 	OMEGA_P ~ InverseWishart1(scaleMatrix=R, degreesOfFreedom=rho) # dof = degrees of freedom
        OMEGA_T ~ InverseGamma1(shape=a_OMEGA_T, scale=b_OMEGA_T)
        
        # prior on "SIGMA"  
      	POP_T ~ Gamma2(shape=a_POP_T, rate=b_POP_T)     
    } # end PRIOR_VARIABLE_DEFINITION

} # end of parameter object

example3321_mdl = mdlObj{

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

   STRUCTURAL_PARAMETERS{
      lPOP_P :: vector
	  POP_T
   }# end STRUCTURAL_PARAMETERS

   VARIABILITY_PARAMETERS{
      OMEGA_P :: matrix
      OMEGA_T
   }# end VARIABILITY_PARAMETERS


   RANDOM_VARIABLE_DEFINITION(level=ID){
      lP ~ MultivariateNormal1(mean=lPOP_P, covarianceMatrix=OMEGA_P) 
      eta_T ~ Normal( mean=0, var=OMEGA_T)
      }# end RANDOM_VARIABLE_DEFINITION

   INDIVIDUAL_VARIABLES{
	  K = exp(lP[1])
	  V = exp(lP[2])
      TAU : {type is linear, trans is ln, pop = POP_T, ranEff=eta_T}
      SD_ADD = 1/sqrt(TAU)  
   }# end INDIVIDUAL_VARIABLES

   MODEL_PREDICTION{ 
      DOSE :: dosingVar
      C=DOSE/V*exp(-K*T)
   }# end MODEL_PREDICTION   
   
   
   RANDOM_VARIABLE_DEFINITION(level=DV){
   		eps ~ Normal( mean=0, var=1)
   	}
   
   OBSERVATION{
      Y : {type is additiveError, additive=SD_ADD, eps=eps, prediction= C}
   }# end OBSERVATION  

} # end of model object


example331_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
        }
}
} # end of task object
	
example331 = mogObj{
	OBJECTS{
		example331_dat : { type is dataObj }
		example3321_mdl : { type is mdlObj }
		example3321_pri : { type is priorObj }
		example331_task : { type is taskObj }
	}
}

