sensminimax: Verifying Optimality of The Minimax and Standardized maximin...
In ICAOD: Optimal Designs for Nonlinear Statistical Models by Imperialist Competitive Algorithm (ICA)
Description
Usage
Arguments
Details
Value
Note
References
Examples
Description
It plots the sensitivity (derivative) function of the minimax or
standardized maximin D-optimal criterion
at a given approximate (continuous) design and also
calculates its efficiency lower bound (ELB) with respect
to the optimality criterion.
For an approximate (continuous) design, when the design space is one or two-dimensional,
the user can visually verify the optimality of the design by observing the
sensitivity plot. Furthermore, the proximity of the design to the optimal design
can be measured by the ELB without knowing the latter.
See, for more details, Masoudi et al. (2017).
Usage
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
sensminimax(
formula,
predvars,
parvars,
family = gaussian(),
x,
w,
lx,
ux,
lp,
up,
fimfunc = NULL,
standardized = FALSE,
localdes = NULL,
sens.control = list(),
sens.minimax.control = list(),
calculate_criterion = TRUE,
crt.minimax.control = list(),
plot_3d = c("lattice", "rgl"),
plot_sens = TRUE,
npar = length(lp),
silent = FALSE,
crtfunc = NULL,
sensfunc = NULL
)
Arguments
formula
A linear or nonlinear model formula.
A symbolic description of the model consists of predictors and the unknown model parameters.
Will be coerced to a formula if necessary.
predvars
A vector of characters. Denotes the predictors in the formula.
parvars
A vector of characters. Denotes the unknown parameters in the formula.
family
A description of the response distribution and the link function to be used in the model.
This can be a family function, a call to a family function or a character string naming the family.
Every family function has a link argument allowing to specify the link function to be applied on the response variable.
If not specified, default links are used. For details see family.
By default, a linear gaussian model gaussian() is applied.
x
Vector of the design (support) points. See 'Details' of sensminimax for models with more than one predictors.
w
Vector of the corresponding design weights for x.
lx
Vector of lower bounds for the predictors. Should be in the same order as predvars.
ux
Vector of upper bounds for the predictors. Should be in the same order as predvars.
lp
Vector of lower bounds for the model parameters. Should be in the same order as parvars or param in the argument fimfunc.
up
Vector of upper bounds for the model parameters. Should be in the same order as parvars or param in the argument fimfunc.
When a parameter is known (has a fixed value), its associated lower and upper bound values in lp and up must be set equal.
fimfunc
A function. Returns the FIM as a matrix. Required when formula is missing. See 'Details' of minimax.
standardized
Maximin standardized design? When standardized = TRUE, the argument localdes must be given.
Defaults to FALSE. See 'Details' of minimax.
localdes
A function that takes the parameter values as inputs and returns the design points and weights of the locally optimal design.
Required when standardized = "TRUE". See 'Details' of minimax.
sens.control
Control Parameters for Calculating the ELB. For details, see sens.control.
sens.minimax.control
Control parameters to construct the answering set required for verify the general equivalence theorem and calculating the ELB. For details, see the function sens.minimax.control.
calculate_criterion
Calculate the optimality criterion? See 'Details' of sensminimax.
crt.minimax.control
Control parameters to calculate the value of the minimax or standardized maximin optimality criterion over the continuous parameter space.
Only applicable when calculate_criterion = TRUE.
For more details, see crt.minimax.control.
plot_3d
Which package should be used to plot the sensitivity (derivative) function for models with two predictors.
Either "rgl" or "lattice" (default).
plot_sens
Plot the sensitivity (derivative) function? Defaults to TRUE.
npar
Number of model parameters. Used when fimfunc is given instead of formula to specify the number of model parameters.
If not specified truly, the sensitivity (derivative) plot may be shifted below the y-axis. When NULL (default), it is set to length(lp).
silent
Do not print anything? Defaults to FALSE.
crtfunc
(Optional) a function that specifies an arbitrary criterion. It must have especial arguments and output. See 'Details' of minimax.
sensfunc
(Optional) a function that specifies the sensitivity function for crtfunc. See 'Details' of minimax.
Details
Let the unknown parameters belong to Θ.
A design ξ* is minimax D-optimal among all designs on χ if and only if there exists a probability measure μ* on
A(ξ*) = {ν belongs to Θ where -log|M(ξ*, ν) = maxima of function -log|M(ξ*, θ)| with respect to θ over Θ} ,
such that the following inequality holds for all x belong to χ
c(x, μ*, ξ*) = integration over A(ξ*) with integrand tr M^-1(ξ*, ν)I(x, ν)μ* d(ν)-p <= 0,
with equality at all support points of ξ*.
Here, p is the number of model parameters. c(x, μ*, ξ*) is called sensitivity or derivative function.
The set A(ξ*) is sometimes called answering set of
ξ* and the measure μ* is a sub-gradient of the
non-differentiable criterion evaluated at M(ξ*,ν).
For the standardized maximin D-optimal designs, the answering set N(ξ*) is
N(ξ*) = {ν belongs to Θ where eff_D(ξ*, ν) = minima of function eff_D(ξ*, θ) with respect to θ over Θ},
where
eff_D(ξ, θ) = (|M(ξ, θ)|/|M(ξ_θ, θ)|)^(1/p) and ξ_θ is the locally D-optimal design with respect to θ.
See 'Details' of sens.minimax.control on how we find the answering set.
The argument x is the vector of design points.
For design points with more than one dimension (the models with more than one predictors),
it is a concatenation of the design points, but dimension-wise.
For example, let the model has three predictors (I, S, Z).
Then, a two-point optimal design has the following points:
{point1 = (I1, S1, Z1), point2 = (I2, S2, Z2)}.
Then, the argument x is equal to
x = c(I1, I2, S1, S2, Z1, Z2).
ELB is a measure of proximity of a design to the optimal design without knowing the latter.
Given a design, let ε be the global maximum
of the sensitivity (derivative) function with respect x where x belong to χ.
ELB is given by
ELB = p/(p + ε),
where p is the number of model parameters. Obviously,
calculating ELB requires finding ε and
another optimization problem to be solved.
The tuning parameters of this optimization can be regulated via the argument sens.minimax.control.
See, for more details, Masoudi et al. (2017).
The criterion value for the minimax D-optimal design is (global maximum over Θ)
max -log|M(ξ, θ)|;
for the standardized maximin D-optimal design is (global minimum over Θ)
inf {|M(ξ, θ)| / |M(ξ_θ, θ)|}^p.
This function confirms the optimality assuming only a continuous parameter space Θ.
Value
an object of class sensminimax that is a list with the following elements:
typeArgument type that is required for print methods.
optimaA matrix that stores all the local optima over the parameter space.
The cost (criterion) values are stored in a column named Criterion_Value.
The last column (Answering_Set)
shows if the optimum belongs to the answering set (1) or not (0). See 'Details' of sens.minimax.control.
Only applicable for minimax or standardized maximin designs.
muProbability measure on the answering set.
Corresponds to the rows of optima for which the associated row in column Answering_Set is equal to 1.
Only applicable for minimax or standardized maximin designs.
max_derivGlobal maximum of the sensitivity (derivative) function (ε in 'Details').
ELBD-efficiency lower bound. Can not be larger than 1. If negative, see 'Note' in sensminimax or sens.minimax.control.
merge_tolMerging tolerance to create the answering set from the set of all local optima. See 'Details' in sens.minimax.control.
Only applicable for minimax or standardized maximin designs.
crtvalCriterion value. Compare it with the column Crtiterion_Value in optima for minimax and standardized maximin designs.
timeUsed CPU time (rough approximation).
Note
Theoretically, ELB can not be larger than 1. But if so, it may have one of the following reasons:
-
max_deriv is not a GLOBAL maximum. Please increase the value of the parameter maxeval in sens.minimax.control to find the global maximum.
The sensitivity function is shifted below the y-axis because
the number of model parameters has not been specified correctly (less value given).
Please specify the correct number of model parameters via argument npar.
Please increase the value of the parameter
n_seg in sens.minimax.control
for models with larger number of parameters or large parameter space to find the true
answering set for minimax and standardized maximin designs. See sens.minimax.control for more details.
References
Masoudi E, Holling H, Wong W.K. (2017). Application of Imperialist Competitive Algorithm to Find Minimax and Standardized Maximin Optimal Designs. Computational Statistics and Data Analysis, 113, 330-345.
Examples
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
##########################
# Power logistic model
##########################
# verifying the minimax D-optimality of a design with points x0 and weights w0
x0 <- c(-4.5515, 0.2130, 2.8075)
w0 <- c(0.4100, 0.3723, 0.2177)
# Power logistic model when s = .2
sensminimax(formula = ~ (1/(1 + exp(-b * (x-a))))^.2,
predvars = "x",
parvars = c("a", "b"),
family = binomial(),
x = x0, w = w0,
lx = -5, ux = 5,
lp = c(0, 1), up = c(3, 1.5))
##############################
# A model with two predictors
##############################
# Verifying the minimax D-optimality of a design for a model with two predictors
# The model is the mixed inhibition model.
# X0 is the vector of four design points that are:
# (3.4614, 0) (4.2801, 3.1426) (30, 0) (30, 4.0373)
x0 <- c(3.4614, 4.2801, 30, 30, 0, 3.1426, 0, 4.0373)
w0 <- rep(1/4, 4)
sensminimax(formula = ~ V*S/(Km * (1 + I/Kic)+ S * (1 + I/Kiu)),
predvars = c("S", "I"),
parvars = c("V", "Km", "Kic", "Kiu"),
family = "gaussian",
x = x0, w = w0,
lx = c(0, 0), ux = c(30, 60),
lp = c(1, 4, 2, 4), up = c(1, 5, 3, 5))
##########################################
# Standardized maximin D-optimal designs
##########################################
# Verifying the standardized maximin D-optimality of a design for
# the loglinear model
# First we should define the function for 'localdes' argument
# The function LDOD takes the parameters and returns the points and
# weights of the locally D-optimal design
LDOD <- function(theta0, theta1, theta2){
## param is the vector of theta = (theta0, theta1, theta2)
lx <- 0 # lower bound of the design space
ux <- 150 # upper bound of the design space
param <- c()
param[1] <- theta0
param[2] <- theta1
param[3] <- theta2
xstar <- (ux+param[3]) * (lx + param[3]) *
(log(ux + param[3]) - log(lx + param[3]))/(ux - lx) - param[3]
return(list(x = c(lx, xstar, ux) , w = rep(1/3, 3)))
}
x0 <- c(0, 4.2494, 17.0324, 149.9090)
w0 <- c(0.3204, 0.1207, 0.2293, 0.3296)
## Not run:
sensminimax(formula = ~theta0 + theta1* log(x + theta2),
predvars = c("x"),
parvars = c("theta0", "theta1", "theta2"),
x = x0, w = w0,
lx = 0, ux = 150,
lp = c(2, 2, 1), up = c(2, 2, 15),
localdes = LDOD,
standardized = TRUE,
sens.minimax.control = list(n_seg = 10))
## End(Not run)
################################################################
# Not necessary!
# The rest of the examples here are only for professional uses.
################################################################
# Imagine you have written your own FIM, say in Rcpp that is faster than
# the FIM created by the formula interface here.
##########################
# Power logistic model
##########################
# For example, th cpp FIM function for the power logistic model is named:
FIM_power_logistic
args(FIM_power_logistic)
# The arguments do not match the standard of the argument 'fimfunc'
# in 'sensminimax'
# So we reparameterize it:
myfim1 <- function(x, w, param)
FIM_power_logistic(x = x, w = w, param =param, s = .2)
args(myfim1)
## Not run:
# Verify minimax D-optimality of a design
sensminimax(fimfunc = myfim1,
x = c(-4.5515, 0.2130, 2.8075),
w = c(0.4100, 0.3723, 0.2177),
lx = -5, ux = 5,
lp = c(0, 1), up = c(3, 1.5))
## End(Not run)
##############################
# A model with two predictors
##############################
# An example of a model with two-predictors: mixed inhibition model
# Fisher information matrix:
FIM_mixed_inhibition
args(FIM_mixed_inhibition)
# We should first reparameterize the FIM to match the standard of the
# argument 'fimfunc'
myfim2 <- function(x, w, param){
npoint <- length(x)/2
S <- x[1:npoint]
I <- x[(npoint+1):(npoint*2)]
out <- FIM_mixed_inhibition(S = S, I = I, w = w, param = param)
return(out)
}
args(myfim2)
## Not run:
# Verifyng minimax D-optimality of a design
sensminimax(fimfunc = myfim2,
x = c(3.4614, 4.2801, 30, 30, 0, 3.1426, 0, 4.0373),
w = rep(1/4, 4),
lx = c(0, 0), ux = c(30, 60),
lp = c(1, 4, 2, 4), up = c(1, 5, 3, 5))
## End(Not run)
#########################################
# Standardized maximin D-optimal designs
#########################################
# An example of a user-written FIM function:
help(FIM_loglin)
# An example of verfying standardaized maximin D-optimality for a design
# Look how we re-define the function LDOD above
LDOD2 <- function(param){
## param is the vector of theta = (theta0, theta1, theta2)
lx <- 0 # lower bound of the design space
ux <- 150 # upper bound of the design space
xstar <- (ux + param[3]) * (lx + param[3]) *
(log(ux + param[3]) - log(lx + param[3]))/(ux - lx) - param[3]
return(list(x = c(lx, xstar, ux) , w = rep(1/3, 3)))
}
args(LDOD2)
sensminimax(fimfunc = FIM_loglin,
x = x0,
w = w0,
lx = 0, ux = 150,
lp = c(2, 2, 1), up = c(2, 2, 15),
localdes = LDOD2,
standardized = TRUE)
###################################
# user-defined optimality criterion
##################################
# When the model is defined by the formula interface
# Checking the A-optimality 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))
}
sensminimax(formula = ~1/(1 + exp(-b * (x-a))), predvars = "x",
parvars = c("a", "b"), family = "binomial",
lp = c(-2, 1), up = c(2, 1.5),
crtfunc = Aopt,
lx = -2, ux = 2,
sensfunc = Aopt_sens,
x = c(-2, .0033, 2), w = c(.274, .452, .274))
ICAOD documentation built on Oct. 23, 2020, 6:40 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Usage Arguments Details Value Note References Examples
Description
It plots the sensitivity (derivative) function of the minimax or standardized maximin D-optimal criterion at a given approximate (continuous) design and also calculates its efficiency lower bound (ELB) with respect to the optimality criterion. For an approximate (continuous) design, when the design space is one or two-dimensional, the user can visually verify the optimality of the design by observing the sensitivity plot. Furthermore, the proximity of the design to the optimal design can be measured by the ELB without knowing the latter. See, for more details, Masoudi et al. (2017).
Usage
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 | sensminimax(
formula,
predvars,
parvars,
family = gaussian(),
x,
w,
lx,
ux,
lp,
up,
fimfunc = NULL,
standardized = FALSE,
localdes = NULL,
sens.control = list(),
sens.minimax.control = list(),
calculate_criterion = TRUE,
crt.minimax.control = list(),
plot_3d = c("lattice", "rgl"),
plot_sens = TRUE,
npar = length(lp),
silent = FALSE,
crtfunc = NULL,
sensfunc = NULL
)
|
Arguments
formula |
A linear or nonlinear model |
predvars |
A vector of characters. Denotes the predictors in the |
parvars |
A vector of characters. Denotes the unknown parameters in the |
family |
A description of the response distribution and the link function to be used in the model.
This can be a family function, a call to a family function or a character string naming the family.
Every family function has a link argument allowing to specify the link function to be applied on the response variable.
If not specified, default links are used. For details see |
x |
Vector of the design (support) points. See 'Details' of |
w |
Vector of the corresponding design weights for |
lx |
Vector of lower bounds for the predictors. Should be in the same order as |
ux |
Vector of upper bounds for the predictors. Should be in the same order as |
lp |
Vector of lower bounds for the model parameters. Should be in the same order as |
up |
Vector of upper bounds for the model parameters. Should be in the same order as |
fimfunc |
A function. Returns the FIM as a |
standardized |
Maximin standardized design? When |
localdes |
A function that takes the parameter values as inputs and returns the design points and weights of the locally optimal design.
Required when |
sens.control |
Control Parameters for Calculating the ELB. For details, see |
sens.minimax.control |
Control parameters to construct the answering set required for verify the general equivalence theorem and calculating the ELB. For details, see the function |
calculate_criterion |
Calculate the optimality criterion? See 'Details' of |
crt.minimax.control |
Control parameters to calculate the value of the minimax or standardized maximin optimality criterion over the continuous parameter space.
Only applicable when |
plot_3d |
Which package should be used to plot the sensitivity (derivative) function for models with two predictors.
Either |
plot_sens |
Plot the sensitivity (derivative) function? Defaults to |
npar |
Number of model parameters. Used when |
silent |
Do not print anything? Defaults to |
crtfunc |
(Optional) a function that specifies an arbitrary criterion. It must have especial arguments and output. See 'Details' of |
sensfunc |
(Optional) a function that specifies the sensitivity function for |
Details
Let the unknown parameters belong to Θ. A design ξ* is minimax D-optimal among all designs on χ if and only if there exists a probability measure μ* on
A(ξ*) = {ν belongs to Θ where -log|M(ξ*, ν) = maxima of function -log|M(ξ*, θ)| with respect to θ over Θ} ,
such that the following inequality holds for all x belong to χ
c(x, μ*, ξ*) = integration over A(ξ*) with integrand tr M^-1(ξ*, ν)I(x, ν)μ* d(ν)-p <= 0,
with equality at all support points of ξ*.
Here, p is the number of model parameters. c(x, μ*, ξ*) is called sensitivity or derivative function.
The set A(ξ*) is sometimes called answering set of
ξ* and the measure μ* is a sub-gradient of the
non-differentiable criterion evaluated at M(ξ*,ν).
For the standardized maximin D-optimal designs, the answering set N(ξ*) is
N(ξ*) = {ν belongs to Θ where eff_D(ξ*, ν) = minima of function eff_D(ξ*, θ) with respect to θ over Θ},
where
eff_D(ξ, θ) = (|M(ξ, θ)|/|M(ξ_θ, θ)|)^(1/p) and ξ_θ is the locally D-optimal design with respect to θ.
See 'Details' of sens.minimax.control on how we find the answering set.
The argument x is the vector of design points.
For design points with more than one dimension (the models with more than one predictors),
it is a concatenation of the design points, but dimension-wise.
For example, let the model has three predictors (I, S, Z).
Then, a two-point optimal design has the following points:
{point1 = (I1, S1, Z1), point2 = (I2, S2, Z2)}.
Then, the argument x is equal to
x = c(I1, I2, S1, S2, Z1, Z2).
ELB is a measure of proximity of a design to the optimal design without knowing the latter. Given a design, let ε be the global maximum of the sensitivity (derivative) function with respect x where x belong to χ. ELB is given by
ELB = p/(p + ε),
where p is the number of model parameters. Obviously,
calculating ELB requires finding ε and
another optimization problem to be solved.
The tuning parameters of this optimization can be regulated via the argument sens.minimax.control.
See, for more details, Masoudi et al. (2017).
The criterion value for the minimax D-optimal design is (global maximum over Θ)
max -log|M(ξ, θ)|;
for the standardized maximin D-optimal design is (global minimum over Θ)
inf {|M(ξ, θ)| / |M(ξ_θ, θ)|}^p.
This function confirms the optimality assuming only a continuous parameter space Θ.
Value
an object of class sensminimax that is a list with the following elements:
typeArgument
typethat is required for print methods.optimaA
matrixthat stores all the local optima over the parameter space. The cost (criterion) values are stored in a column namedCriterion_Value. The last column (Answering_Set) shows if the optimum belongs to the answering set (1) or not (0). See 'Details' ofsens.minimax.control. Only applicable for minimax or standardized maximin designs.muProbability measure on the answering set. Corresponds to the rows of
optimafor which the associated row in columnAnswering_Setis equal to 1. Only applicable for minimax or standardized maximin designs.max_derivGlobal maximum of the sensitivity (derivative) function (ε in 'Details').
ELBD-efficiency lower bound. Can not be larger than 1. If negative, see 'Note' in
sensminimaxorsens.minimax.control.merge_tolMerging tolerance to create the answering set from the set of all local optima. See 'Details' in
sens.minimax.control. Only applicable for minimax or standardized maximin designs.crtvalCriterion value. Compare it with the column
Crtiterion_Valueinoptimafor minimax and standardized maximin designs.timeUsed CPU time (rough approximation).
Note
Theoretically, ELB can not be larger than 1. But if so, it may have one of the following reasons:
-
max_derivis not a GLOBAL maximum. Please increase the value of the parametermaxevalinsens.minimax.controlto find the global maximum. The sensitivity function is shifted below the y-axis because the number of model parameters has not been specified correctly (less value given). Please specify the correct number of model parameters via argument
npar.
Please increase the value of the parameter
n_seg in sens.minimax.control
for models with larger number of parameters or large parameter space to find the true
answering set for minimax and standardized maximin designs. See sens.minimax.control for more details.
References
Masoudi E, Holling H, Wong W.K. (2017). Application of Imperialist Competitive Algorithm to Find Minimax and Standardized Maximin Optimal Designs. Computational Statistics and Data Analysis, 113, 330-345.
Examples
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 | ##########################
# Power logistic model
##########################
# verifying the minimax D-optimality of a design with points x0 and weights w0
x0 <- c(-4.5515, 0.2130, 2.8075)
w0 <- c(0.4100, 0.3723, 0.2177)
# Power logistic model when s = .2
sensminimax(formula = ~ (1/(1 + exp(-b * (x-a))))^.2,
predvars = "x",
parvars = c("a", "b"),
family = binomial(),
x = x0, w = w0,
lx = -5, ux = 5,
lp = c(0, 1), up = c(3, 1.5))
##############################
# A model with two predictors
##############################
# Verifying the minimax D-optimality of a design for a model with two predictors
# The model is the mixed inhibition model.
# X0 is the vector of four design points that are:
# (3.4614, 0) (4.2801, 3.1426) (30, 0) (30, 4.0373)
x0 <- c(3.4614, 4.2801, 30, 30, 0, 3.1426, 0, 4.0373)
w0 <- rep(1/4, 4)
sensminimax(formula = ~ V*S/(Km * (1 + I/Kic)+ S * (1 + I/Kiu)),
predvars = c("S", "I"),
parvars = c("V", "Km", "Kic", "Kiu"),
family = "gaussian",
x = x0, w = w0,
lx = c(0, 0), ux = c(30, 60),
lp = c(1, 4, 2, 4), up = c(1, 5, 3, 5))
##########################################
# Standardized maximin D-optimal designs
##########################################
# Verifying the standardized maximin D-optimality of a design for
# the loglinear model
# First we should define the function for 'localdes' argument
# The function LDOD takes the parameters and returns the points and
# weights of the locally D-optimal design
LDOD <- function(theta0, theta1, theta2){
## param is the vector of theta = (theta0, theta1, theta2)
lx <- 0 # lower bound of the design space
ux <- 150 # upper bound of the design space
param <- c()
param[1] <- theta0
param[2] <- theta1
param[3] <- theta2
xstar <- (ux+param[3]) * (lx + param[3]) *
(log(ux + param[3]) - log(lx + param[3]))/(ux - lx) - param[3]
return(list(x = c(lx, xstar, ux) , w = rep(1/3, 3)))
}
x0 <- c(0, 4.2494, 17.0324, 149.9090)
w0 <- c(0.3204, 0.1207, 0.2293, 0.3296)
## Not run:
sensminimax(formula = ~theta0 + theta1* log(x + theta2),
predvars = c("x"),
parvars = c("theta0", "theta1", "theta2"),
x = x0, w = w0,
lx = 0, ux = 150,
lp = c(2, 2, 1), up = c(2, 2, 15),
localdes = LDOD,
standardized = TRUE,
sens.minimax.control = list(n_seg = 10))
## End(Not run)
################################################################
# Not necessary!
# The rest of the examples here are only for professional uses.
################################################################
# Imagine you have written your own FIM, say in Rcpp that is faster than
# the FIM created by the formula interface here.
##########################
# Power logistic model
##########################
# For example, th cpp FIM function for the power logistic model is named:
FIM_power_logistic
args(FIM_power_logistic)
# The arguments do not match the standard of the argument 'fimfunc'
# in 'sensminimax'
# So we reparameterize it:
myfim1 <- function(x, w, param)
FIM_power_logistic(x = x, w = w, param =param, s = .2)
args(myfim1)
## Not run:
# Verify minimax D-optimality of a design
sensminimax(fimfunc = myfim1,
x = c(-4.5515, 0.2130, 2.8075),
w = c(0.4100, 0.3723, 0.2177),
lx = -5, ux = 5,
lp = c(0, 1), up = c(3, 1.5))
## End(Not run)
##############################
# A model with two predictors
##############################
# An example of a model with two-predictors: mixed inhibition model
# Fisher information matrix:
FIM_mixed_inhibition
args(FIM_mixed_inhibition)
# We should first reparameterize the FIM to match the standard of the
# argument 'fimfunc'
myfim2 <- function(x, w, param){
npoint <- length(x)/2
S <- x[1:npoint]
I <- x[(npoint+1):(npoint*2)]
out <- FIM_mixed_inhibition(S = S, I = I, w = w, param = param)
return(out)
}
args(myfim2)
## Not run:
# Verifyng minimax D-optimality of a design
sensminimax(fimfunc = myfim2,
x = c(3.4614, 4.2801, 30, 30, 0, 3.1426, 0, 4.0373),
w = rep(1/4, 4),
lx = c(0, 0), ux = c(30, 60),
lp = c(1, 4, 2, 4), up = c(1, 5, 3, 5))
## End(Not run)
#########################################
# Standardized maximin D-optimal designs
#########################################
# An example of a user-written FIM function:
help(FIM_loglin)
# An example of verfying standardaized maximin D-optimality for a design
# Look how we re-define the function LDOD above
LDOD2 <- function(param){
## param is the vector of theta = (theta0, theta1, theta2)
lx <- 0 # lower bound of the design space
ux <- 150 # upper bound of the design space
xstar <- (ux + param[3]) * (lx + param[3]) *
(log(ux + param[3]) - log(lx + param[3]))/(ux - lx) - param[3]
return(list(x = c(lx, xstar, ux) , w = rep(1/3, 3)))
}
args(LDOD2)
sensminimax(fimfunc = FIM_loglin,
x = x0,
w = w0,
lx = 0, ux = 150,
lp = c(2, 2, 1), up = c(2, 2, 15),
localdes = LDOD2,
standardized = TRUE)
###################################
# user-defined optimality criterion
##################################
# When the model is defined by the formula interface
# Checking the A-optimality 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))
}
sensminimax(formula = ~1/(1 + exp(-b * (x-a))), predvars = "x",
parvars = c("a", "b"), family = "binomial",
lp = c(-2, 1), up = c(2, 1.5),
crtfunc = Aopt,
lx = -2, ux = 2,
sensfunc = Aopt_sens,
x = c(-2, .0033, 2), w = c(.274, .452, .274))
|
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.