Description Usage Arguments Details Value Note Author(s) References See Also Examples
Gives estimates of model states and random effects η. The
function is intended to be used based on population parameters found
using PSM.estimate
or to check initial values before
parameter estimation.
1 | PSM.smooth(Model, Data, THETA, subsample = 0, trace = 0, etaList = NULL)
|
Model |
Model list.* |
Data |
Data list.* |
THETA |
Vector of population parameters used for the state estimation. |
subsample |
Number of points to estimate states in between measurements. The extra points are linearly spaced. |
trace |
Non-negative integer. If positive, tracing information on the progress of the optimization is produced. Higher values produces more tracing information. |
etaList |
Matrix where each column contains an etimate of
η_i. |
* See description in PSM.estimate.
The function produces three types of estimates.
Only past measurements are used for the state estimate at time t.
Only past and the current measurements are used for the state estimate at time t.
All measurements (both past and future) are used to form the state estimate at time t. This is usually the prefered type of state estimate.
If subsample
>0 then the data is automatically subsampled to
provide estimated of the model states between observation time points.
An unnamed list with one element for each individual. Each element contains the following elements:
Time |
Possibly subsampled time-vector corresponding to the estimated states |
Xs, Ps |
Smoothed state and state co-variance estimate |
Ys |
Response based on smoothed state: Ys = g(Xs). |
Xf, Pf |
Filtered state and state co-variance estimate |
Xp, Pp |
Predicted state and state co-variance estimate |
Yp, R |
Predicted observations and observation variances |
eta |
Estimated eta |
etaSE |
Standard errors of eta |
negLogL |
Value of the negative log-likelihood function at
|
For further details please also read the package vignette pdf-document
by writing vignette("PSM")
in R.
Stig B. Mortensen, Soeren Klim, and Robert Miller
Please visit http://www.imm.dtu.dk/psm or refer to the
help page for PSM
.
PSM
, PSM.estimate
,
PSM.simulate
, PSM.plot
, PSM.template
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 | #detailed examples are provided in the package vignette
#Theophylline data from Boeckmann et al (1994)
#objective: recover the administered doses
library(datasets)
data(Theoph)
#reshape data to PSM format
TheophPSM = list()
for(i in 1:length(unique(Theoph$Subject))){
TheophPSM[[i]] = with(
Theoph[Theoph$Subject == i,],
list(Y = matrix(conc, nrow=1), Time = Time)
)
}
#specify a simple pharmacokinetic model comprised of
#2 state equations and 1 observation equation
#initial value of 1 state eq. varies randomly across individuals
mod = vector(mode="list")
mod$Matrices = function(phi) {
list(
matA=matrix(c(-phi$ka, 0, phi$ka, -phi$ke), nrow=2, ncol=2, byrow=TRUE),
matC=matrix(c(0, 1), nrow=1, ncol=2)
)
}
mod$h = function(eta, theta, covar) {
phi = theta
phi$dose = theta$dose * exp(eta[1])
phi
}
mod$S = function(phi) {
matrix(c(phi$sigma), nrow=1, ncol=1)
}
mod$SIG = function(phi) {
matrix(c(0, 0, 0, phi$omega), nrow=2, ncol=2, byrow=TRUE)
}
mod$X0 = function(Time, phi, U) {
matrix(c(phi$dose, 0), nrow=2, ncol=1)
}
mod$ModelPar = function(THETA) {
list(theta=list(dose = THETA["dose"], ka = THETA["ka"], ke = THETA["ke"],
omega = THETA["omega"], sigma = THETA["sigma"]),
OMEGA=matrix(c(THETA["BSV_dose"]), nrow=1, ncol=1)
)
}
#specify the search space of the fitting algorithm
parM = c(ka = 1.5, ke = 0.1, dose = 10, omega = .3, sigma = 1,
BSV_dose = 0.015)
pars = list(LB=parM*.25, Init=parM, UB=parM*1.75)
#fit model and predict data
fit = PSM.estimate(mod, TheophPSM, pars, trace = 1, fast = TRUE,
control=list(optimizer="optim", maxit=1))
pred = PSM.smooth(mod, TheophPSM, fit$THETA)
#visualize recovery performance
true_dose = tapply(Theoph$conc, Theoph$Subject, mean)
true_dose = true_dose[order(as.numeric(names(true_dose)))]
est_dose = fit$THETA["dose"] * exp(unlist(lapply(pred, function(x) x$eta)))
plot(true_dose, est_dose,
xlab="actually administered dose", ylab= "recovered dose")
abline(lm(est_dose ~ true_dose), lty=2)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.