Finds Locally D-optimal designs for Richards regression model which is defined as E(y) = a/(1+bexp(-λ*x))^h with Var(y) = σ^2, where a, b, λ, h and σ are unknown parameters.
lower bound of design interval, must be greater than or equal to 0.
ub
upper bound of design interval.
user.points
(optional) vector of user design points which calculation of its D-efficiency is aimed. Each element of user.points must be within the design interval.
user.weights
(optional) vector of weights which its elements correspond to user.points elements. The sum of weights should be 1; otherwise they will be normalized.
...
(optional) additional parameters will be passed to function
curve.
prec
(optional)
a number, the maximal precision to be used for D-efficiency calculation, in bite. Must be at least 2 (default 53), see 'Details'.
n.restarts
(optional optimization parameter)
number of solver restarts required in optimization process (default 1), see 'Details'.
n.sim
(optional optimization parameter)
number of random parameters to generate for every restart of solver in optimization process (default 1), see 'Details'.
tol
(optional optimization parameter)
relative tolerance on feasibility and optimality in optimization process (default 1e-8).
rseed
(optional optimization parameter) a seed to initiate the random number generator, else system time will be used.
Details
While D-efficiency is NaN, an increase in prec can be beneficial to achieve a numeric value, however, it can slow down the calculation speed.
Values of n.restarts and n.sim should be chosen according to the length of design interval.
Value
plot of derivative function, see 'Note'.
a list containing the following values:
points
obtained design points
weights
corresponding weights to the obtained design points
det.value
value of Fisher information matrix determinant at the obtained design
user.eff
D-efficeincy of user design, if user.design and user.weights are not NULL.
Note
To verify optimality of obtained design, derivate function
(symmetry of Frechet derivative with respect to the x-axis)
will be plotted on the design interval. Based on the equivalence theorem (Kiefer, 1974), a design is optimal if and only if its derivative function are equal or less than 0 on the design interval. The equality must be achieved just at the obtained points.
Author(s)
Ehsan Masoudi, Majid Sarmad and Hooshang Talebi
References
Masoudi, E., Sarmad, M. and Talebi, H. 2012, An Almost General Code in R to Find Optimal Design, In Proceedings of the 1st ISM International Statistical Conference 2012, 292-297.
Dette, H., Pepelyshev, A. (2008), Efficient Experimental Designs for Sigmoidal Growth Models, Statistical Planning and Inference, 138, 2-17.
Kiefer, J. C. 1974, General equivalence theory for optimum designs (approximate theory), Ann. Statist., 2, 849-879.
See Also
cfisher, cfderiv and eff.
Examples
ldrichards(a = 1, b = 2, lambda = 2, h = 3, lb = 0, ub =3)
# $points: 0.1805017 0.8296549 1.6139494 3.0000000
## usage of n.sim and n.restars
# Various responses for different values of rseed
ldrichards(a = 1, b = 4, lambda = 3, h = 6, lb = 0, ub = 19, rseed = 6)
# $points: 5.022689 11.520735 17.815197 19.000000
ldrichards(a = 1, b = 4, lambda = 3, h = 6, lb = 0, ub = 19, rseed = 7)
# $points: 2.198258 7.557164 18.789277 19.000000
ldrichards(a = 1, b = 4, lambda = 3, h = 6, lb = 0, ub = 19, n.sim = 5, n.restarts = 5)
# (valid response) $points: 0.6562008 1.0485843 1.5894946 19.000000
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(LDOD)
Loading required package: Rsolnp
Loading required package: Rmpfr
Loading required package: gmp
Attaching package: 'gmp'
The following objects are masked from 'package:base':
%*%, apply, crossprod, matrix, tcrossprod
C code of R package 'Rmpfr': GMP using 64 bits per limb
Attaching package: 'Rmpfr'
The following objects are masked from 'package:stats':
dbinom, dnorm, dpois, pnorm
The following objects are masked from 'package:base':
cbind, pmax, pmin, rbind
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/LDOD/ldrichards.Rd_%03d_medium.png", width=480, height=480)
> ### Name: ldrichards
> ### Title: Locally D-optimal designs for Richards model
> ### Aliases: ldrichards
> ### Keywords: optimal design Richards equivalence theorem
>
> ### ** Examples
>
> ldrichards(a = 1, b = 2, lambda = 2, h = 3, lb = 0, ub =3)
Iter: 1 fn: 22.9881 Pars: 0.82965 1.61394 0.18050
Iter: 2 fn: 22.9881 Pars: 0.82966 1.61395 0.18050
solnp--> Completed in 2 iterations
$points
[1] 0.1805017 0.8296555 1.6139499 3.0000000
$weights
[1] 0.25 0.25 0.25 0.25
$det.value
[1] 1.038522e-10
> # $points: 0.1805017 0.8296549 1.6139494 3.0000000
>
> ## usage of n.sim and n.restars
> # Various responses for different values of rseed
>
> ldrichards(a = 1, b = 4, lambda = 3, h = 6, lb = 0, ub = 19, rseed = 6)
Iter: 1 fn: 152.0706 Pars: 11.51924 17.81520 5.02240
Iter: 2 fn: 152.0706 Pars: 11.51924 17.81520 5.02240
solnp--> Completed in 2 iterations
$points
[1] 5.022404 11.519241 17.815197 19.000000
$weights
[1] 0.25 0.25 0.25 0.25
$det.value
[1] 9.04825e-67
Warning message:
inverse of information matrix can not be calculated because of computational singularity at obtained design
> # $points: 5.022689 11.520735 17.815197 19.000000
>
> ldrichards(a = 1, b = 4, lambda = 3, h = 6, lb = 0, ub = 19, rseed = 7)
Iter: 1 fn: 84.1221 Pars: 18.78928 5.83263 1.84200
Iter: 2 fn: 76.9919 Pars: 18.78928 5.83263 1.24473
Iter: 3 fn: 76.9919 Pars: 18.78928 5.83263 1.24473
solnp--> Completed in 3 iterations
$points
[1] 1.244733 5.832633 18.789277 19.000000
$weights
[1] 0.25 0.25 0.25 0.25
$det.value
[1] 3.654677e-34
Warning message:
inverse of information matrix can not be calculated because of computational singularity at obtained design
> # $points: 2.198258 7.557164 18.789277 19.000000
>
> ldrichards(a = 1, b = 4, lambda = 3, h = 6, lb = 0, ub = 19, n.sim = 5, n.restarts = 5)
Iter: 1 fn: 27.1504 Pars: 1.04858 0.65620 1.58950
Iter: 2 fn: 27.1504 Pars: 1.04858 0.65620 1.58949
solnp--> Completed in 2 iterations
Iter: 1 fn: 65.0738 Pars: 1.00674 15.42080 3.74223
solnp--> Completed in 1 iterations
Iter: 1 fn: 59.9396 Pars: 7.73439 0.28567 1.37634
Iter: 2 fn: 59.9396 Pars: 7.73439 0.28567 1.37634
solnp--> Completed in 2 iterations
Iter: 1 fn: 58.8803 Pars: 5.09724 0.84411 3.85059
Iter: 2 fn: 58.8803 Pars: 5.09724 0.84411 3.85059
solnp--> Completed in 2 iterations
Iter: 1 fn: 68.4498 Pars: 2.59760 1.95802 9.89127
Iter: 2 fn: 67.8799 Pars: 2.59760 1.95802 9.89127
Iter: 3 fn: 61.6666 Pars: 2.59760 1.43278 9.89127
Iter: 4 fn: 61.4577 Pars: 2.57706 1.41454 9.89127
Iter: 5 fn: 61.4577 Pars: 2.57706 1.41454 9.89127
solnp--> Completed in 5 iterations
$points
[1] 0.6562008 1.0485843 1.5894949 19.0000000
$weights
[1] 0.25 0.25 0.25 0.25
$det.value
[1] 1.617045e-12
> # (valid response) $points: 0.6562008 1.0485843 1.5894946 19.000000
>
>
>
>
>
> dev.off()
null device
1
>