Description Usage Arguments Details Value References See Also Examples
Calculates an approximation xi^{**} of the T_P-optimal design xi^* for discrimination between a given list of models {eta_i(x,theta_i), i = 1,…,nu}. This procedure is based on the algorithms developed by Holger Dette, Viatcheslav B. Melas and Roman Guchenko in [7]. T_P-optimal design is a probability measure, which maximizes the functional
T_P(xi) = sum_ij p_ij inf_{theta_ij \in Theta_j} int_X [ eta_i(x,theta_i^) - η_j(x,theta_ij) ]^2 xi(dx),
where xi is an arbitrary design on X (it is presumed here, that X is an interval from R), P = {p_ij} is a table of non-negative weights with zeros on the diagonal (comparison table) and theta_i^ are predefined fixed parameters.
It was also shown in [7] that calculation of Bayesian T_P-optimal design, which maximizes more complicated criterion
T_P(xi) = sum_ij p_ij int_{Lambda_i} inf_{theta_ij \in Theta_j} int_X [ eta_i(x,lambda_i^) - η_j(x,theta_ij) ]^2 xi(dx) P_i(d lambda_i),
can be reduced to calculation of ordinary T_P-optimal design, when distributions P_i are discrete. That is why in this case the current function is also suitable for calculation of Bayesian designs.
1 2 3 4 5 6 7 8 9 |
x |
a numeric vector specifying support points from X for initial design. Current algorithm operates under the assumption, that X is an interval from R |
.
w |
a numeric vector specifying weights for initial design. This vector should have the same length as vector of support points. Furthermore, the weights of the design should sum to 1. If this vector is not specified, then the weights are presumed to be equal. |
eta |
a list of models between which proposed optimization should be performed. Every function from this list should be defined in the form of eta_i(x,theta_i), where x is one dimensional variable from X and θ_i is a vector of corresponding model parameters. We will refer to length of this list as nu. |
theta.fix |
a list of fixed model parameters θ_i^' from the functional T_P. This list should have the same length as the list of models. |
theta.var |
an array with two dimensions specifying initial values for parameter vectors theta_ij. The default value here is NULL, which means that initial guess is calculated automatically. |
p |
a nu x nu square table (R-matrix) containing non-negative weights for comparison. The diagonal values of this table should all be zeros. If one want to include comparison of the i'th model with fixed parameters against j'th model with variable parameters into optimization, then he/she should place non-negative weight p_ij into the table; otherwise this weight should be zero. |
x.lb |
a left bound for support points. If it is not specified, then minimal value from input vector x is taken. |
x.rb |
a right bound for support points. If it is not specified, then maximal value from input vector x is taken. |
opt |
a list of options containing such named fields:
|
Firstly, lets define
T_P(xi) = sum_ij p_ij [ eta_i(x,theta_i^) - η_j(x,theta_ij) ]^2, θ_{i,j} = arginf_{theta_{i,j} \in Theta_j} int_{X} [ eta_i(x,θ_i) - eta_j(x, theta_{i,j}) ]^2 xi(dx).
The simplified algorithm schema is as follows:
Let xi_s denotes the design obtained on the s'th iteration of the algorithm. Then
Support of the new design xi_{s+1} consists of all local maximums of function Psi(x,xi_s) on X united with the support of current design xi_s.
Weights of the new design xi_{s+1} are calculated so that the functional T_P(xi) achieves its maximum in the class of all designs with support from previous step.
Object of class “tpopt” which contains the following fields:
the numeric vector of support points from X for resulting approximation of T_P-optimal design.
the numeric vector of weights for resulting approximation of T_P-optimal design. The values of this vector sum to 1.
the numeric vector containing efficiency lower bound values by iteration. See details section for definition.
the numeric vector containing values of functional T_P by iteration.
the list of models, which is exactly the same as one from the arguments list.
the list of fixed model parameters. It goes to the result without any changes too.
the array with two dimensions specifying calculated values for parameter vectors theta_ij according to resulting design.
same as in input.
max.iter from options list.
number of iterations done.
desired efficiency from options list.
overall execution time.
[1] Atkinson A.C., Fedorov V.V. (1975) The design of experiments for discriminating between two rival models. Biometrika, vol. 62(1), pp. 57–70.
[2] Atkinson A.C., Fedorov V.V. (1975) Optimal design: Experiments for discriminating between several models. Biometrika, vol. 62(2), pp. 289–303.
[3] Dette H., Pepelyshev A. (2008) Efficient experimental designs for sigmoidal growth models. Journal of statistical planning and inference, vol. 138, pp. 2–17.
[4] Dette H., Melas V.B., Shpilev P. (2013) Robust T-optimal discriminating designs. Annals of Statistics, vol. 41(4), pp. 1693–1715.
[5] Braess D., Dette H. (2013) Optimal discriminating designs for several competing regression models. Annals of Statistics, vol. 41(2), pp. 897–922.
[6] Braess D., Dette H. (2013) Supplement to “Optimal discriminating designs for several competing regression models”. Annals of Statistics, online supplementary material.
[7] Dette H., Melas V.B., Guchenko R. (2014) Bayesian T-optimal discriminating designs. ArXiv link.
plot.tpopt
, summary.tpopt
, print.tpopt
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 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 | ### Auxiliary libraries for examples
library(mvtnorm)
### EMAX vs MM
#List of models
eta.1 <- function(x, theta.1)
theta.1[1] + theta.1[2] * x / (x + theta.1[3])
eta.2 <- function(x, theta.2)
theta.2[1] * x / (x + theta.2[2])
eta <- list(eta.1, eta.2)
#List of fixed parameters
theta.1 <- c(1, 1, 1)
theta.2 <- c(1, 1)
theta.fix <- list(theta.1, theta.2)
#Comparison table
p <- matrix(
c(
0, 1,
0, 0
), c(2, 2), byrow = TRUE)
#Design estimation
res <- tpopt(x = c(1.2, 1.5, 1.7), eta = eta, theta.fix = theta.fix, p = p,
x.lb = 1, x.rb = 2)
plot(res)
summary(res)
### Sigmoidal second
#List of models
eta.1 <- function(x, theta.1)
theta.1[1] / (1 + theta.1[2] * exp(-theta.1[3] * x)) ^ theta.1[4]
eta.2 <- function(x, theta.2)
theta.2[1] / (1 + theta.2[2] * exp(-theta.2[3] * x))
eta <- list(eta.1, eta.2)
#List of fixed parameters
theta.1 <- c(2, 5, 1, 2)
theta.2 <- c(3, 5, 0.7)
theta.fix <- list(theta.1, theta.2)
#Comparison table
p <- matrix(
c(
0, 1,
0, 0
), c(2, 2), byrow = TRUE)
#Design estimation
res <- tpopt(x = seq(0, 10), eta = eta, theta.fix = theta.fix, p = p)
plot(res)
summary(res)
### Sigmoidal first
#List of models
eta.1 <- function(x, theta.1)
theta.1[1] - theta.1[2] * exp(-theta.1[3] * x ^ theta.1[4])
eta.2 <- function(x, theta.2)
theta.2[1] - theta.2[2] * exp(-theta.2[3] * x)
eta <- list(eta.1, eta.2)
#List of fixed parameters
theta.1 <- c(2, 1, 0.8, 1.5)
theta.2 <- c(2, 1, 1)
theta.fix <- list(theta.1, theta.2)
#Comparision table
p <- matrix(
c(
0, 1,
0, 0
), c(2, 2), byrow = TRUE)
#Design estimation
res <- tpopt(x = seq(0, 10), eta = eta, theta.fix = theta.fix, p = p)
plot(res)
summary(res)
### Sigmoidal first --- Bayes
#List of fixed parameters
sigma <- sqrt(0.3)
theta.1.sigma <- matrix(
c(
sigma^2, 0,
0, sigma^2
), c(2, 2), byrow = TRUE)
grid <- expand.grid(
theta.1[1],
theta.1[2],
seq(theta.1[3] - sigma, theta.1[3] + sigma, length.out = 5),
seq(theta.1[4] - sigma, theta.1[4] + sigma, length.out = 5)
)
eta <- c(replicate(length(grid[,1]), eta.1, simplify = FALSE), eta.2)
theta.fix <- list()
for(i in 1:length(grid[,1]))
theta.fix[[length(theta.fix) + 1]] <- as.numeric(grid[i,])
theta.fix[[length(theta.fix) + 1]] <- theta.2
density.on.grid <- dmvnorm(grid[,3:4], mean = theta.1[3:4], sigma = theta.1.sigma)
density.on.grid <- density.on.grid / sum(density.on.grid)
#Comparison table
p <- rep(0,length(eta))
for(i in 1:length(grid[,1]))
p <- rbind(p, c(rep(0,length(eta) - 1), density.on.grid[i]))
p <- rbind(p, rep(0,length(eta)))
p <- p[-1,]
res <- tpopt(x = seq(0, 10), eta = eta, theta.fix = theta.fix, p = p)
plot(res)
summary(res)
### Dose response study
#List of models
eta.1 <- function(x, theta.1)
theta.1[1] + theta.1[2] * x
eta.2 <- function(x, theta.2)
theta.2[1] + theta.2[2] * x * (theta.2[3] - x)
eta.3 <- function(x, theta.3)
theta.3[1] + theta.3[2] * x / (theta.3[3] + x)
eta.4 <- function(x, theta.4)
theta.4[1] + theta.4[2] / (1 + exp((theta.4[3] - x) / theta.4[4]))
eta <- list(eta.1, eta.2, eta.3, eta.4)
#List of fixed parameters
theta.1 <- c(60, 0.56)
theta.2 <- c(60, 7/2250, 600)
theta.3 <- c(60, 294, 25)
theta.4 <- c(49.62, 290.51, 150, 45.51)
theta.fix <- list(theta.1, theta.2, theta.3, theta.4)
#Comparison table
p <- matrix(
c(
0, 0, 0, 0,
1, 0, 0, 0,
1, 1, 0, 0,
1, 1, 1 ,0
), c(4, 4), byrow = TRUE)
#Design estimation
res <- tpopt(x = seq(0, 500, 100), eta = eta, theta.fix = theta.fix, p = p)
plot(res)
summary(res)
### Dose response study --- Bayes
#List of fixed parameters
sigma <- 37
theta.4.sigma <- matrix(
c(
sigma^2, 0, 0, 0,
0, sigma^2, 0, 0,
0, 0, sigma^2, 0,
0, 0, 0, sigma^2
), c(4, 4), byrow = TRUE)
grid <- expand.grid(
seq(theta.4[1] - sigma, theta.4[1] + sigma, length.out = 3),
seq(theta.4[2] - sigma, theta.4[2] + sigma, length.out = 3),
seq(theta.4[3] - sigma, theta.4[3] + sigma, length.out = 3),
seq(theta.4[4] - sigma, theta.4[4] + sigma, length.out = 3)
)
eta <- c(eta.1, eta.2, eta.3, replicate(length(grid[,1]), eta.4, simplify = FALSE))
theta.fix <- list(theta.1, theta.2, theta.3)
for(i in 1:length(grid[,1]))
theta.fix[[length(theta.fix) + 1]] <- as.numeric(grid[i,])
density.on.grid <- dmvnorm(grid, mean = theta.4, sigma = theta.4.sigma)
density.on.grid <- density.on.grid / sum(density.on.grid)
#Comparison table
p <- rbind(
rep(0, length(eta)),
c(1, rep(0, length(eta) - 1)),
c(1, 1, rep(0,length(eta) - 2))
)
for(i in 1:length(grid[,1]))
p <- rbind(p, c(rep(density.on.grid[i], 3), rep(0, length(eta) - 3)))
#Design estimation
## Not run:
res <- tpopt(x = seq(0, 500, 100), eta = eta, theta.fix = theta.fix, p = p)
## End(Not run)
plot(res)
summary(res)
## Not run:
### Example from [8]
### An example of how case 2 can be computed for example 1 in [8] with tpopt function
#List of models
eta.1 <- function(x, theta.1)
log(theta.1[1] * x + theta.1[2] * x / (x + theta.1[3]))
eta.2 <- function(x, theta.2)
log(theta.2[1] * x / (x + theta.2[2]))
eta <- list(eta.1, eta.2)
#List of fixed parameters
theta.1 <- c(1, 1, 1)
theta.2 <- c(1, 1)
theta.fix <- list(theta.1, theta.2)
#Comparison table
p <- matrix(
c(
0,1,
0,0
), c(length(eta), length(eta)), byrow = TRUE)
#Case 2, method 1
#Design estimation
res <- tpopt(
x = seq(0.1, 5, length.out = 10),
eta = eta, theta.fix = theta.fix, p = p, x.lb = 0.1, x.rb = 5,
opt = list(method = 1)
)
plot(res)
summary(res)
#Case 2, method 2
#Design estimation
res <- tpopt(
x = seq(0.1, 5, length.out = 10),
eta = eta, theta.fix = theta.fix, p = p, x.lb = 0.1, x.rb = 5,
opt = list(method = 2)
)
plot(res)
summary(res)
## End(Not run)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.