Single Dose Oral - PK
Vijay Ivaturi
2017-07-12
## Warning: Installed Rcpp (0.12.11.1) different from Rcpp used to build dplyr (0.12.11.3).
## Please reinstall dplyr to avoid random crashes or undefined behavior.Description
Model: One compartment, oral, linear absorption
A dataset containing simulated PK data. A single oral dose was given and rich pharmacokinetic measurements were taken. Demographic data of Gender, Weight, and Ethinicity were collected.
Data Structure
id- Numerical ID (1-40)time- Time in hours (0-19)mdv- missing dependent variable (1 when concentration is missing and 0 otherwise)evid- event identifier (1 = dose event; 0 = observation)dv- Plasma concentration in ug/Lamt- dose of drug in micrograms- Covariates:
sex: Sex (binary: 0/1) - 0=Male, 1=Femalewt: Weight in kg (continuous)etn: Ethnicity (categorical: 0/1/2) 0=Caucasian, 1=Black, 2=Hispanic
Model Structure
Notes -
- \(Dose\) (units - microgram) is given at time \(t = 0\) directly into the \(Depot\) compartment
- \(Ka\) is the first-order rate constant for absorption of the drug from the \(Depot\) to the \(Central\) compartment
- \(Central\) compartment is the observation compartment where concentrations of the drug \(C_p(t)\) are measured (units - microgram/Liter).
- Volume of distribution parameter \(V\) (units = Liter) is a proportionality constant that scales the amount of drug in the \(Central\) compartment (mg) to the observed concentration (mg/L) - \(Concentration = \frac{Dose}{V}\)
- \(CL\) (units - Liter/time) is the parameter that relates the rate of elimination (units - microgram/time) to the measured concentration. \[ Rate of Elimination = CL \times C_p \]
The population values of \(CL\), \(TVCL\) is allometrically scaled by body weight wt and dependent on sex as follows: \[ CL = TVCL \times (\frac{WEIGHT}{70})^{0.75} \times {\theta}^{SEX} \]
- \(Out\) is collection compartment without any exit rates. Such a compartment can be use to accumulate drug over time and is commonly referred to as area under the curve or \(AUC\).
\[ AUC = \int\limits_0^\infty C_p\, \mathrm{d}t\]
ODE
\[ \begin{array}{lcl} \frac{ dDepot }{ dt } & = &- Ka\times Depot \\ \frac{ dCentral }{ dt } & = &Ka\times Depot- \frac{CL}{V}\times Central \\\end{array} \]
Exploratory plots
data1 %>% filter(amt==0) %>%
mutate(sex =ifelse(sex==0,"Male","Female")) %>%
ggplot(aes(x=time, y=dv, group=id))+
geom_line()+
scale_y_log10()+
base_theme()+
facet_wrap(~sex)+
labs(title = "Individual subject concentration-time profiles",
x= "Time after dose (hrs)",
y = "Concentration (ug/L)")
Base Model - run100.mod (without covariates)
Analytical Solution Based
In the model file below, the analytical solution is invoked by the line $SUBROUTINE ADVAN2 TRANS2 where ADVAN2 deals with models of one compartment oral absorption with a closed form of the following structure
\[ Cp = \frac{F\times Dose \times Ka}{V \times (Ka-kel)}\times[{e}^{-kel \times t} \times {e}^{-ka \times t} ] \] where F stands for bioavailability. In this example, F is assumed to be 1.
\(CL\) can also be derived from \(k\), \(V\) where, \[ CL = kel \times V\]
$EST line provides the estimation method, where METHOD=1 stands for FOCE and INTERACTION makes the method equivalent to FOCEI
cat(readLines("../modeling/run100.mod"), sep = '\n')Warning in readLines("../modeling/run100.mod"): incomplete final line found
on '../modeling/run100.mod';; 2. Description: PK model 1 cmt base
;; x1. Author: Vijay
$PROBLEM PK model 1 cmt base
$INPUT ID TIME MDV EVID DV AMT SEX WT ETN
$DATA ../data/data1.csv IGNORE=@
$SUBROUTINES ADVAN2 TRANS2
$PK
ET=1
IF(ETN.EQ.3) ET=1.3
KA = THETA(1)
CL = THETA(2)*EXP(ETA(1))
V = THETA(3)*EXP(ETA(2))
SC=V
$THETA
(0, 2) ; KA
(0, 20) ; CL
(0, 100) ; V
$OMEGA
0.05 ; iiv CL
0.2 ; iiv V
$ERROR
IPRED = F
IRES = DV-IPRED
PROP=SQRT(SIGMA(1,1))*IPRED
ADD=SQRT(SIGMA(2,2))
; COV = SIGMA(1,2); if covariance is required between two epsilons
; SD=SQRT(PROP*PROP + ADD*ADD + 2*COV) ; if covariance is required between two epsilons
SD=SQRT(PROP*PROP + ADD*ADD)
IWRES = IRES/SD
Y = IPRED+IPRED*EPS(1) + EPS(2)
$SIGMA
0.02 ; proportional error
0.1 ; additive error
$EST METHOD=1 INTERACTION MAXEVAL=9999 SIG=3 PRINT=5 NOABORT POSTHOC MSFO=msf100
$COV PRINT=E
;$SIM (1234) NSUBPROBLEMS=1 ONLYSIM
$TABLE ID TIME DV MDV EVID IPRED IWRES CWRES ONEHEADER NOPRINT FILE=sdtab100
$TABLE ID CL V ETA1 ETA2 ONEHEADER NOPRINT FILE=patab100
$TABLE ID SEX ETN ONEHEADER NOPRINT FILE=catab100
$TABLE ID WT ONEHEADER NOPRINT FILE=cotab100
$TABLE ID CL V SEX ETN WT ETA1 ETA2 ONEHEADER NOPRINT FILE=mytab100The results of fitting a one compartmental oral model with FOCEI are given below. The output provides information regarding the general model fit such as, whether the estimation algorithm successfully minimized the objective function value, were ther any zero gradients etc. The estimation run time and the correspondind covariance step (to calculate the standard errors) time are provided. This is followed by the two most important blocks, the Objective Function value of the likelihood function and the last table of the output which provides the estimates of the various fixed and random effect parameters.
cd ../modeling/.
sumo run100.lst
cd ../lab-notebook## -----------------------------------------------------------------------
##
## run100.lst
##
## Successful minimization [ OK ]
## No rounding errors [ OK ]
## No zero gradients [ OK ]
## No final zero gradients [ OK ]
## Hessian not reset [ OK ]
## No parameter near boundary [ OK ]
## Covariance step [ OK ]
## No large relative standard errors found [ OK ]
## Condition number [ OK ]
## Correlations [ OK ]
##
## Total run time for model (hours:min:sec): 0:00:04
## Estimation time for subproblem, sum over $EST (seconds): 1.65
## Covariance time for subproblem, sum over $EST (seconds): 1.05
##
## Objective function value: 2097.4206
##
## Condition number: 11.91
##
## ETA shrinkage: 3.7615E+00 1.0000E-10
## EPS shrinkage: 3.8821E+00 5.0875E+00
##
## Number of observation records: 760
## Number of individuals: 40
##
## THETA OMEGA SIGMA
## KA 2.272 (0.03566) iiv CL 0.3691 (0.08386) proportional error 0.11 (0.08745)
## CL 42.18 (0.06163) iiv V 0.3799 (0.08066) additive error 1.32 (0.07523)
## V 469.2 (0.06032)
##
## The relative standard errors for omega and sigma are reported on the approximate
## standard deviation scale (SE/variance estimate)/2.
## -----------------------------------------------------------------------ODE Solution Based
Below is the representation of the same model above in a differential equation setting where $ODE specifies the differential equations. The only other change to the model file above is the specification of ADVAN6 which takes ODE equations and the tolerance, TOL is set to 6 (10-6) and $MODEL which specifies the number of compartments and the names of the compartments.
cat(readLines("../modeling/run100_ode.mod"), sep = '\n')Warning in readLines("../modeling/run100_ode.mod"): incomplete final line
found on '../modeling/run100_ode.mod';; 1. Based on: 100
;; 2. Description: PK model 1 cmt base
;; x1. Author: Vijay
$PROBLEM PK model 1 cmt base
$INPUT ID TIME MDV EVID DV AMT SEX WT ETN
$DATA ../data/data1.csv IGNORE=@
$SUBROUTINE ADVAN6 TRANS=1 TOL=6
$MODEL NCOMP=2
COMP=(GUT)
COMP=(CENTRAL)
$PK
ET=1
IF(ETN.EQ.3) ET=1.3
KA = THETA(1)
CL = THETA(2)*EXP(ETA(1))
V = THETA(3)*EXP(ETA(2))
SC=V
$DES
C=A(2)/V
DADT(1)=-KA*A(1)
DADT(2)=KA*A(1)-CL*C
$THETA
(0, 2) ; KA
(0, 20) ; CL
(0, 100) ; V
$OMEGA
0.05 ; iiv CL
0.2 ; iiv V
$ERROR
IPRED = A(2)/V
IRES = DV-IPRED
PROP=SQRT(SIGMA(1,1))*IPRED
ADD=SQRT(SIGMA(2,2))
; COV = SIGMA(1,2); if covariance is required between two epsilons
; SD=SQRT(PROP*PROP + ADD*ADD + 2*COV) ; if covariance is required between two epsilons
SD=SQRT(PROP*PROP + ADD*ADD)
IWRES = IRES/SD
Y = IPRED+IPRED*EPS(1) + EPS(2)
$SIGMA
0.02 ; proportional error
0.1 ; additive error
$EST METHOD=1 INTERACTION MAXEVAL=9999 SIG=3 PRINT=5 NOABORT POSTHOC MSFO=msf100
$COV PRINT=E
;$SIM (1234) NSUBPROBLEMS=1 ONLYSIM
$TABLE ID TIME DV MDV EVID IPRED IWRES CWRES ONEHEADER NOPRINT FILE=sdtab100_ode
$TABLE ID CL V ETA1 ETA2 ONEHEADER NOPRINT FILE=patab100_ode
$TABLE ID SEX ETN ONEHEADER NOPRINT FILE=catab100_ode
$TABLE ID WT ONEHEADER NOPRINT FILE=cotab100_ode
$TABLE ID CL V SEX ETN WT ETA1 ETA2 ONEHEADER NOPRINT FILE=mytab100_odeAs you can see below, the results by the ODE are identical to the analytical solution
cd ../modeling/.
sumo run100_ode.lst
cd ../lab-notebook## -----------------------------------------------------------------------
##
## run100_ode.lst
##
## Successful minimization [ OK ]
## No rounding errors [ OK ]
## No zero gradients [ OK ]
## No final zero gradients [ OK ]
## Hessian not reset [ OK ]
## No parameter near boundary [ OK ]
## Covariance step [ OK ]
## No large relative standard errors found [ OK ]
## Condition number [ OK ]
## Correlations [ OK ]
##
## Total run time for model (hours:min:sec): 0:00:20
## Estimation time for subproblem, sum over $EST (seconds): 10.53
## Covariance time for subproblem, sum over $EST (seconds): 7.13
##
## Objective function value: 2097.4206
## Delta ofv (from run100): 0.0000
##
## Condition number: 11.85
##
## ETA shrinkage: 3.7615E+00 1.0000E-10
## EPS shrinkage: 3.8821E+00 5.0875E+00
##
## Number of observation records: 760
## Number of individuals: 40
##
## THETA OMEGA SIGMA
## KA 2.272 (0.03566) iiv CL 0.3691 (0.08386) proportional error 0.11 (0.08721)
## CL 42.18 (0.06163) iiv V 0.3799 (0.08064) additive error 1.32 (0.07509)
## V 469.2 (0.06032)
##
## The relative standard errors for omega and sigma are reported on the approximate
## standard deviation scale (SE/variance estimate)/2.
## -----------------------------------------------------------------------Final Model - run105.mod (with covariates)
Analytical Solution Based
The model specified below is the final model where the \(CL\) is dependent of \(WT\) and \(SEX\).
cat(readLines("../modeling/run105.mod"), sep = '\n');; 1. Based on: 102
;; 2. Description: PK model 1 cmt base WT-CL allom SEX-CL
;; x1. Author: Vijay
$PROBLEM PK model 1 cmt base
$INPUT ID TIME MDV EVID DV AMT SEX WT ETN
$DATA ../data/data1.csv IGNORE=@
$SUBROUTINES ADVAN2 TRANS2
$PK
ET=1
IF(ETN.EQ.3) ET=1.3
KA = THETA(1)
CL = THETA(2)*((WT/70)**0.75)*( THETA(4)**SEX) * EXP(ETA(1))
V = THETA(3)*EXP(ETA(2))
SC=V
$THETA
(0, 2) ; KA
(0, 3) ; CL
(0, 10) ; V2
(0,1) ; SEX - CL
$OMEGA
0.05 ; iiv CL
0.2 ; iiv V2
$ERROR
IPRED = F
IRES = DV-IPRED
PROP=SQRT(SIGMA(1,1))*IPRED
ADD=SQRT(SIGMA(2,2))
; COV = SIGMA(1,2); if covariance is required between two epsilons
; SD=SQRT(PROP*PROP + ADD*ADD + 2*COV) ; if covariance is required between two epsilons
SD=SQRT(PROP*PROP + ADD*ADD)
IWRES = IRES/SD
Y = IPRED+IPRED*EPS(1) + EPS(2)
$SIGMA
0.02 ; proportional error
1 ; additive error
$EST METHOD=1 INTERACTION MAXEVAL=9999 SIG=3 PRINT=5 NOABORT POSTHOC MSFO=msf105
$COV PRINT=E
;$SIM (1234) NSUBPROBLEMS=1 ONLYSIM
$TABLE ID TIME DV MDV EVID IPRED IWRES ONEHEADER NOPRINT FILE=sdtab105
$TABLE ID CL V ETA1 ETA2 ONEHEADER NOPRINT FILE=patab105
$TABLE ID SEX ETN ONEHEADER NOPRINT FILE=catab105
$TABLE ID WT ONEHEADER NOPRINT FILE=cotab105
$TABLE ID CL V SEX ETN WT ETA1 ETA2 ONEHEADER NOPRINT FILE=mytab105Below are the results of the analytical solution for the final model
cd ../modeling/.
sumo run105.lst
cd ../lab-notebook## -----------------------------------------------------------------------
##
## run105.lst
##
## Successful minimization [ OK ]
## No rounding errors [ OK ]
## No zero gradients [ OK ]
## No final zero gradients [ OK ]
## Hessian not reset [ OK ]
## No parameter near boundary [ OK ]
## Covariance step [ OK ]
## No large relative standard errors found [ OK ]
## Condition number [ OK ]
## Correlations [ OK ]
##
## Total run time for model (hours:min:sec): 0:00:05
## Estimation time for subproblem, sum over $EST (seconds): 1.96
## Covariance time for subproblem, sum over $EST (seconds): 1.4
##
## Objective function value: 2053.8304
## Delta ofv (from run102): -35.0958
##
## Condition number: 14.39
##
## EPS shrinkage: 3.4665E+00 4.6819E+00
##
## Number of observation records: 760
## Number of individuals: 40
##
## THETA OMEGA SIGMA
## KA 2.268 (0.03573) iiv CL 0.1911 (0.1516) proportional error 0.1105 (0.08668)
## CL 74.17 (0.03831) iiv V2 0.3791 (0.08068) additive error 1.317 (0.07509)
## V2 468.6 (0.06035)
## SEX - CL 0.5876 (0.06594)
##
## The relative standard errors for omega and sigma are reported on the approximate
## standard deviation scale (SE/variance estimate)/2.
## -----------------------------------------------------------------------ODE Solution Based
The same final model as above but written in the ODE form
cat(readLines("../modeling/run105_ode.mod"), sep = '\n');; 1. Based on: 102
;; 2. Description: PK model 1 cmt base WT-CL allom SEX-CL
;; x1. Author: Vijay
$PROBLEM PK model 1 cmt base
$INPUT ID TIME MDV EVID DV AMT SEX WT ETN
$DATA ../data/data1.csv IGNORE=@
$SUBROUTINE ADVAN6 TRANS=1 TOL=6
$MODEL NCOMP=2
COMP=(GUT)
COMP=(CENTRAL)
$PK
ET=1
IF(ETN.EQ.3) ET=1.3
KA = THETA(1)
CL = THETA(2)*((WT/70)**0.75)*( THETA(4)**SEX) * EXP(ETA(1))
V = THETA(3)*EXP(ETA(2))
SC=V
$DES
C=A(2)/V
DADT(1)=-KA*A(1)
DADT(2)=KA*A(1)-CL*C
$THETA
(0, 2) ; KA
(0, 3) ; CL
(0, 10) ; V2
(0,1) ; SEX - CL
$OMEGA
0.05 ; iiv CL
0.2 ; iiv V2
$ERROR
IPRED = A(2)/V
IRES = DV-IPRED
PROP=SQRT(SIGMA(1,1))*IPRED
ADD=SQRT(SIGMA(2,2))
; COV = SIGMA(1,2); if covariance is required between two epsilons
; SD=SQRT(PROP*PROP + ADD*ADD + 2*COV) ; if covariance is required between two epsilons
SD=SQRT(PROP*PROP + ADD*ADD)
IWRES = IRES/SD
Y = IPRED+IPRED*EPS(1) + EPS(2)
$SIGMA
0.02 ; proportional error
1 ; additive error
$EST METHOD=1 INTERACTION MAXEVAL=9999 SIG=3 PRINT=5 NOABORT POSTHOC MSFO=msf105
$COV PRINT=E
;$SIM (1234) NSUBPROBLEMS=1 ONLYSIM
$TABLE ID TIME DV MDV EVID IPRED IWRES ONEHEADER NOPRINT FILE=sdtab105_ode
$TABLE ID CL V ETA1 ETA2 ONEHEADER NOPRINT FILE=patab105_ode
$TABLE ID SEX ETN ONEHEADER NOPRINT FILE=catab105_ode
$TABLE ID WT ONEHEADER NOPRINT FILE=cotab105_ode
$TABLE ID CL V SEX ETN WT ETA1 ETA2 ONEHEADER NOPRINT FILE=mytab105_odeThe results of the ODE model are identical to the results from the analytical solution
cd ../modeling/.
sumo run105_ode.lst
cd ../lab-notebook## -----------------------------------------------------------------------
##
## run105_ode.lst
##
## Successful minimization [ OK ]
## No rounding errors [ OK ]
## No zero gradients [ OK ]
## No final zero gradients [ OK ]
## Hessian not reset [ OK ]
## No parameter near boundary [ OK ]
## Covariance step [ OK ]
## No large relative standard errors found [ OK ]
## Condition number [ OK ]
## Correlations [ OK ]
##
## Total run time for model (hours:min:sec): 0:00:24
## Estimation time for subproblem, sum over $EST (seconds): 13.03
## Covariance time for subproblem, sum over $EST (seconds): 9.44
##
## Objective function value: 2053.8304
## Delta ofv (from run102): -35.0958
##
## Condition number: 14.24
##
## EPS shrinkage: 3.4707E+00 4.6863E+00
##
## Number of observation records: 760
## Number of individuals: 40
##
## THETA OMEGA SIGMA
## KA 2.268 (0.03573) iiv CL 0.1911 (0.1512) proportional error 0.1105 (0.08655)
## CL 74.17 (0.0383) iiv V2 0.3791 (0.08068) additive error 1.317 (0.07486)
## V2 468.6 (0.06035)
## SEX - CL 0.5876 (0.06593)
##
## The relative standard errors for omega and sigma are reported on the approximate
## standard deviation scale (SE/variance estimate)/2.
## -----------------------------------------------------------------------Base Model evaluation
Some simple diagnostic plots required for base model evaluation (without covariates)
setwd("~/PKPDexamples/modeling/")
runno <- "100"
xpdb100 <- xpose.data(runno)##
## Looking for NONMEM table files.
## Reading sdtab100
## Reading patab100
## Reading catab100
## Reading cotab100
## Reading mytab100
## Table files read.
## Reading run100.phi
##
## Looking for NONMEM simulation table files.
## No simulated table files read.fit100 <- xpdb100@Data
setwd("~/PKPDexamples/lab-notebook/")fit100 %>%
filter(WRES!=0) %>%
ggplot(aes(x=PRED, y=DV))+
geom_point(shape=1,color="blue",size=3)+
stat_smooth(method="lm",se=FALSE,color="red",size=1.3)+
geom_abline(intercept=0,slope=1, color="black",size=1)+
base_theme(axis_title_x=16, axis_text_x = 14)+
labs(title ="PK Base model GOF",
caption = "source:run100.lst",
x="Population predictions (ng/mL)", y="Observed Plasma Concentrations (ng/mL)")+
theme(plot.title = element_text(size=12,face="bold"),
plot.caption = element_text(face="italic"))
Population Predictions vs Observations
fit100 %>%
filter(WRES!=0) %>%
ggplot(aes(x=IPRED, y=DV))+
geom_point(shape=1,color="blue",size=3)+
stat_smooth(method="lm",se=FALSE,color="red",size=1.3)+
geom_abline(intercept=0,slope=1, color="black",size=1)+
base_theme(axis_title_x=16, axis_text_x = 14)+
labs(title ="PK Base model GOF",
caption = "source:run100.lst",
x="Individual predictions (ng/mL)", y="Observed Plasma Concentrations (ng/mL)")+
theme(plot.title = element_text(size=12,face="bold"),
plot.caption = element_text(face="italic"))
Individual Predictions vs Observations
fit100 %>%
filter(WRES!=0) %>%
ggplot(aes(y=IWRES, x=TIME))+
geom_point(shape=1,color="blue",size=3)+
stat_smooth(method="loess",se=FALSE,color="red",size=1.3)+
geom_hline(yintercept=0, color="black",size=1)+
base_theme(axis_title_x=18, axis_text_x = 16)+
scale_x_continuous(breaks=seq(0,60,10))+
labs(title ="PK Base model GOF",
caption = "source:run100.lst",
y="Individual Weighted Residuals", x="Time after first dose (Days)")+
theme(plot.title = element_text(size=12,face="bold"),
plot.caption = element_text(face="italic"))
Individual Weighted Residuals vs Time
fit100 %>%
filter(WRES!=0) %>%
ggplot(aes(y=IWRES, x=IPRED))+
geom_point(shape=1,color="blue",size=3)+
stat_smooth(method="loess",se=FALSE,color="red",size=1.3)+
geom_hline(yintercept=0, color="black",size=1)+
base_theme(axis_title_x=16, axis_text_x = 14)+
labs(title ="PK Base model GOF",
caption = "source:run100.lst",
y="Individual Weighted Residuals", x="Individual Population Predictions (ng/mL)")+
theme(plot.title = element_text(size=12,face="bold"),
plot.caption = element_text(face="italic"))
Individual Weighted Residuals vs IPRED
EBE vs covariate plots for the PK Model
ebes <- as_data_frame(fit100) %>%
filter(!duplicated(ID)) %>%
mutate(SEX =ifelse(SEX==0,"Male","Female")) %>%
mutate(ETN =ifelse(ETN==0,"Caucasian",
ifelse(ETN == 1, "Black", "Hispanic"))) %>%
select(id=ID, starts_with("ETA"), WT, SEX, ETN)
ebelongcont <- ebes %>%
gather(key=eta, value=etaval, ETA1:ETA2, -id ) %>%
gather(key=cov,value=coval,WT,-id) %>%
arrange(eta) %>%
mutate(etaval=as.numeric(etaval))
ebelongcat <- ebes %>%
gather(key=eta, value=etaval, ETA1:ETA2, -id ) %>%
gather(key=cov,value=coval,SEX:ETN,-id) %>%
arrange(eta) %>%
mutate(etaval=as.numeric(etaval),
coval = factor(coval))cont_ebe_plots <- function(df,etaname){
p1 <-ggplot(data=df,aes(y=etaval,x=coval))+
geom_point()+
facet_wrap(~cov, scales = "free")+
# xlim(-0.5,0.5)+
stat_smooth(method = "lm")+
geom_hline(yintercept=0,col="red",linetype="dashed")+
base_theme()+
labs(title ="Base model - EBE versus covariate plots",
caption = "source:run100.lst",
y="ETA Value", x="Covariate Value")+
theme(plot.title = element_text(size=12,face="bold"),
plot.caption = element_text(face="italic"))
return(p1)
}
ebelongcont %>% split(.$eta) %>% lapply(cont_ebe_plots)## $ETA1
##
## $ETA2
cat_ebe_plots <- function(df,etaname){
p1 <-ggplot(data=df,aes(x=coval,y=etaval))+
geom_boxplot()+coord_flip()+
facet_wrap(~cov, scales="free_y")+
geom_hline(yintercept = 0, color="red", linetype="dashed")+
base_theme()+
labs(title ="PK Base model - EBE versus covariate plots",
caption = "source:run100.lst",
y="ETA Value", x="Covariate Value")+
theme(plot.title = element_text(size=12,face="bold"),
plot.caption = element_text(face="italic"))
# scale_x_continuous((breaks=c(-2,-1.5,-1,-0.5,0,0.5,1,1.5,2))
return(p1)
}
ebelongcat %>% split(.$eta) %>% lapply(cat_ebe_plots)## $ETA1
##
## $ETA2
Final Model evaluation
Some simple diagnostic plots required for base model evaluation (with covariates)
setwd("~/PKPDexamples/modeling/")
runno <- "105"
xpdb105 <- xpose.data(runno)##
## Looking for NONMEM table files.
## Reading sdtab105
## Reading patab105
## Reading catab105
## Reading cotab105
## Reading mytab105
## Table files read.
## Reading run105.phi
##
## Looking for NONMEM simulation table files.
## No simulated table files read.fit105 <- xpdb105@Data
setwd("~/PKPDexamples/lab-notebook/")fit105 %>%
filter(WRES!=0) %>%
ggplot(aes(x=PRED, y=DV))+
geom_point(shape=1,color="blue",size=3)+
stat_smooth(method="lm",se=FALSE,color="red",size=1.3)+
geom_abline(intercept=0,slope=1, color="black",size=1)+
base_theme(axis_title_x=16, axis_text_x = 14)+
labs(title ="PK Final model GOF",
caption = "source:run105.lst",
x="Population predictions (ng/mL)", y="Observed Plasma Concentrations (ng/mL)")+
theme(plot.title = element_text(size=12,face="bold"),
plot.caption = element_text(face="italic"))
Population Predictions vs Observations
fit105 %>%
filter(WRES!=0) %>%
ggplot(aes(x=IPRED, y=DV))+
geom_point(shape=1,color="blue",size=3)+
stat_smooth(method="lm",se=FALSE,color="red",size=1.3)+
geom_abline(intercept=0,slope=1, color="black",size=1)+
base_theme(axis_title_x=16, axis_text_x = 14)+
labs(title ="PK Final model GOF",
caption = "source:run105.lst",
x="Individual predictions (ng/mL)", y="Observed Plasma Concentrations (ng/mL)")+
theme(plot.title = element_text(size=12,face="bold"),
plot.caption = element_text(face="italic"))
Individual Predictions vs Observations
fit105 %>%
filter(WRES!=0) %>%
ggplot(aes(y=IWRES, x=TIME))+
geom_point(shape=1,color="blue",size=3)+
stat_smooth(method="loess",se=FALSE,color="red",size=1.3)+
geom_hline(yintercept=0, color="black",size=1)+
base_theme(axis_title_x=18, axis_text_x = 16)+
scale_x_continuous(breaks=seq(0,60,10))+
labs(title ="PK Final model GOF",
caption = "source:run105.lst",
y="Individual Weighted Residuals", x="Time after first dose (Days)")+
theme(plot.title = element_text(size=12,face="bold"),
plot.caption = element_text(face="italic"))
Individual Weighted Residuals vs Time
fit105 %>%
filter(WRES!=0) %>%
ggplot(aes(y=IWRES, x=IPRED))+
geom_point(shape=1,color="blue",size=3)+
stat_smooth(method="loess",se=FALSE,color="red",size=1.3)+
geom_hline(yintercept=0, color="black",size=1)+
base_theme(axis_title_x=16, axis_text_x = 14)+
labs(title ="PK Final model GOF",
caption = "source:run105.lst",
y="Individual Weighted Residuals", x="Individual Population Predictions (ng/mL)")+
theme(plot.title = element_text(size=12,face="bold"),
plot.caption = element_text(face="italic"))
Individual Weighted Residuals vs IPRED
EBE vs covariate plots for the PK Model
ebes <- as_data_frame(fit105) %>%
filter(!duplicated(ID)) %>%
mutate(SEX =ifelse(SEX==0,"Male","Female")) %>%
mutate(ETN =ifelse(ETN==0,"Caucasian",
ifelse(ETN == 1, "Black", "Hispanic"))) %>%
select(id=ID, starts_with("ETA"), WT, SEX, ETN)
ebelongcont <- ebes %>%
gather(key=eta, value=etaval, ETA1:ETA2, -id ) %>%
gather(key=cov,value=coval,WT,-id) %>%
arrange(eta) %>%
mutate(etaval=as.numeric(etaval))
ebelongcat <- ebes %>%
gather(key=eta, value=etaval, ETA1:ETA2, -id ) %>%
gather(key=cov,value=coval,SEX:ETN,-id) %>%
arrange(eta) %>%
mutate(etaval=as.numeric(etaval),
coval = factor(coval))cont_ebe_plots <- function(df,etaname){
p1 <-ggplot(data=df,aes(y=etaval,x=coval))+
geom_point()+
facet_wrap(~cov, scales = "free")+
# xlim(-0.5,0.5)+
stat_smooth(method = "lm")+
geom_hline(yintercept=0,col="red",linetype="dashed")+
base_theme()+
labs(title ="Final model - EBE versus covariate plots",
caption = "source:run105.lst",
y="ETA Value", x="Covariate Value")+
theme(plot.title = element_text(size=12,face="bold"),
plot.caption = element_text(face="italic"))
return(p1)
}
ebelongcont %>% split(.$eta) %>% lapply(cont_ebe_plots)## $ETA1
##
## $ETA2
cat_ebe_plots <- function(df,etaname){
p1 <-ggplot(data=df,aes(x=coval,y=etaval))+
geom_boxplot()+coord_flip()+
facet_wrap(~cov, scales="free_y")+
geom_hline(yintercept = 0, color="red", linetype="dashed")+
base_theme()+
labs(title ="PK Final model - EBE versus covariate plots",
caption = "source:run105.lst",
y="ETA Value", x="Covariate Value")+
theme(plot.title = element_text(size=12,face="bold"),
plot.caption = element_text(face="italic"))
# scale_x_continuous((breaks=c(-2,-1.5,-1,-0.5,0,0.5,1,1.5,2))
return(p1)
}
ebelongcat %>% split(.$eta) %>% lapply(cat_ebe_plots)## $ETA1
##
## $ETA2
Output of final model
cat(readLines("../modeling/run105.lst"), sep = '\n')Sun Jul 2 07:19:20 EDT 2017
;; 1. Based on: 102
;; 2. Description: PK model 1 cmt base WT-CL allom SEX-CL
;; x1. Author: Vijay
$PROBLEM PK model 1 cmt base
$INPUT ID TIME MDV EVID DV AMT SEX WT ETN
$DATA data1.csv IGNORE=@
$SUBROUTINE ADVAN2 TRANS2
$PK
ET=1
IF(ETN.EQ.3) ET=1.3
KA = THETA(1)
CL = THETA(2)*((WT/70)**0.75)*( THETA(4)**SEX) * EXP(ETA(1))
V = THETA(3)*EXP(ETA(2))
SC=V
$THETA (0,2) ; KA
(0,3) ; CL
(0,10) ; V2
(0,1) ; SEX - CL
$OMEGA 0.05 ; iiv CL
0.2 ; iiv V2
$ERROR
IPRED = F
IRES = DV-IPRED
PROP=SQRT(SIGMA(1,1))*IPRED
ADD=SQRT(SIGMA(2,2))
; COV = SIGMA(1,2); if covariance is required between two epsilons
; SD=SQRT(PROP*PROP + ADD*ADD + 2*COV) ; if covariance is required between two epsilons
SD=SQRT(PROP*PROP + ADD*ADD)
IWRES = IRES/SD
Y = IPRED+IPRED*EPS(1) + EPS(2)
$SIGMA 0.02 ; proportional error
1 ; additive error
$ESTIMATION METHOD=1 INTERACTION MAXEVAL=9999 SIG=3 PRINT=5 NOABORT
POSTHOC MSFO=msf105
$COVARIANCE PRINT=E
;$SIM (1234) NSUBPROBLEMS=1 ONLYSIM
$TABLE ID TIME DV MDV EVID IPRED IWRES ONEHEADER NOPRINT
FILE=sdtab105
$TABLE ID CL V ETA1 ETA2 ONEHEADER NOPRINT FILE=patab105
$TABLE ID SEX ETN ONEHEADER NOPRINT FILE=catab105
$TABLE ID WT ONEHEADER NOPRINT FILE=cotab105
$TABLE ID CL V SEX ETN WT ETA1 ETA2 ONEHEADER NOPRINT
FILE=mytab105
NM-TRAN MESSAGES
WARNINGS AND ERRORS (IF ANY) FOR PROBLEM 1
(WARNING 2) NM-TRAN INFERS THAT THE DATA ARE POPULATION.
(WARNING 79) SIGMA IS USED ON THE RIGHT. WITH A SUBSEQUENT RUN, IF AN
INITIAL ESTIMATE OF A DIAGONAL BLOCK OF SIGMA IS TO BE COMPUTED BY
NONMEM, THAT BLOCK WILL BE SET TO AN IDENTITY MATRIX DURING THAT
COMPUTATION. THIS COULD LEAD TO AN ARITHMETIC EXCEPTION.*
* THE MAXIMUM NUMBER OF WARNINGS OF ONE OR MORE TYPES WAS REACHED.
IT IS POSSIBLE THAT SOME WARNING MESSAGES WERE SUPPRESSED.
Note: The following floating-point exceptions are signalling: IEEE_INVALID_FLAG IEEE_DIVIDE_BY_ZERO
License Registered to: University of Maryland - Baltimore
Expiration Date: 14 JUN 2018
Current Date: 2 JUL 2017
Days until program expires : 347
1NONLINEAR MIXED EFFECTS MODEL PROGRAM (NONMEM) VERSION 7.3.0
ORIGINALLY DEVELOPED BY STUART BEAL, LEWIS SHEINER, AND ALISON BOECKMANN
CURRENT DEVELOPERS ARE ROBERT BAUER, ICON DEVELOPMENT SOLUTIONS,
AND ALISON BOECKMANN. IMPLEMENTATION, EFFICIENCY, AND STANDARDIZATION
PERFORMED BY NOUS INFOSYSTEMS.
PROBLEM NO.: 1
PK model 1 cmt base
0DATA CHECKOUT RUN: NO
DATA SET LOCATED ON UNIT NO.: 2
THIS UNIT TO BE REWOUND: NO
NO. OF DATA RECS IN DATA SET: 800
NO. OF DATA ITEMS IN DATA SET: 9
ID DATA ITEM IS DATA ITEM NO.: 1
DEP VARIABLE IS DATA ITEM NO.: 5
MDV DATA ITEM IS DATA ITEM NO.: 3
0INDICES PASSED TO SUBROUTINE PRED:
4 2 6 0 0 0 0 0 0 0 0
0LABELS FOR DATA ITEMS:
ID TIME MDV EVID DV AMT SEX WT ETN
0(NONBLANK) LABELS FOR PRED-DEFINED ITEMS:
CL V IPRED IWRES
0FORMAT FOR DATA:
(E3.0,E5.0,2E2.0,E21.0,E4.0,E2.0,E5.0,E2.0)
TOT. NO. OF OBS RECS: 760
TOT. NO. OF INDIVIDUALS: 40
0LENGTH OF THETA: 4
0DEFAULT THETA BOUNDARY TEST OMITTED: NO
0OMEGA HAS SIMPLE DIAGONAL FORM WITH DIMENSION: 2
0DEFAULT OMEGA BOUNDARY TEST OMITTED: NO
0SIGMA HAS SIMPLE DIAGONAL FORM WITH DIMENSION: 2
0DEFAULT SIGMA BOUNDARY TEST OMITTED: NO
0INITIAL ESTIMATE OF THETA:
LOWER BOUND INITIAL EST UPPER BOUND
0.0000E+00 0.2000E+01 0.1000E+07
0.0000E+00 0.3000E+01 0.1000E+07
0.0000E+00 0.1000E+02 0.1000E+07
0.0000E+00 0.1000E+01 0.1000E+07
0INITIAL ESTIMATE OF OMEGA:
0.5000E-01
0.0000E+00 0.2000E+00
0INITIAL ESTIMATE OF SIGMA:
0.2000E-01
0.0000E+00 0.1000E+01
0COVARIANCE STEP OMITTED: NO
EIGENVLS. PRINTED: YES
SPECIAL COMPUTATION: NO
COMPRESSED FORMAT: NO
SIGDIGITS ETAHAT (SIGLO): -1
SIGDIGITS GRADIENTS (SIGL): -1
RELATIVE TOLERANCE (TOL): -1
ABSOLUTE TOLERANCE-ADVAN 9,13 ONLY (ATOL): -1
EXCLUDE COV FOR FOCE (NOFCOV): NO
RESUME COV ANALYSIS (RESUME): NO
0TABLES STEP OMITTED: NO
NO. OF TABLES: 5
SEED NUMBER (SEED): 11456
RANMETHOD:
MC SAMPLES (ESEED): 300
WRES SQUARE ROOT TYPE: EIGENVALUE
0-- TABLE 1 --
PRINTED: NO
HEADER: YES
FILE TO BE FORWARDED: NO
FORMAT: S1PE11.4
LFORMAT:
RFORMAT:
0USER-CHOSEN ITEMS:
ID TIME DV MDV EVID IPRED IWRES
0-- TABLE 2 --
PRINTED: NO
HEADER: YES
FILE TO BE FORWARDED: NO
FORMAT: S1PE11.4
LFORMAT:
RFORMAT:
0USER-CHOSEN ITEMS:
ID CL V ETA1 ETA2
0-- TABLE 3 --
PRINTED: NO
HEADER: YES
FILE TO BE FORWARDED: NO
FORMAT: S1PE11.4
LFORMAT:
RFORMAT:
0USER-CHOSEN ITEMS:
ID SEX ETN
0-- TABLE 4 --
PRINTED: NO
HEADER: YES
FILE TO BE FORWARDED: NO
FORMAT: S1PE11.4
LFORMAT:
RFORMAT:
0USER-CHOSEN ITEMS:
ID WT
0-- TABLE 5 --
PRINTED: NO
HEADER: YES
FILE TO BE FORWARDED: NO
FORMAT: S1PE11.4
LFORMAT:
RFORMAT:
0USER-CHOSEN ITEMS:
ID CL V SEX ETN WT ETA1 ETA2
1DOUBLE PRECISION PREDPP VERSION 7.3.0
ONE COMPARTMENT MODEL WITH FIRST-ORDER ABSORPTION (ADVAN2)
0MAXIMUM NO. OF BASIC PK PARAMETERS: 3
0BASIC PK PARAMETERS (AFTER TRANSLATION):
ELIMINATION RATE (K) IS BASIC PK PARAMETER NO.: 1
ABSORPTION RATE (KA) IS BASIC PK PARAMETER NO.: 3
TRANSLATOR WILL CONVERT PARAMETERS
CLEARANCE (CL) AND VOLUME (V) TO K (TRANS2)
0COMPARTMENT ATTRIBUTES
COMPT. NO. FUNCTION INITIAL ON/OFF DOSE DEFAULT DEFAULT
STATUS ALLOWED ALLOWED FOR DOSE FOR OBS.
1 DEPOT OFF YES YES YES NO
2 CENTRAL ON NO YES NO YES
3 OUTPUT OFF YES NO NO NO
1
ADDITIONAL PK PARAMETERS - ASSIGNMENT OF ROWS IN GG
COMPT. NO. INDICES
SCALE BIOAVAIL. ZERO-ORDER ZERO-ORDER ABSORB
FRACTION RATE DURATION LAG
1 * * * * *
2 4 * * * *
3 * - - - -
- PARAMETER IS NOT ALLOWED FOR THIS MODEL
* PARAMETER IS NOT SUPPLIED BY PK SUBROUTINE;
WILL DEFAULT TO ONE IF APPLICABLE
0DATA ITEM INDICES USED BY PRED ARE:
EVENT ID DATA ITEM IS DATA ITEM NO.: 4
TIME DATA ITEM IS DATA ITEM NO.: 2
DOSE AMOUNT DATA ITEM IS DATA ITEM NO.: 6
0PK SUBROUTINE CALLED WITH EVERY EVENT RECORD.
PK SUBROUTINE NOT CALLED AT NONEVENT (ADDITIONAL OR LAGGED) DOSE TIMES.
0ERROR SUBROUTINE CALLED WITH EVERY EVENT RECORD.
1
#TBLN: 1
#METH: First Order Conditional Estimation with Interaction
ESTIMATION STEP OMITTED: NO
ANALYSIS TYPE: POPULATION
CONDITIONAL ESTIMATES USED: YES
CENTERED ETA: NO
EPS-ETA INTERACTION: YES
LAPLACIAN OBJ. FUNC.: NO
NO. OF FUNCT. EVALS. ALLOWED: 9999
NO. OF SIG. FIGURES REQUIRED: 3
INTERMEDIATE PRINTOUT: YES
ESTIMATE OUTPUT TO MSF: YES
ABORT WITH PRED EXIT CODE 1: NO
IND. OBJ. FUNC. VALUES SORTED: NO
NUMERICAL DERIVATIVE
FILE REQUEST (NUMDER): NONE
MAP (ETAHAT) ESTIMATION METHOD (OPTMAP): 0
ETA HESSIAN EVALUATION METHOD (ETADER): 0
INITIAL ETA FOR MAP ESTIMATION (MCETA): 0
SIGDIGITS FOR MAP ESTIMATION (SIGLO): 100
GRADIENT SIGDIGITS OF
FIXED EFFECTS PARAMETERS (SIGL): 100
EXCLUDE TITLE (NOTITLE): NO
EXCLUDE COLUMN LABELS (NOLABEL): NO
NOPRIOR SETTING (NOPRIOR): OFF
NOCOV SETTING (NOCOV): OFF
DERCONT SETTING (DERCONT): OFF
ABSOLUTE TOLERANCE-ADVAN 9,13 ONLY(ATOL):-100
FINAL ETA RE-EVALUATION (FNLETA): ON
EXCLUDE NON-INFLUENTIAL (NON-INFL.) ETAS
IN SHRINKAGE (ETASTYPE): NO
NON-INFL. ETA CORRECTION (NONINFETA): OFF
FORMAT FOR ADDITIONAL FILES (FORMAT): S1PE12.5
PARAMETER ORDER FOR OUTPUTS (ORDER): TSOL
ADDITIONAL CONVERGENCE TEST (CTYPE=4)?: NO
EM OR BAYESIAN METHOD USED: NONE
THE FOLLOWING LABELS ARE EQUIVALENT
PRED=PREDI
RES=RESI
WRES=WRESI
IWRS=IWRESI
IPRD=IPREDI
IRS=IRESI
MONITORING OF SEARCH:
0ITERATION NO.: 0 OBJECTIVE VALUE: 8830.58029680930 NO. OF FUNC. EVALS.: 8
CUMULATIVE NO. OF FUNC. EVALS.: 8
NPARAMETR: 2.0000E+00 3.0000E+00 1.0000E+01 1.0000E+00 5.0000E-02 2.0000E-01 2.0000E-02 1.0000E+00
PARAMETER: 1.0000E-01 1.0000E-01 1.0000E-01 1.0000E-01 1.0000E-01 1.0000E-01 1.0000E-01 1.0000E-01
GRADIENT: -9.5465E+02 -1.4437E+03 -1.5699E+03 -4.4719E+02 -3.6312E+03 -6.1770E+03 -2.6099E+03 -7.9008E+02
0ITERATION NO.: 5 OBJECTIVE VALUE: 2455.46902797787 NO. OF FUNC. EVALS.: 9
CUMULATIVE NO. OF FUNC. EVALS.: 53
NPARAMETR: 1.8963E+00 9.5171E+00 1.8815E+01 1.6104E+00 4.9924E+00 1.3899E+01 3.4904E-03 3.8894E+00
PARAMETER: 4.6747E-02 1.2545E+00 7.3205E-01 5.7648E-01 2.4018E+00 2.2206E+00 -7.7285E-01 7.7913E-01
GRADIENT: -3.6479E+02 -1.8691E+01 -1.9157E+00 -8.5963E+00 3.8821E+01 2.0986E+01 -2.3993E+01 9.6543E+01
0ITERATION NO.: 10 OBJECTIVE VALUE: 2281.24355686989 NO. OF FUNC. EVALS.: 22
CUMULATIVE NO. OF FUNC. EVALS.: 111
NPARAMETR: 2.1773E+00 7.4813E+01 7.1590E+01 4.2487E-01 7.4494E-01 9.0126E+00 7.1261E-03 2.1641E+00
PARAMETER: 1.8495E-01 3.3164E+00 2.0684E+00 -7.5598E-01 1.4506E+00 2.0040E+00 -4.1598E-01 4.8601E-01
GRADIENT: -9.6301E+01 -1.8239E+01 -1.6583E+01 -1.8782E+01 6.7273E+01 4.7591E+01 -9.3981E+01 -7.4464E+01
0ITERATION NO.: 15 OBJECTIVE VALUE: 2139.32354446459 NO. OF FUNC. EVALS.: 15
CUMULATIVE NO. OF FUNC. EVALS.: 190
NPARAMETR: 2.2042E+00 1.0763E+02 5.2848E+02 4.2671E-01 6.9276E-01 3.1587E-01 1.2377E-02 1.8987E+00
PARAMETER: 1.9723E-01 3.6801E+00 4.0674E+00 -7.5166E-01 1.4143E+00 3.2851E-01 -1.3994E-01 4.2059E-01
GRADIENT: -6.5776E+01 2.2172E+01 3.1251E+01 3.9410E+00 6.5663E+01 3.8647E+01 3.0921E+01 3.7747E+01
0ITERATION NO.: 20 OBJECTIVE VALUE: 2059.90616947527 NO. OF FUNC. EVALS.: 15
CUMULATIVE NO. OF FUNC. EVALS.: 266
NPARAMETR: 2.2493E+00 7.9818E+01 4.7782E+02 5.8430E-01 5.8584E-02 1.5647E-01 1.3291E-02 1.6651E+00
PARAMETER: 2.1749E-01 3.3811E+00 3.9666E+00 -4.3735E-01 1.7922E-01 -2.2734E-02 -1.0431E-01 3.5494E-01
GRADIENT: -1.0065E+01 8.1530E+01 1.0910E+01 4.4609E+01 1.6098E+01 6.4079E+00 2.3683E+01 9.3095E+00
0ITERATION NO.: 25 OBJECTIVE VALUE: 2053.83194992076 NO. OF FUNC. EVALS.: 15
CUMULATIVE NO. OF FUNC. EVALS.: 341
NPARAMETR: 2.2659E+00 7.4145E+01 4.6844E+02 5.8807E-01 3.6541E-02 1.4406E-01 1.2250E-02 1.7269E+00
PARAMETER: 2.2484E-01 3.3074E+00 3.9468E+00 -4.3090E-01 -5.6788E-02 -6.4054E-02 -1.4510E-01 3.7317E-01
GRADIENT: -1.5411E+00 7.1997E-02 -8.9898E-02 3.5627E-01 8.2479E-03 1.9281E-01 3.4253E-01 -4.1103E-01
0ITERATION NO.: 30 OBJECTIVE VALUE: 2053.83042655628 NO. OF FUNC. EVALS.: 12
CUMULATIVE NO. OF FUNC. EVALS.: 413
NPARAMETR: 2.2676E+00 7.4170E+01 4.6858E+02 5.8761E-01 3.6515E-02 1.4371E-01 1.2208E-02 1.7332E+00
PARAMETER: 2.2556E-01 3.3078E+00 3.9471E+00 -4.3169E-01 -5.7148E-02 -6.5261E-02 -1.4681E-01 3.7499E-01
GRADIENT: -4.1394E-03 1.2620E-02 1.0149E-02 3.1581E-03 1.7556E-03 2.0821E-03 -1.0978E-02 -7.8332E-03
#TERM:
0MINIMIZATION SUCCESSFUL
NO. OF FUNCTION EVALUATIONS USED: 413
NO. OF SIG. DIGITS IN FINAL EST.: 3.7
ETABAR IS THE ARITHMETIC MEAN OF THE ETA-ESTIMATES,
AND THE P-VALUE IS GIVEN FOR THE NULL HYPOTHESIS THAT THE TRUE MEAN IS 0.
ETABAR: 1.3928E-03 -1.6503E-03
SE: 2.6256E-02 5.9356E-02
N: 40 40
P VAL.: 9.5770E-01 9.7782E-01
ETAshrink(%): 1.1992E+01 1.0000E-10
EBVshrink(%): 1.2634E+01 7.4468E-01
EPSshrink(%): 3.4665E+00 4.6819E+00
#TERE:
Elapsed estimation time in seconds: 1.96
Elapsed covariance time in seconds: 1.40
1
************************************************************************************************************************
******************** ********************
******************** FIRST ORDER CONDITIONAL ESTIMATION WITH INTERACTION ********************
#OBJT:************** MINIMUM VALUE OF OBJECTIVE FUNCTION ********************
******************** ********************
************************************************************************************************************************
#OBJV:******************************************** 2053.830 **************************************************
1
************************************************************************************************************************
******************** ********************
******************** FIRST ORDER CONDITIONAL ESTIMATION WITH INTERACTION ********************
******************** FINAL PARAMETER ESTIMATE ********************
******************** ********************
************************************************************************************************************************
THETA - VECTOR OF FIXED EFFECTS PARAMETERS *********
TH 1 TH 2 TH 3 TH 4
2.27E+00 7.42E+01 4.69E+02 5.88E-01
OMEGA - COV MATRIX FOR RANDOM EFFECTS - ETAS ********
ETA1 ETA2
ETA1
+ 3.65E-02
ETA2
+ 0.00E+00 1.44E-01
SIGMA - COV MATRIX FOR RANDOM EFFECTS - EPSILONS ****
EPS1 EPS2
EPS1
+ 1.22E-02
EPS2
+ 0.00E+00 1.73E+00
1
OMEGA - CORR MATRIX FOR RANDOM EFFECTS - ETAS *******
ETA1 ETA2
ETA1
+ 1.91E-01
ETA2
+ 0.00E+00 3.79E-01
SIGMA - CORR MATRIX FOR RANDOM EFFECTS - EPSILONS ***
EPS1 EPS2
EPS1
+ 1.10E-01
EPS2
+ 0.00E+00 1.32E+00
1
************************************************************************************************************************
******************** ********************
******************** FIRST ORDER CONDITIONAL ESTIMATION WITH INTERACTION ********************
******************** STANDARD ERROR OF ESTIMATE ********************
******************** ********************
************************************************************************************************************************
THETA - VECTOR OF FIXED EFFECTS PARAMETERS *********
TH 1 TH 2 TH 3 TH 4
8.10E-02 2.84E+00 2.83E+01 3.87E-02
OMEGA - COV MATRIX FOR RANDOM EFFECTS - ETAS ********
ETA1 ETA2
ETA1
+ 1.11E-02
ETA2
+ ......... 2.32E-02
SIGMA - COV MATRIX FOR RANDOM EFFECTS - EPSILONS ****
EPS1 EPS2
EPS1
+ 2.12E-03
EPS2
+ ......... 2.60E-01
1
OMEGA - CORR MATRIX FOR RANDOM EFFECTS - ETAS *******
ETA1 ETA2
ETA1
+ 2.90E-02
ETA2
+ ......... 3.06E-02
SIGMA - CORR MATRIX FOR RANDOM EFFECTS - EPSILONS ***
EPS1 EPS2
EPS1
+ 9.58E-03
EPS2
+ ......... 9.89E-02
1
************************************************************************************************************************
******************** ********************
******************** FIRST ORDER CONDITIONAL ESTIMATION WITH INTERACTION ********************
******************** COVARIANCE MATRIX OF ESTIMATE ********************
******************** ********************
************************************************************************************************************************
TH 1 TH 2 TH 3 TH 4 OM11 OM12 OM22 SG11 SG12 SG22
TH 1
+ 6.56E-03
TH 2
+ -2.46E-02 8.07E+00
TH 3
+ -1.08E-01 -2.35E+01 8.00E+02
TH 4
+ 2.33E-04 -6.52E-02 1.59E-01 1.50E-03
OM11
+ -2.19E-04 -3.50E-03 4.24E-02 9.27E-05 1.23E-04
OM12
+ ......... ......... ......... ......... ......... .........
OM22
+ 1.65E-04 2.72E-03 -2.73E-03 5.50E-05 -4.58E-05 ......... 5.38E-04
SG11
+ -3.12E-05 -1.16E-04 1.89E-02 -1.11E-05 8.29E-06 ......... -5.94E-06 4.48E-06
SG12
+ ......... ......... ......... ......... ......... ......... ......... ......... .........
SG22
+ 4.31E-03 -7.84E-02 -6.81E-01 9.52E-04 -1.24E-03 ......... 4.25E-04 -4.34E-04 ......... 6.78E-02
1
************************************************************************************************************************
******************** ********************
******************** FIRST ORDER CONDITIONAL ESTIMATION WITH INTERACTION ********************
******************** CORRELATION MATRIX OF ESTIMATE ********************
******************** ********************
************************************************************************************************************************
TH 1 TH 2 TH 3 TH 4 OM11 OM12 OM22 SG11 SG12 SG22
TH 1
+ 8.10E-02
TH 2
+ -1.07E-01 2.84E+00
TH 3
+ -4.69E-02 -2.93E-01 2.83E+01
TH 4
+ 7.43E-02 -5.92E-01 1.45E-01 3.87E-02
OM11
+ -2.44E-01 -1.11E-01 1.35E-01 2.16E-01 1.11E-02
OM12
+ ......... ......... ......... ......... ......... .........
OM22
+ 8.79E-02 4.13E-02 -4.16E-03 6.12E-02 -1.78E-01 ......... 2.32E-02
SG11
+ -1.82E-01 -1.92E-02 3.16E-01 -1.36E-01 3.54E-01 ......... -1.21E-01 2.12E-03
SG12
+ ......... ......... ......... ......... ......... ......... ......... ......... .........
SG22
+ 2.04E-01 -1.06E-01 -9.25E-02 9.44E-02 -4.31E-01 ......... 7.05E-02 -7.87E-01 ......... 2.60E-01
1
************************************************************************************************************************
******************** ********************
******************** FIRST ORDER CONDITIONAL ESTIMATION WITH INTERACTION ********************
******************** INVERSE COVARIANCE MATRIX OF ESTIMATE ********************
******************** ********************
************************************************************************************************************************
TH 1 TH 2 TH 3 TH 4 OM11 OM12 OM22 SG11 SG12 SG22
TH 1
+ 1.68E+02
TH 2
+ 5.25E-01 2.18E-01
TH 3
+ 1.91E-02 3.39E-03 1.63E-03
TH 4
+ -1.69E+01 9.30E+00 -6.46E-02 1.18E+03
OM11
+ 2.64E+02 -7.32E-01 -2.24E-01 -9.42E+02 1.18E+04
OM12
+ ......... ......... ......... ......... ......... .........
OM22
+ -2.48E+01 -1.90E+00 -1.18E-01 -1.83E+02 8.54E+02 ......... 1.98E+03
SG11
+ 2.38E+02 7.52E+01 -1.30E+01 7.22E+03 -1.21E+03 ......... 3.32E+03 7.55E+05
SG12
+ ......... ......... ......... ......... ......... ......... ......... ......... .........
SG22
+ -3.11E+00 6.02E-01 -6.68E-02 2.47E+01 1.97E+02 ......... 2.52E+01 4.63E+03 ......... 4.77E+01
1
************************************************************************************************************************
******************** ********************
******************** FIRST ORDER CONDITIONAL ESTIMATION WITH INTERACTION ********************
******************** EIGENVALUES OF COR MATRIX OF ESTIMATE ********************
******************** ********************
************************************************************************************************************************
1 2 3 4 5 6 7 8
1.60E-01 3.53E-01 5.63E-01 8.59E-01 9.21E-01 1.05E+00 1.79E+00 2.30E+00
#CPUT: Total CPU Time in Seconds, 3.528
Stop Time:
Sun Jul 2 07:19:25 EDT 2017