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("statmod-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('statmod')
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "digammaf"
> 
> ### * digammaf
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: Digamma
> ### Title: Digamma generalized linear model family
> ### Aliases: Digamma canonic.digamma d2cumulant.digamma
> ###   unitdeviance.digamma cumulant.digamma meanval.digamma varfun.digamma
> ### Keywords: models
> 
> ### ** Examples
> 
> # Test for log-linear dispersion trend in gamma regression
> y <- rchisq(20,df=1)
> x <- 1:20
> out.gam <- glm(y~x,family=Gamma(link="log"))
> d <- residuals(out.gam)^2
> out.dig <- glm(d~x,family=Digamma(link="log"))
> summary(out.dig,dispersion=2)

Call:
glm(formula = d ~ x, family = Digamma(link = "log"))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.4971  -1.1866  -0.5492   0.1782   1.8488  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  0.89184    0.63193   1.411    0.158
x           -0.02461    0.05302  -0.464    0.642

(Dispersion parameter for Digamma family taken to be 2)

    Null deviance: 36.144  on 19  degrees of freedom
Residual deviance: 35.766  on 18  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 8

> 
> 
> 
> cleanEx(); ..nameEx <- "gauss.quad"
> 
> ### * gauss.quad
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: gauss.quad
> ### Title: Gaussian Quadrature
> ### Aliases: gauss.quad
> ### Keywords: math
> 
> ### ** Examples
> 
> out <- gauss.quad(10,"laguerre",alpha=5)
> sum(out$weights * out$nodes) / gamma(6)
[1] 6
> #  mean of gamma distribution with alpha=6
> 
> 
> 
> cleanEx(); ..nameEx <- "gauss.quad.prob"
> 
> ### * gauss.quad.prob
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: gauss.quad.prob
> ### Title: Gaussian Quadrature with Probability Distributions
> ### Aliases: gauss.quad.prob
> ### Keywords: math
> 
> ### ** Examples
> 
> out <- gauss.quad.prob(10,"normal")
> sum(out$weights * out$nodes^4)
[1] 3
> #  the 4th moment of the standard normal is 3
> 
> out <- gauss.quad.prob(32,"gamma",alpha=5)
> sum(out$weights * log(out$nodes))
[1] 1.506118
> #  the expected value of log(X) where X is gamma is digamma(alpha)
> 
> 
> 
> cleanEx(); ..nameEx <- "glmgam"
> 
> ### * glmgam
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: glmgam.fit
> ### Title: Gamma Generalized Linear Model with Identity Link
> ### Aliases: glmgam.fit
> ### Keywords: regression
> 
> ### ** Examples
> 
> y <- rgamma(10,shape=5)
> X <- cbind(1,1:10)
> fit <- glmgam.fit(X,y,trace=TRUE)
Iter = 0 , Dev = 1.561868  Beta 6.127364 -0.1990615 
Iter = 1 , Dev = 1.560572  Beta 6.129285 -0.2057179 
Iter = 2 , Dev = 1.559795  Beta 6.154427 -0.211027 
Iter = 3 , Dev = 1.557784  Beta 6.260201 -0.2254861 
Iter = 4 , Dev = 1.557186  Beta 6.340586 -0.2367725 
Iter = 5 , Dev = 1.557158  Beta 6.357378 -0.239427 
Iter = 6 , Dev = 1.557157  Beta 6.360288 -0.2399164 
> 
> 
> 
> cleanEx(); ..nameEx <- "growthcurve"
> 
> ### * growthcurve
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: growthcurve
> ### Title: Compare Groups of Growth Curves
> ### Aliases: compareGrowthCurves compareTwoGrowthCurves
> ### Keywords: regression
> 
> ### ** Examples
> 
> # A example with only one time
> data(PlantGrowth)
> compareGrowthCurves(PlantGrowth$group,as.matrix(PlantGrowth$weight))
ctrl trt1  1.19 
ctrl trt2  -2.13 
trt1 trt2  -3.01 
  Group1 Group2      Stat P.Value adj.P.Value
1   ctrl   trt1  1.191260    0.27        0.27
2   ctrl   trt2 -2.134020    0.05        0.10
3   trt1   trt2 -3.010099    0.00        0.00
> # Can make p-values more accurate by nsim=10000
> 
> 
> 
> cleanEx(); ..nameEx <- "hommel.test"
> 
> ### * hommel.test
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: hommel.test
> ### Title: Test Multiple Comparisons Using Hommel's Method
> ### Aliases: hommel.test
> ### Keywords: htest
> 
> ### ** Examples
> 
> p <- sort(runif(100))[1:10]
> cbind(p,p.adjust(p,"hommel"),hommel.test(p))
               p            
 [1,] 0.01339033 0.1205130 1
 [2,] 0.02333120 0.1255551 1
 [3,] 0.05893438 0.1255551 1
 [4,] 0.06178627 0.1255551 1
 [5,] 0.07067905 0.1255551 1
 [6,] 0.08424691 0.1255551 1
 [7,] 0.09946616 0.1255551 1
 [8,] 0.10794363 0.1255551 1
 [9,] 0.12169192 0.1255551 1
[10,] 0.12555510 0.1255551 1
> 
> 
> 
> cleanEx(); ..nameEx <- "invgauss"
> 
> ### * invgauss
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: invgauss
> ### Title: Inverse Gaussian Distribution
> ### Aliases: InverseGaussian dinvgauss pinvgauss qinvgauss rinvgauss
> ### Keywords: distribution
> 
> ### ** Examples
> 
> y <- rinvgauss(20,1,2) # generate vector of 20 random numbers
> p <- pinvgauss(y,1,2) # p should be uniform
> 
> 
> 
> cleanEx(); ..nameEx <- "limdil"
> 
> ### * limdil
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: limdil
> ### Title: Limiting Dilution Analysis
> ### Aliases: limdil
> ### Keywords: regression
> 
> ### ** Examples
> 
> Dose <- c(50,100,200,400,800)
> Responses <- c(2,6,9,15,21)
> Tested <- c(24,24,24,24,24)
> limdil(Responses,Dose,Tested,test.unit.slope=TRUE)
$CI
   Lower Estimate    Upper 
537.7392 403.0632 302.1166 

$test.unit.slope
    Chisq   P.value 
0.0666993 0.7962047 

> 
> 
> 
> cleanEx(); ..nameEx <- "logmdigamma"
> 
> ### * logmdigamma
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: logmdigamma
> ### Title: Log Minus Digamma Function
> ### Aliases: logmdigamma
> ### Keywords: math
> 
> ### ** Examples
> 
> log(10^15) - digamma(10^15) # returns 0
[1] 0
> logmdigamma(10^15) # returns value correct to 15 figures
[1] 5e-16
> 
> 
> 
> cleanEx(); ..nameEx <- "matvec"
> 
> ### * matvec
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: matvec
> ### Title: Multiply a Matrix by a Vector
> ### Aliases: matvec vecmat
> ### Keywords: array algebra
> 
> ### ** Examples
> 
> A <- matrix(1:12,3,4)
> A
     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12
> matvec(A,c(1,2,3,4))
     [,1] [,2] [,3] [,4]
[1,]    1    8   21   40
[2,]    2   10   24   44
[3,]    3   12   27   48
> vecmat(c(1,2,3),A)
     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    4   10   16   22
[3,]    9   18   27   36
> 
> 
> 
> cleanEx(); ..nameEx <- "meanT"
> 
> ### * meanT
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: meanT
> ### Title: Mean t-Statistic Between Two Groups of Growth Curves
> ### Aliases: meanT
> ### Keywords: regression
> 
> ### ** Examples
> 
> y1 <- matrix(rnorm(4*3),4,3)
> y2 <- matrix(rnorm(4*3),4,3)
> meanT(y1,y2)
[1] 0.1232123
> 
> data(PlantGrowth)
> compareGrowthCurves(PlantGrowth$group,as.matrix(PlantGrowth$weight))
ctrl trt1  1.19 
ctrl trt2  -2.13 
trt1 trt2  -3.01 
  Group1 Group2      Stat P.Value adj.P.Value
1   ctrl   trt1  1.191260    0.17        0.17
2   ctrl   trt2 -2.134020    0.06        0.12
3   trt1   trt2 -3.010099    0.00        0.00
> # Can make p-values more accurate by nsim=10000
> 
> 
> 
> cleanEx(); ..nameEx <- "power"
> 
> ### * power
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: power.fisher.test
> ### Title: Power of Fisher's Exact Test for Comparing Proportions
> ### Aliases: power.fisher.test
> ### Keywords: htest
> 
> ### ** Examples
> 
> power.fisher.test(0.5,0.9,20,20) # 70
[1] 0.75
> 
> 
> 
> cleanEx(); ..nameEx <- "qresiduals"
> 
> ### * qresiduals
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: qresiduals
> ### Title: Randomized Quantile Residuals
> ### Aliases: qresiduals qresid qres.binom qres.pois qres.nbinom qres.gamma
> ###   qres.invgauss qres.tweedie qres.default
> ### Keywords: regression
> 
> ### ** Examples
> 
> #  Poisson example: quantile residuals show no granularity
> y <- rpois(20,lambda=4)
> x <- 1:20
> fit <- glm(y~x, family=poisson)
> qr <- qresiduals(fit)
> qqnorm(qr)
> abline(0,1)
> 
> #  Gamma example:
> #  Quantile residuals are nearly normal while usual resids are not
> y <- rchisq(20, df=1)
> fit <- glm(y~1, family=Gamma)
> qr <- qresiduals(fit, dispersion=2)
> qqnorm(qr)
> abline(0,1)
> 
> 
> 
> cleanEx(); ..nameEx <- "ranblock"
> 
> ### * ranblock
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: ranblock
> ### Title: Randomized Block Mixed Linear Model
> ### Aliases: randomizedBlock randomizedBlockFit
> ### Keywords: regression
> 
> ### ** Examples
> 
> #  Compare with first data example from Venable and Ripley (2002),
> #  Chapter 10, "Linear Models"
> library(MASS)
> data(petrol)
> out <- randomizedBlock(Y~SG+VP+V10+EP, random=No, data=petrol)
> cbind(varcomp=out$varcomp,se=out$se.varcomp)
          varcomp       se
Residual 3.505151 1.081110
Block    2.087317 1.908131
> 
> 
> 
> cleanEx(); ..nameEx <- "remlscor"
> 
> ### * remlscor
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: remlscore
> ### Title: REML for Heteroscedastic Regression
> ### Aliases: remlscore
> ### Keywords: regression
> 
> ### ** Examples
> 
> data(welding)
> attach(welding)
> y <- Strength
> # Reproduce results from Table 1 of Smyth (2002)
> X <- cbind(1,(Drying+1)/2,(Material+1)/2)
> colnames(X) <- c("1","B","C")
> Z <- cbind(1,(Material+1)/2,(Method+1)/2,(Preheating+1)/2)
> colnames(Z) <- c("1","C","H","I")
> out <- remlscore(y,X,Z)
> cbind(Estimate=out$gamma,SE=out$se.gam)
                        SE
[1,] -3.15886017 0.8313270
[2,] -2.73542576 0.8224828
[3,] -0.08588713 0.8351308
[4,]  3.33238821 0.8250499
> 
> 
> 
> cleanEx(); ..nameEx <- "remlscorgamma"
> 
> ### * remlscorgamma
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: remlscoregamma
> ### Title: Approximate REML for gamma regression with structured dispersion
> ### Aliases: remlscoregamma
> ### Keywords: regression
> 
> ### ** Examples
> 
> data(welding)
> attach(welding)
> y <- Strength
> X <- cbind(1,(Drying+1)/2,(Material+1)/2)
> colnames(X) <- c("1","B","C")
> Z <- cbind(1,(Material+1)/2,(Method+1)/2,(Preheating+1)/2)
> colnames(Z) <- c("1","C","H","I")
> out <- remlscoregamma(y,X,Z)
> 
> 
> 
> cleanEx(); ..nameEx <- "sage.test"
> 
> ### * sage.test
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: sage.test
> ### Title: Compare Two SAGE Libraries
> ### Aliases: sage.test
> ### Keywords: htest
> 
> ### ** Examples
> 
> sage.test(c(0,5,10),c(0,30,50),n1=10000,n2=15000)
[1] 1.0000000000 0.0015428182 0.0001689180
> #  Exact equivalents
> fisher.test(matrix(c(0,0,10000-0,15000-0),2,2))$p.value
[1] 1
> fisher.test(matrix(c(5,30,10000-5,15000-30),2,2))$p.value
[1] 0.001531718
> fisher.test(matrix(c(10,50,10000-10,15000-50),2,2))$p.value
[1] 0.0001660018
> 
> 
> 
> cleanEx(); ..nameEx <- "tweedie"
> 
> ### * tweedie
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: tweedie
> ### Title: Tweedie Generalized Linear Models
> ### Aliases: tweedie
> ### Keywords: regression
> 
> ### ** Examples
> 
> y <- rgamma(20,shape=5)
> x <- 1:20
> # Fit a poisson generalized linear model with identity link
> glm(y~x,family=tweedie(var.power=1,link.power=1))

Call:  glm(formula = y ~ x, family = tweedie(var.power = 1, link.power = 1)) 

Coefficients:
(Intercept)            x  
    5.35647     -0.05365  

Degrees of Freedom: 19 Total (i.e. Null);  18 Residual
Null Deviance:	    14.88 
Residual Deviance: 14.47 	AIC: NA 
> 
> # Fit an inverse-Gaussion glm with log-link
> glm(y~x,family=tweedie(var.power=3,link.power=0)) 

Call:  glm(formula = y ~ x, family = tweedie(var.power = 3, link.power = 0)) 

Coefficients:
(Intercept)            x  
    1.68057     -0.01100  

Degrees of Freedom: 19 Total (i.e. Null);  18 Residual
Null Deviance:	    1.017 
Residual Deviance: 0.9993 	AIC: NA 
> 
> 
> 
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  2.22 0.1 2.33 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')