Learn R Programming

LDOD (version 1.0)

cfisher: Auto-constructing Fisher Information matrix

Description

Auto-constructs Fisher information matrix for nonlinear and generalized linear models as two Rfunctions.

Usage

cfisher(ymean, yvar, ndpoints, prec = 53)

Arguments

ymean
a character string, formula of $E(y)$ with specific satndard: characters b1, b2, b3, ...symbolize model parameters and x1, x2, x3, ...symbolize explanatory variables. See 'Exa
yvar
a character string, formula of $Var(y)$ with specific standard as ymean. See 'Details' and 'Examples'.
ndpoints
number of design points.
prec
(optional) a number, the maximal precision to be used for D-efficiency calculation, in bite. Must be at least $2$ (default $53$).

Value

  • a list containing two closures:
  • fima function in which its arguments are vector of design points ($x$), vector of corresponding weights ($w$) and vector of parameters ($\beta$) and its output is Fisher information matrix.
  • fim.mpfra function in which its arguments are vector of design points ($x$), vector of corresponding weights ($w$) and vector of parameters ($\beta$) and its output is Fisher information matrix of class 'mpfrMatrix'.

Details

If response variables have the same constant variance, for example $\sigma^2$, then yvar must be $1$.

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.

Examples

Run this code
## Logistic dose response model
ymean <- "(1/(exp(-b2 * (x1 - b1)) + 1))"
yvar <- "(1/(exp(-b2 * (x1 - b1)) + 1)) * (1 - (1/(exp(-b2 * (x1 - b1)) + 1)))"
res <- cfisher(ymean, yvar, ndpoints = 2, prec = 54)

# res$fim is Fisher information matrix for a two-points design 
res$fim(x = c(x11 = 2, x21 = 3), w = c(w1 = .5, w2 = .5), b = c(b1 = .9, b2 = .8))

# res$fim is Fisher information matrix for a two-points design with 54 precision
res$fim.mpfr(x = c(x11 = 2, x21 = 3), w = c(w1 = .5, w2 = .5), b = c(b1 = .9, b2 = .8))

# Fisher information matrix for model:
fim<- cfisher(ymean, yvar, ndpoints = 1, prec = 54)
res$fim(x = c(x11 = 2), w = c(w1 = 1), b = c(b1 = .9, b2 = .8))

## posison with E(y) = exp(b1 + b2 * x1 + b3 * x1^2 + b4 * x2 +b5 * x2^2 + b6 * x1 * x2)
ymean <- yvar <- "exp(b1 + b2 * x1 + b3 * x1^2 + b4 * x2 +b5 * x2^2 + b6 * x1 * x2)"
fim <- cfisher(ymean, yvar, ndpoints = 6, prec = 54)

# res$fim is Fisher information matrix for a six-points design 
res$fim(x = c(1:12), w = rep(1/6, 6), b = c(1:6)) ## NAN

# res$fim.mpfr is Fisher information matrix for a six-points design with 53 precision
res$fim.mpfr(x = c(1:12), w = rep(1/6, 6), b = c(1:6))

## Linear regression with two indeoendent varibales (the design points are two-dimensional)
ymean <- "b1 + b2 * x1 + b3 * x2"
yvar = "1"
res <- cfisher(ymean, yvar, ndpoints = 3, prec = 54)
res$fim(x = c(1:6), w = c(.3, .3, .3))
res$fim.mpfr(x = c(1:6), w = c(.3, .3, .3))

## Logistic model:
ymean <- "1/(exp(-b1 - b2 * x1) + 1)"
yvar <- "(1/(exp(-b1 - b2 * x1) + 1)) * (1 - (1/(exp(-b1 - b2 * x1) + 1)))"
cfisher(ymean, yvar, ndpoints = 2, prec = 54)

## Poisson model:
ymean <- yvar <-  "exp(b1 + b2 * x1)"
cfisher(ymean, yvar, ndpoints = 2, prec = 54)

## Poisson dose response model:
ymean <- yvar <- "b1 * exp(-b2 * x1)"
cfisher(ymean, yvar, ndpoints = 2, prec = 54)

## Inverse Quadratic model:
ymean <- "(b1 * x1)/(b2 + x1 + b3 * (x1)^2)"
yvar <- "1"
cfisher(ymean, yvar, ndpoints = 3, prec = 54)
#
ymean <- "x1/(b1 + b2 * x1 + b3 * (x1)^2)"
yvar <- "1"
cfisher(ymean, yvar, ndpoints = 3, prec = 54)

## Weibull model:
ymean <- "b1 - b2 * exp(-b3 * x1^b4)"
yvar <- "1"
cfisher(ymean, yvar, ndpoints = 4, prec = 54)

## Richards model:
ymean <- "b1/(1 + b2 * exp(-b3 * x1))^b4"
yvar <- "1"
cfisher(ymean, yvar, ndpoints = 4, prec = 54)

## Michaelis-Menten model:
ymean <- "(b1 * x1)/(1 + b2 * x1)"
yvar <- "1"
cfisher(ymean, yvar, ndpoints = 2, prec = 54)
#
ymean <- "(b1 * x1)/(b2 + x1)"
yvar <- "1"
cfisher(ymean, yvar, ndpoints = 2, prec = 54)
#
ymean <- "x1/(b1 + b2 * x1)"
yvar <- "1"
cfisher(ymean, yvar, ndpoints = 2, prec = 54)

## log-linear model
ymean <- "b1 + b2 * log(x1 + b3)"
yvar <- "1"
cfisher(ymean, yvar, ndpoints = 3, prec = 54)

## Exponential model:
ymean <- "b1 + b2 * exp(x1/b3)"
yvar <- "1"
cfisher(ymean, yvar, ndpoints = 3, prec = 54)

## Emax model:
ymean  <- "b1 + (b2 * x1)/(x1 + b3)"
yvar <- "1"
cfisher(ymean, yvar, ndpoints = 3, prec = 54)

## Negative binomial model Y ~ NB(E(Y), theta) where E(Y) = b1*exp(-b2*x1):
theta = 5
ymean <- "b1 * exp(-b2 * x1)"
yvar <- paste("b1 * exp(-b2 * x1) * (1 + (1/", theta, ") * b1 * exp(-b2 * x1))", sep = "")
cfisher(ymean, yvar, ndpoints = 3, prec = 54)

Run the code above in your browser using DataLab