# NOT RUN {
#################################
# Exponential growth model
################################
# See how we set stopping rule by adjusting 'stop_rule', 'checkfreq' and 'stoptol'
# It calls the 'senslocally' function every checkfreq = 50 iterations to
# calculate the ELB. if ELB is greater than stoptol = .95, then the algoithm stops.
# initializing by one iteration
res1 <- locally(formula = ~a + exp(-b*x), predvars = "x", parvars = c("a", "b"),
lx = 0, ux = 1, inipars = c(1, 10),
iter = 1, k = 2,
ICA.control= ICA.control(rseed = 100,
stop_rule = "equivalence",
checkfreq = 20, stoptol = .95))
# }
# NOT RUN {
# update the algorithm
res1 <- update(res1, 150)
#stops at iteration 21 because ELB is greater than .95
# }
# NOT RUN {
### fixed x, lx and ux are only required for equivalence theorem
# }
# NOT RUN {
res1.1 <- locally(formula = ~a + exp(-b*x), predvars = "x", parvars = c("a", "b"),
lx = 0, ux = 1, inipars = c(1, 10),
iter = 100,
x = c(.25, .5, .75),
ICA.control= ICA.control(rseed = 100))
plot(res1.1)
# we can not have an optimal design using this x
# }
# NOT RUN {
################################
## two parameter logistic model
################################
res2 <- locally(formula = ~1/(1 + exp(-b *(x - a))),
predvars = "x", parvars = c("a", "b"),
family = binomial(), lx = -3, ux = 3,
inipars = c(1, 3), iter = 1, k = 2,
ICA.control= list(rseed = 100, stop_rule = "equivalence",
checkfreq = 50, stoptol = .95))
# }
# NOT RUN {
res2 <- update(res2, 100)
# stops at iteration 51
# }
# NOT RUN {
################################
# A model with two predictors
################################
# mixed inhibition model
# }
# NOT RUN {
res3 <- locally(formula = ~ V*S/(Km * (1 + I/Kic)+ S * (1 + I/Kiu)),
predvars = c("S", "I"),
parvars = c("V", "Km", "Kic", "Kiu"),
family = gaussian(),
lx = c(0, 0), ux = c(30, 60),
k = 4,
iter = 300,
inipars = c(1.5, 5.2, 3.4, 5.6),
ICA.control= list(rseed = 100, stop_rule = "equivalence",
checkfreq = 50, stoptol = .95))
# stops at iteration 100
# }
# NOT RUN {
# }
# NOT RUN {
# fixed x
res3.1 <- locally(formula = ~ V*S/(Km * (1 + I/Kic)+ S * (1 + I/Kiu)),
predvars = c("S", "I"),
parvars = c("V", "Km", "Kic", "Kiu"),
family = gaussian(),
lx = c(0, 0), ux = c(30, 60),
iter = 100,
x = c(20, 4, 20, 4, 10, 0, 0, 30, 3, 2),
inipars = c(1.5, 5.2, 3.4, 5.6),
ICA.control= list(rseed = 100))
# }
# NOT RUN {
###################################
# user-defined optimality criterion
##################################
# When the model is defined by the formula interface
# A-optimal design for the 2PL model.
# the criterion function must have argument x, w fimfunc and the parameters defined in 'parvars'.
# use 'fimfunc' as a function of the design points x, design weights w and
# the 'parvars' parameters whenever needed.
Aopt <-function(x, w, a, b, fimfunc){
sum(diag(solve(fimfunc(x = x, w = w, a = a, b = b))))
}
## the sensitivtiy function
# xi_x is a design that put all its mass on x in the definition of the sensitivity function
# x is a vector of design points
Aopt_sens <- function(xi_x, x, w, a, b, fimfunc){
fim <- fimfunc(x = x, w = w, a = a, b = b)
M_inv <- solve(fim)
M_x <- fimfunc(x = xi_x, w = 1, a = a, b = b)
sum(diag(M_inv %*% M_x %*% M_inv)) - sum(diag(M_inv))
}
res4 <- locally(formula = ~1/(1 + exp(-b * (x-a))), predvars = "x",
parvars = c("a", "b"), family = "binomial",
lx = -3, ux = 3, inipars = c(1, 1.25),
iter = 1, k = 2,
crtfunc = Aopt,
sensfunc = Aopt_sens,
ICA.control = list(checkfreq = Inf))
# }
# NOT RUN {
res4 <- update(res4, 50)
# }
# NOT RUN {
# When the FIM of the model is defined directly via the argument 'fimfunc'
# the criterion function must have argument x, w fimfunc and param.
# use 'fimfunc' as a function of the design points x, design weights w
# and param whenever needed.
Aopt2 <-function(x, w, param, fimfunc){
sum(diag(solve(fimfunc(x = x, w = w, param = param))))
}
## the sensitivtiy function
# xi_x is a design that put all its mass on x in the definition of the sensitivity function
# x is a vector of design points
Aopt_sens2 <- function(xi_x, x, w, param, fimfunc){
fim <- fimfunc(x = x, w = w, param = param)
M_inv <- solve(fim)
M_x <- fimfunc(x = xi_x, w = 1, param = param)
sum(diag(M_inv %*% M_x %*% M_inv)) - sum(diag(M_inv))
}
res4.1 <- locally(fimfunc = FIM_logistic,
lx = -3, ux = 3, inipars = c(1, 1.25),
iter = 1, k = 2,
crtfunc = Aopt2,
sensfunc = Aopt_sens2,
ICA.control = list(checkfreq = Inf))
# }
# NOT RUN {
res4.1 <- update(res4.1, 50)
plot(res4.1)
# }
# NOT RUN {
# locally c-optimal design
# example from Chaloner and Larntz (1989) Figure 3
c_opt <-function(x, w, a, b, fimfunc){
gam <- log(.95/(1-.95))
M <- fimfunc(x = x, w = w, a = a, b = b)
c <- matrix(c(1, -gam * b^(-2)), nrow = 1)
B <- t(c) %*% c
sum(diag(B %*% solve(M)))
}
c_sens <- function(xi_x, x, w, a, b, fimfunc){
gam <- log(.95/(1-.95))
M <- fimfunc(x = x, w = w, a = a, b = b)
M_inv <- solve(M)
M_x <- fimfunc(x = xi_x, w = 1, a = a, b = b)
c <- matrix(c(1, -gam * b^(-2)), nrow = 1)
B <- t(c) %*% c
sum(diag(B %*% M_inv %*% M_x %*% M_inv)) - sum(diag(B %*% M_inv))
}
res4.2 <- locally(formula = ~1/(1 + exp(-b * (x-a))), predvars = "x",
parvars = c("a", "b"), family = "binomial",
lx = -1, ux = 1, inipars = c(0, 7),
iter = 1, k = 2,
crtfunc = c_opt, sensfunc = c_sens,
ICA.control = list(rseed = 1, checkfreq = Inf))
# }
# NOT RUN {
res4.2 <- update(res4.2, 100)
# }
# NOT RUN {
# }
Run the code above in your browser using DataLab