R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.1.1  (2005-06-20), ISBN 3-900051-07-0

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 a HTML browser interface to help.
Type 'q()' to quit R.

> ### * <HEADER>
> ###
> attach(NULL, name = "CheckExEnv")
> assign(".CheckExEnv", as.environment(2), pos = length(search())) # base
> ## add some hooks to label plot pages for base and grid graphics
> setHook("plot.new", ".newplot.hook")
> setHook("persp", ".newplot.hook")
> setHook("grid.newpage", ".gridplot.hook")
> 
> assign("cleanEx",
+        function(env = .GlobalEnv) {
+ 	   rm(list = ls(envir = env, all.names = TRUE), envir = env)
+            RNGkind("default", "default")
+ 	   set.seed(1)
+    	   options(warn = 1)
+ 	   delayedAssign("T", stop("T used instead of TRUE"),
+ 		  assign.env = .CheckExEnv)
+ 	   delayedAssign("F", stop("F used instead of FALSE"),
+ 		  assign.env = .CheckExEnv)
+ 	   sch <- search()
+ 	   newitems <- sch[! sch %in% .oldSearch]
+ 	   for(item in rev(newitems))
+                eval(substitute(detach(item), list(item=item)))
+ 	   missitems <- .oldSearch[! .oldSearch %in% sch]
+ 	   if(length(missitems))
+ 	       warning("items ", paste(missitems, collapse=", "),
+ 		       " have been removed from the search path")
+        },
+        env = .CheckExEnv)
> assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now
> assign("ptime", proc.time(), env = .CheckExEnv)
> grDevices::postscript("asypow-Examples.ps")
> assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv)
> options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly"))
> options(warn = 1)    
> library('asypow')
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "asypow.n"
> 
> ### * asypow.n
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: asypow.n
> ### Title: Asymptotic Sample Size
> ### Aliases: asypow.n
> ### Keywords: htest
> 
> ### ** Examples
> 
> # Three Sample Poisson Example :
> #    Three independent Poisson processes produce events at
> #    mean rates of 1, 2 and 3 per day. For how many days
> #    must the processes be observed to have an 80% chance
> #    of detecting that the means are different at an
> #    overall significance level of 0.05?
> # Step 1 : Find the information matrix
> pois.mean <- c(1,2,3)
> info.pois <- info.poisson.kgroup(pois.mean, group.size=3)
> # Step 2: Create the constraints matrix
> constraints <- matrix(c(2,1,2,2,2,3),ncol=3,byrow=TRUE)
> # Step 3: Find the noncentrality parameter and
> #         degrees of freedom for the test
> poisson.object <- asypow.noncent(pois.mean, info.pois, constraints)
> # Step 4: Compute sample size needed for 
> #         0.8 power with significance level 0.05
> n.pois <- asypow.n(poisson.object, 0.8, 0.05)
> # Step 5: Divide the sample size by 3 (the number of processes)
> #         to get the number of days required.
> n.days <- n.pois/3
> print(n.days)
         [,1]
[1,] 8.831797
> 
> 
> 
> cleanEx(); ..nameEx <- "asypow.noncent"
> 
> ### * asypow.noncent
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: asypow.noncent
> ### Title: Asymptotic Noncentrality Parameter
> ### Aliases: asypow.noncent
> ### Keywords: htest
> 
> ### ** Examples
> 
> # Three Sample Poisson Example :
> #    Three independent Poisson processes produce events at
> #    mean rates of 1, 2 and 3 per day. 
> # Find the information matrix
> pois.mean <- c(1,2,3)
> info.pois <- info.poisson.kgroup(pois.mean,group.size=3)
> # Create the constraints matrix
> constraints <- matrix(c(2,1,2,2,2,3),ncol=3,byrow=TRUE)
> # Calculate noncentrality parameter, degrees of freedom and parameter
> # value estimates under the null hypothesis for the test.
> poisson.object <- asypow.noncent(pois.mean,info.pois,constraints)
> 
> 
> 
> cleanEx(); ..nameEx <- "asypow.power"
> 
> ### * asypow.power
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: asypow.power
> ### Title: Asymptotic Power
> ### Aliases: asypow.power
> ### Keywords: htest
> 
> ### ** Examples
> 
> # Single Group Binomial Example:
> #    Consider testing the null hypothesis that the binomial
> #    probability is p = .2 with a sample size of 47 and
> #    signficance level of 0.05. What is the power of the
> #    test if p is actually .4?
> # Step 1: Find the information matrix
> info.binom <- info.binomial.kgroup(.4)
> # Step 2: Create the constraints matrix
> constraints <- c(1, 1, .2)
> # Step 3: Find the noncentrality parameter and
> #         degrees of freedom for the test
> binom.object <- asypow.noncent(.4, info.binom, constraints)
> # Step 4: Compute the power of a test with
> #         sample size of 47 and a significance level 0.05
> power.binom <- asypow.power(binom.object, 47, 0.05)
> print(power.binom)
[1] 0.799223
> 
> 
> 
> cleanEx(); ..nameEx <- "asypow.sig"
> 
> ### * asypow.sig
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: asypow.sig
> ### Title: Asymptotic Significance
> ### Aliases: asypow.sig
> ### Keywords: htest
> 
> ### ** Examples
> 
> # Single Group Binomial Example:
> #    Consider testing the null hypothesis that the binomial
> #    probability is p = .2 when the actual probability is .4.
> #    What significance level corresponding to a sample
> #    size of 47 and power of .8?
> # Step 1: Find the information matrix
> info.binom <- info.binomial.kgroup(.4)
> # Step 2: Create the constraints matrix
> constraints <- c(1, 1, .2)
> # Step 3: Find the noncentrality parameter and
> #         degrees of freedom for the test
> binom.object <- asypow.noncent(.4, info.binom, constraints)
> # Step 4: Compute the power of a test with
> #         sample size of 47 and a significance level 0.05
> sig.binom <- asypow.sig(binom.object, 47, 0.8)
> print(sig.binom)
[1] 0.05032492
> 
> 
> 
> cleanEx(); ..nameEx <- "info.binomial.design"
> 
> ### * info.binomial.design
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: info.binomial.design
> ### Title: Expected Information Matrix for a Binomial Design
> ### Aliases: info.binomial.design
> ### Keywords: htest
> 
> ### ** Examples
> 
> # Find the information matrix for a 2 group
> # logistic binomial Design with a quadratic
> # combination of covariate x and different
> # sample sizes at each point
> abc <- rbind(c(1.2, .9, .3),c(0.33, .21, .05))
> covar <- c(1, 2, 3,  4,5)
> sample.size <- rbind(c(10,11,12,10,9), c(8,7,10,8,9))
> info.binom <- info.binomial.design(model="quadratic", link="logistic", 
+                                 theta = abc, xpoints = covar,
+                                 natx=sample.size)
> print(info.binom)
            [,1]       [,2]       [,3]      [,4]      [,5]      [,6]
[1,] 0.009996141 0.01206818 0.01660400 0.0000000 0.0000000 0.0000000
[2,] 0.012068181 0.01660400 0.02689832 0.0000000 0.0000000 0.0000000
[3,] 0.016603997 0.02689832 0.05134918 0.0000000 0.0000000 0.0000000
[4,] 0.000000000 0.00000000 0.00000000 0.0664285 0.1664047 0.5268157
[5,] 0.000000000 0.00000000 0.00000000 0.1664047 0.5268157 1.9279663
[6,] 0.000000000 0.00000000 0.00000000 0.5268157 1.9279663 7.7111671
> 
> 
> 
> cleanEx(); ..nameEx <- "info.binomial.kgroup"
> 
> ### * info.binomial.kgroup
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: info.binomial.kgroup
> ### Title: Expected Information Matrix for Single or Multiple Group
> ###   Binomial
> ### Aliases: info.binomial.kgroup
> ### Keywords: htest
> 
> ### ** Examples
> 
> # Find the information matrix for a 2 sample binomial with
> # probability of events .2 and .4 and sample sizes 10 and 11
> info.binom <- info.binomial.kgroup(c(.2,.4), c(10,11))
> print(info.binom)
         [,1]     [,2]
[1,] 2.976190 0.000000
[2,] 0.000000 2.182540
> 
> 
> 
> cleanEx(); ..nameEx <- "info.expsurv.design"
> 
> ### * info.expsurv.design
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: info.expsurv.design
> ### Title: Expected Information Matrix for a Clinical Trial with
> ###   Exponential Survival Design
> ### Aliases: info.expsurv.design
> ### Keywords: htest
> 
> ### ** Examples
> 
> # Find the information matrix for a clinical trial
> # with hazard w(x) = -0.848 + 0.7*x which lasts
> # three years and has 10 x values equally spaced
> # between -3 and 3 with equal sample sizes.
> ab <- c(-.848, .7)
> covar <- seq(-3, 3, length=10)
> LL <- 3
> info.expsurv <- info.expsurv.design(theta = ab, L = LL, xpoints = covar)
> print(info.expsurv)
          [,1]      [,2]
[1,] 0.4661670 0.5601766
[2,] 0.5601766 1.7666561
> 
> 
> 
> cleanEx(); ..nameEx <- "info.expsurv.kgroup"
> 
> ### * info.expsurv.kgroup
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: info.expsurv.kgroup
> ### Title: Expected Information Matrix for a Single or Multiple Group
> ###   Clinical Trial with Exponential Survival
> ### Aliases: info.expsurv.kgroup
> ### Keywords: htest
> 
> ### ** Examples
> 
> # Find the information matrix for a clinical trial of
> # length 3 with hazard 1
> info.expsurv <- info.expsurv.kgroup(1, 3)
> print(info.expsurv)
[1] 0.6832624
> 
> 
> 
> cleanEx(); ..nameEx <- "info.mvlogistic"
> 
> ### * info.mvlogistic
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: info.mvlogistic
> ### Title: Expected Information Matrix for a Multivariate Logistic Model
> ### Aliases: info.mvlogistic
> ### Keywords: htest
> 
> ### ** Examples
> 
> # Find the information matrix for a multivariate
> # logistic design with variables x, y and z
> #    Define coefficient matrix so that 
> #    u = 1 + .5*x + .7*y + .9*z
> coef <- c(1, .5, .7, .9)
> #    Define the design matrix so that there are 10 design points
> intercept <- rep(1, 10)
> x <- rnorm(10)
> y <- rnorm(10)
> z <- rnorm(10)
> design <- cbind(intercept, x, y, z)
> #    Use info.mvlogistic to find the information matrix for
> #    this design
> info.xyz <- info.mvlogistic(coef, design)
> print(info.xyz)
            intercept           x           y           z
intercept  0.15628912  0.01354240  0.01456889 -0.04739962
x          0.01354240  0.09660460 -0.02036044 -0.07892256
y          0.01456889 -0.02036044  0.14113215  0.04648130
z         -0.04739962 -0.07892256  0.04648130  0.12620910
> 
> 
> 
> cleanEx(); ..nameEx <- "info.mvloglin"
> 
> ### * info.mvloglin
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: info.mvloglin
> ### Title: Expected Information Matrix for a Multivariate Log-Linear Model
> ### Aliases: info.mvloglin
> ### Keywords: htest
> 
> ### ** Examples
> 
> # Find the information matrix for a multivariate
> # log-linear design with variables x, y and z
> #    Define coefficient matrix so that 
> #    u = .1 + .2*x + .3*y + .3*z
> coef <- c(.1, .2, .3, .4)
> #    Define the design matrix so that there are 10 design points
> intercept <- rep(1, 10)
> x <- seq(.1, .2, length=10)
> y <- seq(.25, .3, length=10)
> z <- seq(.2, .3, length=10)
> design <- cbind(intercept, x, y, z)
> #    Use info.mvloglin to find the information matrix for
> #    this design
> info.xyz <- info.mvloglin(coef, design)
> print(info.xyz)
          intercept          x          y          z
intercept 4.7237667 0.34641982 0.43039106 0.29130408
x         0.3464198 0.02659991 0.03196129 0.02196045
y         0.4303911 0.03196129 0.03934650 0.02674042
z         0.2913041 0.02196045 0.02674042 0.01826283
> 
> 
> 
> cleanEx(); ..nameEx <- "info.ordinal.design"
> 
> ### * info.ordinal.design
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: info.ordinal.design
> ### Title: Expected Information Matrix for an Ordinal Design
> ### Aliases: info.ordinal.design
> ### Keywords: htest
> 
> ### ** Examples
> 
> # Find the information matrix for an ordinal design
> # with one group and equal sample sizes.
> # Assume 4 categories and use a logistic
> # line and quadratic model. Let
> # u[1] = 1 + 2.5*x
> # u[2] = 2 + 2.5*x
> # u[3] = 3 + 2.5*x
> # Use values x = -3,0,3
> theta <- c(1, 2, 3, 2.5)
> covar <- c(-3, 0, 3)
> info.ord <- info.ordinal.design(theta = theta, xpoints = covar)
> print(info.ord)
              [,1]          [,2]        [,3]          [,4]
[1,]  0.1045761734 -4.678128e-02  0.00000000  1.972695e-04
[2,] -0.0467812776  7.871164e-02 -0.02416497 -3.843805e-05
[3,]  0.0000000000 -2.416497e-02  0.02956737 -1.082200e-02
[4,]  0.0001972695 -3.843805e-05 -0.01082200  3.320982e-02
> 
> 
> 
> cleanEx(); ..nameEx <- "info.ordinal.kgroup"
> 
> ### * info.ordinal.kgroup
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: info.ordinal.kgroup
> ### Title: Expected Information Matrix for Single or Multiple Group Ordinal
> ###   Observations
> ### Aliases: info.ordinal.kgroup
> ### Keywords: htest
> 
> ### ** Examples
> 
> # Find the information matrix for a 2 group ordinal
> # model with 4 categories.
> p1 <- c(.1, .2, .3)  # Probabilities for group 1
> p2 <- c(.2, .5, .7)  # Probabilities for group 2
> p <- rbind(p1,p2)
> ngrps <- c(.4, .6)  # Percentage of data in each group
> info.ord <- info.ordinal.kgroup(p, ngrps)
> print(info.ord)
     [,1] [,2]      [,3] [,4] [,5] [,6]
[1,]    8   -4  0.000000    0    0    0
[2,]   -4    8 -4.000000    0    0    0
[3,]    0   -4  4.571429    0    0    0
[4,]    0    0  0.000000    5   -2    0
[5,]    0    0  0.000000   -2    5   -3
[6,]    0    0  0.000000    0   -3    5
> 
> 
> 
> cleanEx(); ..nameEx <- "info.poisson.design"
> 
> ### * info.poisson.design
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: info.poisson.design
> ### Title: Expected Information Matrix for a Poisson Design
> ### Aliases: info.poisson.design
> ### Keywords: htest
> 
> ### ** Examples
> 
> # Find the information matrix for a 2 group
> # logistic Poisson design with a quadratic
> # combination of covariate x and different
> # sample sizes at each point
> abc <- rbind(c(1.2,.9,.3), c(0.33,.21,.05))
> covar <- c(1, 2, 3, 4, 5)
> sample.size <- rbind(c(10,11,12,10,9), c(8,7,10,8,9))
> info.poiss <- info.poisson.design(model="quadratic",
+                                 theta = abc, xpoints = covar,
+                                 natx=sample.size)
> print(info.poiss)
           [,1]      [,2]     [,3]      [,4]      [,5]       [,6]
[1,]   53410.61  265266.5  1319435  0.000000   0.00000    0.00000
[2,]  265266.53 1319435.2  6570257  0.000000   0.00000    0.00000
[3,] 1319435.22 6570257.1 32745500  0.000000   0.00000    0.00000
[4,]       0.00       0.0        0  2.720434  10.92864   47.81822
[5,]       0.00       0.0        0 10.928636  47.81822  218.55869
[6,]       0.00       0.0        0 47.818219 218.55869 1024.97081
> 
> 
> 
> cleanEx(); ..nameEx <- "info.poisson.kgroup"
> 
> ### * info.poisson.kgroup
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: info.poisson.kgroup
> ### Title: Expected Information Matrix for Single or Multiple Group Poisson
> ### Aliases: info.poisson.kgroup
> ### Keywords: htest
> 
> ### ** Examples
> 
> # Find the information matrix for a 3 sample Poisson with
> # means 1, 2 and 3 and equal sample sizes
> info.pois <- info.poisson.kgroup(c(1,2,3))
> print(info.pois)
          [,1]      [,2]      [,3]
[1,] 0.3333333 0.0000000 0.0000000
[2,] 0.0000000 0.1666667 0.0000000
[3,] 0.0000000 0.0000000 0.1111111
> 
> 
> 
> cleanEx(); ..nameEx <- "info.reparam"
> 
> ### * info.reparam
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: info.reparam
> ### Title: Reparameterize Expected Information Matrix
> ### Aliases: info.reparam
> ### Keywords: htest
> 
> ### ** Examples
> 
> # A logistic model posits that the probability of response
> # is a logtistic function of a + b*x.
> # Consider the value of x that produces 50%
> # response, x = -a/b. Since -a/b is not one of the parameters
> # of the model, we must reparameterize to
> # roe[1] = -a/b
> # roe[2] = b
> dg <- function(theta) {
+ # theta is a vector of length 2 containing c(a,b)
+ # dg <- [d{roe[1]}/d{a} d{roe[1]}/d{b}
+ #        d{roe[2]}/d{a} d{roe[2]}/d{b}]
+   a <- theta[1]
+   b <- theta[2]
+   return(matrix(c(-1/b,a/b^2,0,1), nrow=2, ncol=2, byrow=TRUE))
+ }
> # Let a = -0.9 and b = .7
> theta <- c(-.9, .7)
> # assign a set of covariate values
> covar <- c(0.3, .9, 1.3, 2.5)
> # Use info.binomial.design to calculate the information
> # matrix under the original parameterization
> info.orig <- info.binomial.design(model="linear", link="logistic", 
+                                   theta=theta, xpoints=covar)
> # Get the information matrix of the reparameterized model
> info.new <- info.reparam(theta, info.orig, dg)
> print(info.new)
            [,1]        [,2]
[1,] 0.113645711 0.009742927
[2,] 0.009742927 0.140506899
> 
> 
> 
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  0.12 0.02 0.14 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')