$PROB negative binomial model characterizing epilepsy seizures frequency

$DATA epilepsy_NB_simulated.csv IGNORE=@
$INPUT 
ID ;number unique to each individual
OCC ;occasions (in months) 1=days 1 to 27, 2=days 28 to 55 etc
DAYS ;visit days based on start of trial (negative)
TIME ;time (in days) based on start of observation (starts at 0)
DV ;number of seizures per day
PDV ;previous count
AGE ;age (in yrs)
SEX ;1 = male 2= female
RACE ;0 to 6 

$PRED
;Overdispersion factor
 OVDP=THETA(1)*EXP(ETA(1))
 IF(PDV.EQ.1) OVDP=THETA(2)*EXP(ETA(2))

;Mean count
 LAMB=THETA(3)*EXP(ETA(3))
 IF(PDV.EQ.1) LAMB=THETA(4)*EXP(ETA(4))

;Components of NB distribution
 AGM1=DV+(1/OVDP)
 AGM2=1/OVDP

 LGAM1=LOG(SQRT(2*3.1415))+((AGM1)-0.5)*LOG((AGM1))-(AGM1)+LOG(1+1/(12*(AGM1)))
 LGAM2=LOG(SQRT(2*3.1415))+((AGM2)-0.5)*LOG((AGM2))-(AGM2)+LOG(1+1/(12*(AGM2)))

 LTRM1=(LOG(1/(1+OVDP*LAMB)))*(1/OVDP)
 LTRM2=(LOG(LAMB/(LAMB+1/OVDP)))*(DV)

;Stirling's formula approximating factorial
              LDVFAC=0
 IF (DV.GT.0) LDVFAC=LOG(SQRT(2*3.1415))+(DV+0.5)*LOG(DV)-DV+LOG(1+1/(12*DV))

;-2 log likelihood
 LYY=LGAM1-LDVFAC-LGAM2+LTRM1+LTRM2

 Y = -2*LYY

$THETA  (0,0.612)   ; OVDP pdv == 0
$THETA  (0,0.173)  ; OVDP pdv == 1
$THETA  (0,0.409)   ; Lambda pdv == 0
$THETA  (0,0.653)   ; Lambda pdv == 1
$OMEGA BLOCK(2)  2.4 ; IIV_OVDP_0
                 2.7 3.84     ; IIV_OVDP_1
$OMEGA BLOCK(2)  0.707  ; IIV_LAMBDA_0
                 0.636 0.884  ; IIV_LAMBDA_1

$ESTIMATION MAXEVAL=0 METHOD=COND LAPLACE -2LL
;$COV