Re: Condition number
Dear Bill,
You can see how NONMEM calculates in R with nmw package.
install.packages("nmw")
DataAll = Theoph
colnames(DataAll) = c("ID", "BWT", "DOSE", "TIME", "DV")
DataAll[,"ID"] = as.numeric(as.character(DataAll[,"ID"]))
nTheta = 3
nEta = 3
nEps = 2
THETAinit = c(2, 50, 0.1)
OMinit = matrix(c(0.2, 0.1, 0.1, 0.1, 0.2, 0.1, 0.1, 0.1, 0.2), nrow=nEta,
ncol=nEta)
SGinit = diag(c(0.1, 0.1))
LB = rep(0, nTheta) # Lower bound
UB = rep(1000000, nTheta) # Upper bound
FGD = deriv(~DOSE/(TH2*exp(ETA2))*TH1*exp(ETA1)/(TH1*exp(ETA1) -
TH3*exp(ETA3))*
(exp(-TH3*exp(ETA3)*TIME)-exp(-TH1*exp(ETA1)*TIME)),
c("ETA1","ETA2","ETA3"),
function.arg=c("TH1", "TH2", "TH3", "ETA1", "ETA2", "ETA3",
"DOSE", "TIME"),
func=TRUE, hessian=TRUE)
H = deriv(~F + F*EPS1 + EPS2, c("EPS1", "EPS2"), function.arg=c("F",
"EPS1", "EPS2"), func=TRUE)
PRED = function(THETA, ETA, DATAi)
{
FGDres = FGD(THETA[1], THETA[2], THETA[3], ETA[1], ETA[2], ETA[3],
DOSE=320, DATAi[,"TIME"])
Gres = attr(FGDres, "gradient")
Hres = attr(H(FGDres, 0, 0), "gradient")
if (e$METHOD == "LAPL") {
Dres = attr(FGDres, "hessian")
Res = cbind(FGDres, Gres, Hres, Dres[,1,1], Dres[,2,1], Dres[,2,2],
Dres[,3,])
colnames(Res) = c("F", "G1", "G2", "G3", "H1", "H2", "D11", "D21",
"D22", "D31", "D32", "D33")
} else {
Res = cbind(FGDres, Gres, Hres)
colnames(Res) = c("F", "G1", "G2", "G3", "H1", "H2")
}
return(Res)
}
####### First Order Approximation Method
InitStep(DataAll, THETAinit=THETAinit, OMinit=OMinit, SGinit=SGinit, LB=LB,
UB=UB,
Pred=PRED, METHOD="ZERO")
(EstRes = EstStep()) # 4 sec
(CovRes = CovStep()) # 2 sec
PostHocEta() # Using e$FinalPara from EstStep()
TabStep() # Commented out for the CRAN CPU time
######## First Order Conditional Estimation with Interaction Method
InitStep(DataAll, THETAinit=THETAinit, OMinit=OMinit, SGinit=SGinit, LB=LB,
UB=UB,
Pred=PRED, METHOD="COND")
(EstRes = EstStep()) # 2 min
(CovRes = CovStep()) # 1 min
get("EBE", envir=e)
TabStep()
######## Laplacian Approximation with Interacton Method
InitStep(DataAll, THETAinit=THETAinit, OMinit=OMinit, SGinit=SGinit, LB=LB,
UB=UB,
Pred=PRED, METHOD="LAPL")
(EstRes = EstStep()) # 4 min
(CovRes = CovStep()) # 1 min
get("EBE", envir=e)
TabStep()
Best regards,
Kyun-Seop Bae
Quoted reply history
On Thu, 1 Dec 2022 at 04:38, Bill Denney <[email protected]>
wrote:
> Hi everyone,
>
>
>
> This has been a great discussion!
>
>
>
> Bob: I’d like to clarify something that Pete, Matt, Ken, and Leonid were
> discussing about how the covariance matrix is calculated. I believe that
> NONMEM rescales the values for estimation and then reverses the rescaling
> for reporting. Is the covariance matrix calculated on the rescaled values
> or on the final parameter estimate values?
>
>
>
> Thanks,
>
>
>
> Bill
>
>
>
> *From:* [email protected] <[email protected]> *On
> Behalf Of *Bauer, Robert
> *Sent:* Wednesday, November 30, 2022 1:53 PM
> *To:* '[email protected]' <[email protected]>
> *Subject:* RE: [NMusers] Condition number
>
>
>
> Hello all:
>
> Report of non-positive definiteness or negative eigenvalues, are reported
> during the analysis of the R matrix (decomposition and inversion), which
> occurs before the correlation matrix is constructed. Often, this is
> caused by numerical imprecision. If the R matrix step fails, the $COV step
> fails to produce a final variance-covariance matrix, and of course, does
> not produce a correlation matrix. If the R matrix inversion step succeeds,
> the variance-covariance matrix and its correlation matrix are produced, and
> the correlation matrix is then assessed for its eigenvalues. So, both the
> R matrix (first step) and correlation matrix (second step) are decomposed
> and assessed.
>
>
>
> Robert J. Bauer, Ph.D.
>
> Senior Director
>
> Pharmacometrics R&D
>
> ICON Early Phase
>
> 731 Arbor way, suite 100
>
> Blue Bell, PA 19422
>
> Office: (215) 616-6428
>
> Mobile: (925) 286-0769
>
> [email protected]
>
> www.iconplc.com
>
>
>
>
>
>