# Theoph dataset is available to you directly from R # You will need to load the three following packages library(lattice) library(nlme) library(nlmeODE) # Plot the data: one curve per subject xyplot(conc ~ Time, groups = Subject, data = Theoph, type = "l", ylab = "Concentration (mg / L)", xlab = "Time (hr)") ### Fit nonlinear mixed model with analytical solution ### fitME = nlme(conc ~ SSfol(Dose, Time, lKe, lKa, lCl), # R has a built-in function for this compartmental model, see ?SSfol random = pdDiag(lKa + lCl ~ 1), # diagonal cov structure for the two random effects data = Theoph) # data # Parameter estimates summary(fitME) # Model fit plot(augPred(fitME, level = 0:1)) ### Fit nonlinear mixed model with numerical solution ### # Create data with variables for nlme ode model TheophODE = Theoph # Dose is at time zero. Set dose at other times to zero TheophODE$Dose[TheophODE$Time!= 0 ] = 0 # Cmt is the dosing compartment. For this example, compartment 1 TheophODE$Cmt = rep(1, dim(TheophODE)[1]) # Create Model Object OneComp = list( DiffEq=list( # DiffEq requires a list of formulas containing the ode's dy1dt = ~ -ka*y1, dy2dt = ~ ka*y1-ke*y2), ObsEq=list( # ObsEq requires a list of formulas with states that are observed and scaling parameters c1 = ~ 0, c2 = ~ y2/CL*ke), Parms=c("ka","ke","CL"), # A vector with the names of the parameters used in DiffEq, ObsEq, and Init States=c("y1","y2"), # A vector with the names of the states in DiffEq Init=list(0, 0) # A list with the same length as States specifying the initial states of the system # It is a logical vector specifying whether initial state estimates should be obtained # See ?nlmeODE for more info # No theophylline in the two compartments before drug administration ) # Create the nonlinear ode mixed model TheophModel = nlmeODE(OneComp, TheophODE, LogParms = TRUE) # By default, parameters are on a log scale. If want to change, use LogParms = FALSE # Fit model using nlme fitMEode = nlme(conc ~ TheophModel(ka, ke, CL, Time, Subject), # use the created ode model object data = TheophODE, # data fixed = ka + ke + CL ~ 1, # fixed effects random = pdDiag(ka + CL ~ 1), # random effects start = c(ka = -2.52, ke =0.40, CL = -3.25)) # starting values # Parameter Estimates summary(fitMEode) # Model fit plot(augPred(fitMEode, level=0:1))