;; 1. Initial estimation of all data
$PROBLEM published CV
$INPUT ID ID2 DV TRT TIME DAYS VISIT MDV1 MDV2 MDV3 MDV4 MDV5
$DATA diabet.csv IGNORE=@
$PRED
TVBAS= THETA(1)/10
PHBAS= LOG(TVBAS/(1-TVBAS))+ETA(1)
BAS = 10*EXP(PHBAS)/(1+EXP(PHBAS))
PEF = EXP(LOG(THETA(2)/(1- THETA(2)))+ETA(2))/(1+EXP(LOG(THETA(2)/(1- THETA(2)))+ETA(2)))
PHA = THETA(3)*EXP(ETA(3))
PIMAX = (1- PEF*(1-EXP(-(LOG(2)/PHA*DAYS))))
IF(TRT.EQ.1.AND.VISIT.GT.1) THEN
PEF = EXP(LOG(THETA(2)/(1- THETA(2)))+ETA(2)+THETA(4)+ETA(4))/(1+EXP(LOG(THETA(2)/(1- THETA(2)))+ETA(2)+THETA(4)+ETA(4)))
PHA = THETA(3)*EXP(ETA(3))
PIMAX = (1- PEF*(1-EXP(-(LOG(2)/PHA*DAYS))))
ENDIF
TVIPR = BAS*PIMAX/10
PHIPR = LOG(TVIPR/(1-TVIPR))
IPR = 11*EXP(PHIPR+EPS(1))/(1+EXP(PHIPR+EPS(1)))-0.5
Y = IPR
$THETA (4,6.25507,10) ; theta1
$THETA (0,0.225109,1) ; theta2
$THETA (0,26.8887) ; theta3
$THETA 0.140631 ; theta4
$OMEGA 0.506692 ; 1. OM BASELINE
$OMEGA 3.06564 ; 2. OM PLACEBO EFFECT
$OMEGA 0.478743 ; 3. OM PLACEBO HALF TIME
$OMEGA 0 FIX ; 4. OM Offset
$SIGMA 0.162231 ; variance prop res error, initial estimate
$ESTIMATION METHOD=1 INTER MAXEVAL=9999 PRINT=1 MSFO=msf1
$COVARIANCE
$TABLE ID TIME VISIT IPR PRED DV TRT NOAPPEND NOPRINT FILE=sdtab11 Estimate parameters from initial dataset
We bgin from an initial model and dataset, shown below.
Code
dat <- readr::read_csv(file.path('.','diabet.csv'),
show_col_types = FALSE)
knitr::kable(head(dat,n=21),booktabs=TRUE) |>
kableExtra::kable_styling("striped",full_width = F) | ID | ID2 | DV | TRT | TIME | DAYS | VISIT | MDV1 | MDV2 | MDV3 | MDV4 | MDV5 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 20804 | 20804 | 5 | 0 | 0.0714286 | 0.5 | 1 | 0 | 0 | 0 | 0 | 0 |
| 20804 | 20804 | 6 | 0 | 2.9285714 | 20.5 | 2 | 0 | 0 | 0 | 0 | 0 |
| 20804 | 20804 | 5 | 0 | 5.7857143 | 40.5 | 3 | 0 | 0 | 0 | 0 | 0 |
| 20804 | 20804 | 5 | 0 | 8.6428571 | 60.5 | 4 | 0 | 0 | 0 | 0 | 1 |
| 20804 | 20804 | 5 | 0 | 11.5000000 | 80.5 | 5 | 0 | 0 | 0 | 1 | 1 |
| 20804 | 20804 | 5 | 0 | 14.3571429 | 100.5 | 6 | 0 | 0 | 1 | 1 | 1 |
| 20804 | 20804 | 5 | 0 | 17.2142857 | 120.5 | 7 | 0 | 1 | 1 | 1 | 1 |
| 20903 | 20903 | 7 | 0 | 0.0714286 | 0.5 | 1 | 0 | 0 | 0 | 0 | 0 |
| 20903 | 20903 | 6 | 0 | 2.9285714 | 20.5 | 2 | 0 | 0 | 0 | 0 | 0 |
| 20903 | 20903 | 6 | 0 | 5.7857143 | 40.5 | 3 | 0 | 0 | 0 | 0 | 0 |
| 20903 | 20903 | 6 | 0 | 8.6428571 | 60.5 | 4 | 0 | 0 | 0 | 0 | 1 |
| 20903 | 20903 | 4 | 0 | 11.5000000 | 80.5 | 5 | 0 | 0 | 0 | 1 | 1 |
| 20903 | 20903 | 4 | 0 | 14.3571429 | 100.5 | 6 | 0 | 0 | 1 | 1 | 1 |
| 20903 | 20903 | 4 | 0 | 17.2142857 | 120.5 | 7 | 0 | 1 | 1 | 1 | 1 |
| 21001 | 21001 | 5 | 0 | 0.0714286 | 0.5 | 1 | 0 | 0 | 0 | 0 | 0 |
| 21001 | 21001 | 5 | 0 | 2.9285714 | 20.5 | 2 | 0 | 0 | 0 | 0 | 0 |
| 21001 | 21001 | 6 | 0 | 5.7857143 | 40.5 | 3 | 0 | 0 | 0 | 0 | 0 |
| 21001 | 21001 | 7 | 0 | 8.6428571 | 60.5 | 4 | 0 | 0 | 0 | 0 | 1 |
| 21001 | 21001 | 7 | 0 | 11.5000000 | 80.5 | 5 | 0 | 0 | 0 | 1 | 1 |
| 21001 | 21001 | 7 | 0 | 14.3571429 | 100.5 | 6 | 0 | 0 | 1 | 1 | 1 |
| 21001 | 21001 | 8 | 0 | 17.2142857 | 120.5 | 7 | 0 | 1 | 1 | 1 | 1 |
We estimate parameters:
PsN
execute run1.mod Examine the results:
PsN
sumo run1.lst-----------------------------------------------------------------------
run1.lst
Termination problems [ ERROR ]
Rounding errors [ ERROR ]
No zero gradients [ OK ]
No final zero gradients [ OK ]
Hessian not reset [ OK ]
No parameter near boundary [ OK ]
Covariance step [ ERROR ]
Total run time for model (hours:min:sec): 0:00:01
Estimation time for subproblem, sum over $EST (seconds): 0.29
Objective function value: 1253.3016
Number of observation records: 798
Number of individuals: 114
THETA OMEGA SIGMA
theta1 6.255 (........) 1. OM BASELINE 0.7118 (........) variance prop res error, initial estimate 0.4028 (........)
theta2 0.2251 (........) 2. OM PLACEBO EFFECT 1.751 (........)
theta3 26.89 (........) 3. OM PLACEBO HALF TIME 0.6919 (........)
theta4 0.1406 (........)
The relative standard errors for omega and sigma are reported on the approximate
standard deviation scale (SE/variance estimate)/2.
-----------------------------------------------------------------------We see that there are termination problems. If we look at the output file we see that the error is “rounding errors”. I get the same error when adding SIGDIGITS=8 in the $ESTIMATION line. This needs to be fixed before proceeding.
2 Initial simulation model
Can we get a successful minimization if we simulate data and then estimate from the same model?
;; 1. Initial estimation of all data
$PROBLEM published CV
$INPUT ID ID2 DV TRT TIME DAYS VISIT MDV1 MDV2 MDV3 MDV4 MDV5
$DATA diabet.csv IGNORE=@
$PRED
TVBAS= THETA(1)/10
PHBAS= LOG(TVBAS/(1-TVBAS))+ETA(1)
BAS = 10*EXP(PHBAS)/(1+EXP(PHBAS))
PEF = EXP(LOG(THETA(2)/(1- THETA(2)))+ETA(2))/(1+EXP(LOG(THETA(2)/(1- THETA(2)))+ETA(2)))
PHA = THETA(3)*EXP(ETA(3))
PIMAX = (1- PEF*(1-EXP(-(LOG(2)/PHA*DAYS))))
IF(TRT.EQ.1.AND.VISIT.GT.1) THEN
PEF = EXP(LOG(THETA(2)/(1- THETA(2)))+ETA(2)+THETA(4)+ETA(4))/(1+EXP(LOG(THETA(2)/(1- THETA(2)))+ETA(2)+THETA(4)+ETA(4)))
PHA = THETA(3)*EXP(ETA(3))
PIMAX = (1- PEF*(1-EXP(-(LOG(2)/PHA*DAYS))))
ENDIF
TVIPR = BAS*PIMAX/10
PHIPR = LOG(TVIPR/(1-TVIPR))
IPR = 11*EXP(PHIPR+EPS(1))/(1+EXP(PHIPR+EPS(1)))-0.5
Y = IPR
$THETA (4,6.25507,10) ; theta1
$THETA (0,0.225109,1) ; theta2
$THETA (0,26.8887) ; theta3
$THETA 0.140631 ; theta4
$OMEGA 0.506692 ; 1. OM BASELINE
$OMEGA 3.06564 ; 2. OM PLACEBO EFFECT
$OMEGA 0.478743 ; 3. OM PLACEBO HALF TIME
$OMEGA 0 FIX ; 4. OM Offset
$SIGMA 0.162231 ; variance prop res error, initial estimate
$SIMULATION (979879) ONLYSIM
$TABLE ID ID2 DV TRT TIME DAYS VISIT MDV1 MDV2 MDV3 MDV4 MDV5 NOAPPEND NOPRINT FILE=sim.datEstimation model: run2.mod ;; 1. Based on run1.mod
$PROBLEM published CV
$INPUT ID ID2 DV TRT TIME DAYS VISIT MDV1 MDV2 MDV3 MDV4 MDV5
$DATA sim.dat IGNORE=@
$PRED
TVBAS= THETA(1)/10
PHBAS= LOG(TVBAS/(1-TVBAS))+ETA(1)
BAS = 10*EXP(PHBAS)/(1+EXP(PHBAS))
PEF = EXP(LOG(THETA(2)/(1- THETA(2)))+ETA(2))/(1+EXP(LOG(THETA(2)/(1- THETA(2)))+ETA(2)))
PHA = THETA(3)*EXP(ETA(3))
PIMAX = (1- PEF*(1-EXP(-(LOG(2)/PHA*DAYS))))
IF(TRT.EQ.1.AND.VISIT.GT.1) THEN
PEF = EXP(LOG(THETA(2)/(1- THETA(2)))+ETA(2)+THETA(4)+ETA(4))/(1+EXP(LOG(THETA(2)/(1- THETA(2)))+ETA(2)+THETA(4)+ETA(4)))
PHA = THETA(3)*EXP(ETA(3))
PIMAX = (1- PEF*(1-EXP(-(LOG(2)/PHA*DAYS))))
ENDIF
TVIPR = BAS*PIMAX/10
PHIPR = LOG(TVIPR/(1-TVIPR))
IPR = 11*EXP(PHIPR+EPS(1))/(1+EXP(PHIPR+EPS(1)))-0.5
Y = IPR
$THETA (4,6.25507,10) ; theta1
$THETA (0,0.225109,1) ; theta2
$THETA (0,26.8887) ; theta3
$THETA 0.140631 ; theta4
$OMEGA 0.506692 ; 1. OM BASELINE
$OMEGA 3.06564 ; 2. OM PLACEBO EFFECT
$OMEGA 0.478743 ; 3. OM PLACEBO HALF TIME
$OMEGA 0 FIX ; 4. OM Offset
$SIGMA 0.162231 ; variance prop res error, initial estimate
$ESTIMATION METHOD=1 INTER MAXEVAL=9999 PRINT=1 MSFO=msf2
$COVARIANCE
$TABLE ID TIME VISIT IPR PRED DV TRT NOAPPEND NOPRINT FILE=sdtab2We simulate one dataset:
PsN
execute sim.mod Then we estimate using the same model:
PsN
execute run1.mod And look at the results:
PsN
sumo run2.lst-----------------------------------------------------------------------
run2.lst
Termination problems [ ERROR ]
No rounding errors [ OK ]
No zero gradients [ OK ]
No final zero gradients [ OK ]
Hessian not reset [ OK ]
No parameter near boundary [ OK ]
Covariance step [ ERROR ]
Total run time for model (hours:min:sec): 0:00:35
Estimation time for subproblem, sum over $EST (seconds): 34.12
Objective function value: 1184.9561
Number of observation records: 798
Number of individuals: 114
THETA OMEGA SIGMA
theta1 5.962 (........) 1. OM BASELINE 0.6756 (........) variance prop res error, initial estimate 0.3871 (........)
theta2 0.1707 (........) 2. OM PLACEBO EFFECT 1.296 (........)
theta3 27.89 (........) 3. OM PLACEBO HALF TIME 0.1146 (........)
theta4 0.8413 (........)
The relative standard errors for omega and sigma are reported on the approximate
standard deviation scale (SE/variance estimate)/2.
-----------------------------------------------------------------------Here again we have termination problems. This time the error is “DUE TO MAX. NO. OF FUNCTION EVALUATIONS EXCEEDED”.
This needs to be fixed before proceeding.