# Prior_Example334 : Parameters V, K are correlated (through IIV)
# 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 from MultivariateNormal1(mean=[0,0], covarianceMatrix=OMEGA_P)
#   - Additive residual error model
#########
## PRIOR SPECIFICATION:
## Population mean parameters are independent:
## ln(POP_K) ~ Normal1(ln(MU_POP_K), VAR_POP_K)
## ln(POP_V) ~ Normal1(ln(MU_POP_V), VAR_POP_V)
##
## Population between subject variabilities are correlated:
## OMEGA_P ~ Wishart2(R_mat, rho)
##
## 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
     
}

example334_pri = priorObj{
  PRIOR_PARAMETERS{
        # prior on "THETA"
		MU_POP_K = 0.1
		MU_POP_V = 8
        VAR_POP_K = 0.3
        VAR_POP_V = 0.3

        # prior on "OMEGA"
		a_OMEGA_K = 0.1
		b_OMEGA_K = 0.1
		a_OMEGA_V = 0.1
		b_OMEGA_V = 0.1
		a_OMEGA_T = 0.3
		b_OMEGA_T = 0.3
		
		# prior on "SIGMA"
        a_POP_T = 0.1
		b_POP_T = 0.1
 
        R_mat = [[1,0.1;0.1,1]]
        rho=2 
} # end PRIOR_PARAMETERS

   PRIOR_VARIABLE_DEFINITION{
	  # prior on "THETA"     
      lMU_POP_K = ln(MU_POP_K)
      lMU_POP_V = ln(MU_POP_V)		
	  lPOP_K ~ Normal(mean=lMU_POP_K, var=VAR_POP_K)
	  lPOP_V ~ Normal(mean=lMU_POP_V, var=VAR_POP_V)
	  POP_K = exp(lPOP_K)
	  POP_V = exp(lPOP_V)	
  
	  # prior on "OMEGA"	
	  invOMEGA_P ~ Wishart2(inverseScaleMatrix=R_mat, degreesOfFreedom=rho)
      OMEGA_P=inverse(invOMEGA_P)
      TAU_T ~ Gamma2(shape=a_OMEGA_T, rate=b_OMEGA_T)
      OMEGA_T=1/TAU_T 

      # prior on "SIGMA"  
      POP_T ~ Gamma2(shape=a_POP_T, rate=b_POP_T)            
    } # end PRIOR_VARIABLE_DEFINITION

} # end of parameter object

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

   STRUCTURAL_PARAMETERS{
      POP_K
      POP_V
      POP_T
   }# end STRUCTURAL_PARAMETERS

   VARIABILITY_PARAMETERS{
      invOMEGA_P :: matrix
      OMEGA_T
   }# end VARIABILITY_PARAMETERS


  RANDOM_VARIABLE_DEFINITION(level=ID){
      DIST_K_V ~ MultivariateNormal2(mean=[ln(POP_K),ln(POP_V)], precisionMatrix=invOMEGA_P) 
      eta_T ~ Normal( mean=0, var=OMEGA_T)
   }# end RANDOM_VARIABLE_DEFINITION

   INDIVIDUAL_VARIABLES{
	  K = exp(DIST_K_V[1]) 
	  V = exp(DIST_K_V[2]) 
      TAU : {type is linear, trans is ln, pop = POP_V, 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
	
example334 = mogObj{
	OBJECTS{
		example331_dat : { type is dataObj }
		example332_mdl : { type is mdlObj }
		example334_pri : { type is priorObj }
		example331_task : { type is taskObj }
	}
}
