.packageName <- "statmod"
#  SPECIAL FUNCTIONS

logmdigamma <- function(x)
{
#  log(x) - digamma(x)
#  Saves computation of log(x) and avoids subtractive cancellation in digamma(x) when x is large
#  Gordon Smyth, smyth@wehi.edu.au
#  19 Jan 98.  Last revised 9 Dec 2002.
#
	z <- x
	if(any(omit <- is.na(z) | Re(z) <= 0)) {
		ps <- z
		ps[omit] <- NA
		if(any(!omit)) ps[!omit] <- Recall(z[!omit])
		return(ps)
	}
	if(any(small <- Mod(z) < 5)) {
		ps <- z
		x <- z[small]
		ps[small] <- log(x/(x+5)) + Recall(x+5) + 1/x + 1/(x+1) + 1/(x+2) + 1/(x+3) + 1/(x+4)
		if(any(!small)) ps[!small] <- Recall(z[!small])
		return(ps)
	}
	x <- 1/z^2
	tail <- ((x * (-1/12 + ((x * (1/120 + ((x * (-1/252 + ((
		x * (1/240 + ((x * (-1/132 + ((x * (691/32760 + (
		(x * (-1/12 + (3617 * x)/8160)))))))))))))))))))))
	1/(2 * z) - tail
}
Digamma <- function(link = "log") {
#	Digamma generalized linear model family
#	Gordon Smyth, smyth@wehi.edu.au
#  3 July 1998.  Last revised 9 Dec 2002.
#
#	improve on the link deparsing code in quasi()
	linkarg <- substitute(link)
	if (is.expression(linkarg) || is.call(linkarg)) {
		linkname <- deparse(linkarg)
	} else if(is.character(linkarg)) { 
		linkname <- linkarg
		link <- make.link(linkarg)
	} else if(is.numeric(linkarg)) {
		linkname <- paste("power(",linkarg,")",sep="")
		link <- make.link(linkarg)
	} else {
		linkname <- deparse(linkarg) 
		link <- make.link(linkname)
	}
	validmu <- function(mu) all(mu>0)
	dev.resids <- function(y, mu, wt) wt * unitdeviance.digamma(y,mu)
	initialize <- expression({
		if (any(y <= 0)) stop(paste("Non-positive values not", "allowed for the Digamma family"))
		n <- rep(1, nobs)
		mustart <- y
	})
	aic <- function(y, n, mu, wt, dev) NA
	structure(list(
		family = "Digamma",
		variance = varfun.digamma, dev.resids = dev.resids, aic = aic,
		link = linkname,
		linkfun = link$linkfun, linkinv = link$linkinv, mu.eta = link$mu.eta,
		valideta = link$valideta, validmu = validmu, initialize = initialize, 
		class = "family"))
}

cumulant.digamma <- function(theta)
#	Cumulant function for the Digamma family
#	GKS  3 July 98
	2*( theta*(log(-theta)-1) + lgamma(-theta) )

meanval.digamma <- function(theta)
#	Mean value function for the Digamma family
#	GKS  3 July 98
	2*( log(-theta) - digamma(-theta) )

d2cumulant.digamma <- function(theta)
#	2nd derivative of cumulant function for Digamma family
#	GKS  3 July 98
	2*( 1/theta + trigamma(-theta) )

canonic.digamma <- function(mu) {
#	Canonical mapping for Digamma family
#	Solve meanval.digamma(theta) = mu for theta
#	GKS  3 July 98
#
#	Starting value from -log(-theta) =~ log(mu)
	mlmt <- log(mu)	
	theta <- -exp(-mlmt)

	for (i in 1:3) {
		mu1 <- meanval.digamma(theta)
		v <- d2cumulant.digamma(theta)
		deriv <- -v/mu1*theta
		mlmt <- mlmt - log(mu1/mu)/deriv	
		theta <- -exp(-mlmt)
	}
	theta
}

varfun.digamma <- function(mu) {
#	Variance function for Digamma family
#	GKS  3 July 98
#
	theta <- canonic.digamma(mu)
	2*( 1/theta + trigamma(-theta) )
}

unitdeviance.digamma <- function(y,mu) {
#	Unit deviance for Digamma family
#	GKS  3 July 98
#
	thetay <- canonic.digamma(y)
	theta <- canonic.digamma(mu)
	2*( y*(thetay-theta) - (cumulant.digamma(thetay)-cumulant.digamma(theta)) )
}
#  NUMERICAL INTEGRATION

gauss.quad <- function(n,kind="legendre",alpha=0,beta=0) {
#	Calculate nodes and weights for Guassian quadrature.
#	Adapted from Netlib routine gaussq.f
#	Gordon Smyth, Walter and Eliza Hall Institute, smyth@wehi.edu.au
#	4 Sept 2002

	kind <- match.arg(kind,c("legendre","chebyshev1","chebyshev2","hermite","jacobi","laguerre"))
	i <- 1:n
	i1 <- 1:(n-1)
	switch(kind, legendre={
		muzero <- 2
		a <- rep(0,n)
		b <- i1/sqrt(4*i1^2-1)
	}, chebyshev1={
		muzero <- pi
		a <- rep(0,n)
		b <- rep(0.5,n-1)
		b[1] <- sqrt(0.5)
	}, chebyshev2={
		muzero <- pi/2
		a <- rep(0,n)
		b <- rep(0.5,n-1)
	}, hermite={
		muzero <- sqrt(pi)
		a <- rep(0,n)
		b <- sqrt(i1/2)
	}, jacobi={
		ab <- alpha+beta
		muzero <- 2^(ab+1)*gamma(alpha+1)*gamma(beta+1)/gamma(ab+2)
		a <- i
		a[1] <- (beta-alpha)/(ab+2)
		i2 <- 2:n
		abi <- ab+2*i2
		a[i2] <- (beta^2-alpha^2)/(abi-2)/abi
		b <- i1
		b[1] <- sqrt(4*(alpha+1)*(beta+1)/(ab+2)^2/(ab+3))
		i2 <- 2:(n-1)
		abi <- ab+2*i2
		b[i2] <- sqrt(4*i2*(i2+alpha)*(i2+beta)*(i2+ab)/(abi^2-1)/abi^2)
	}, laguerre={
		a <- 2*i-1+alpha
		b <- sqrt(i1*(i1+alpha))
		muzero <- gamma(alpha+1)
	})
	A <- rep(0,n*n)
	A[(n+1)*(i-1)+1] <- a
	A[(n+1)*(i1-1)+2] <- b
	A[(n+1)*i1] <- b
	dim(A) <- c(n,n)
	vd <- eigen(A,symmetric=TRUE)
	w <- rev(as.vector( vd$vectors[1,] ))
	w <- muzero * w^2
	x <- rev( vd$values )
	list(nodes=x,weights=w)
}

gauss.quad.prob <- function(n,dist="uniform",l=0,u=1,mu=0,sigma=1,alpha=1,beta=1) {
#	Calculate nodes and weights for Guassian quadrature using probability densities.
#	Adapted from Netlib routine gaussq.f
#	Gordon Smyth, Walter and Eliza Hall Institute, smyth@wehi.edu.au
#	4 Sept 2002

	dist <- match.arg(dist,c("uniform","beta1","beta2","normal","beta","gamma"))
	if(dist=="beta" && alpha==0.5 && beta==0.5) dist <- "beta1"
	if(dist=="beta" && alpha==1.5 && beta==1.5) dist <- "beta2"
	i <- 1:n
	i1 <- 1:(n-1)
	switch(dist, uniform={
		a <- rep(0,n)
		b <- i1/sqrt(4*i1^2-1)
	}, beta1={
		a <- rep(0,n)
		b <- rep(0.5,n-1)
		b[1] <- sqrt(0.5)
	}, beta2={
		a <- rep(0,n)
		b <- rep(0.5,n-1)
	}, normal={
		a <- rep(0,n)
		b <- sqrt(i1/2)
	}, beta={
		ab <- alpha+beta
		a <- i
		a[1] <- (alpha-beta)/ab
		i2 <- 2:n
		abi <- ab-2+2*i2
		a[i2] <- ((alpha-1)^2-(beta-1)^2)/(abi-2)/abi
		b <- i1
		b[1] <- sqrt(4*alpha*beta/ab^2/(ab+1))
		i2 <- 2:(n-1)
		abi <- ab-2+2*i2
		b[i2] <- sqrt(4*i2*(i2+alpha-1)*(i2+beta-1)*(i2+ab-2)/(abi^2-1)/abi^2)
	}, gamma={
		a <- 2*i+alpha-2
		b <- sqrt(i1*(i1+alpha-1))
	})
	A <- rep(0,n*n)
	A[(n+1)*(i-1)+1] <- a
	A[(n+1)*(i1-1)+2] <- b
	A[(n+1)*i1] <- b
	dim(A) <- c(n,n)
	vd <- eigen(A,symmetric=TRUE)
	w <- rev(as.vector( vd$vectors[1,] ))^2
	x <- rev( vd$values )
	switch(dist,
		uniform = x <- l+(u-l)*(x+1)/2,
		beta1 = x <- (x+1)/2,
		beta2 = x <- (x+1)/2,
		normal = x <- mu + sqrt(2)*sigma*x,
		beta = x <- (x+1)/2,
		gamma = x <- beta*x)
	list(nodes=x,weights=w)
}
#  GLMGAM.R

glmgam.fit <- function(X,y,start=NULL,tol=1e-6,maxit=50,trace=FALSE) {
#  Fit gamma generalized linear model with identity link
#  by Levenberg damped Fisher scoring
#  Gordon Smyth
#  12 Mar 2003.  Last revised 20 March 2004.

#  check input
X <- as.matrix(X)
n <- nrow(X)
p <- ncol(X)
if(p > n) stop("More columns than rows in X")
y <- as.vector(y)
if(n != length(y)) stop("length(y) not equal to nrow(X)")
if(n == 0) return(list(coefficients=numeric(0),fitted.values=numeric(0),deviance=numeric(0)))
if(!(all(is.finite(y)) || all(is.finite(X)))) stop("All values must be finite and non-missing")
if(any(y < 0)) stop("y must be non-negative")
maxy <- max(y)
if(maxy==0) return(list(coefficients=rep(0,p),fitted.values=rep(0,n),deviance=NA))
y1 <- pmax(y,maxy*1e-3)

#  starting values
if(is.null(start)) {
	fit <- lm.fit(X,y)
	beta <- fit$coefficients
	mu <- fit$fitted.values
	if(any(mu < 0)) {
		fit <- lm.wfit(X,y,1/y1^2)
		beta <- fit$coefficients
		mu <- fit$fitted.values
	}
	if(any(mu < 0)) {
		fit <- lm.fit(X,rep(mean(y),n))
		beta <- fit$coefficients
		mu <- fit$fitted.values
	}
	if(any(mu < 0)) {
		samesign <- apply(X>0,2,all) | apply(X<0,2,all)
		if(any(samesign)) {
			i <- (1:p)[samesign][1]
			beta <- rep(0,p)
			beta[i] <- lm.wfit(X[,i,drop=FALSE],y,1/y1^2)$coef
			mu <- X[,i] * beta[i]
		} else
			return(list(coefficients=rep(0,p),fitted.values=rep(0,n),deviance=Inf))
	}
} else {
	beta <- start
	mu <- X %*% beta
}
if(any(mu<0)) stop("Starting values give negative fitted values")

deviance.gamma <- function(y,mu) {
	if(any(mu<0)) return(Inf)
	o <- (y < 1e-15) & (mu < 1e-15)
	if(any(o)) {
		if(all(o)) {
			dev <- 0
		} else {
			y1 <- y[!o]
			mu1 <- y[!o]
			dev <- 2*sum( (y1-mu1)/mu1 - log(y1/mu1) )
		}
	} else {
		dev <- 2*sum( (y-mu)/mu - log(y/mu) )
	}
}

dev <- deviance.gamma(y,mu)

# Scoring iteration with Levenberg damping
iter <- 0
if(trace) cat("Iter =",iter,", Dev =",dev," Beta",beta,"\n")
repeat {
	iter <- iter+1

	# information matrix
	v <- mu^2
	v <- pmax(v,max(v)/10^3)
	XVX <- crossprod(X,vecmat(1/v,X))
	maxinfo <- max(diag(XVX))
	if(iter==1) {
		lambda <- abs(mean(diag(XVX)))/p
		I <- diag(p)
	}

	# score vector
	dl <- crossprod(X,(y-mu)/v)

	# Levenberg damping
	betaold <- beta
	devold <- dev
	lev <- 0
	repeat {
		lev <- lev+1

		# trial step
		R <- chol(XVX + lambda*I)
		dbeta <- backsolve(R,backsolve(R,dl,transpose=TRUE))
		beta <- betaold + dbeta
		mu <- X %*% beta
		dev <- deviance.gamma(y,mu)
		if(dev <= devold || dev/max(mu) < 1e-15) break

		# exit if too much damping
		if(lambda/maxinfo > 1e15) {
			beta <- betaold
			warning("Too much damping - convergence tolerance not achievable")
			break
		}

		# step not successful so increase damping
		lambda <- 2*lambda
		if(trace) cat("Damping increased to",lambda,"\n")
	}

	# iteration output
	if(trace) cat("Iter =",iter,", Dev =",dev," Beta",beta,"\n")

	# keep exiting if too much damping
	if(lambda/maxinfo > 1e15) break

	# decrease damping if successful at first try
	if(lev==1) lambda <- lambda/10

	# test for convergence
	if( crossprod(dl,dbeta) < tol || dev/max(mu) < 1e-15) break

	# test for iteration limit
	if(iter > maxit) {
		iter <- maxit+1
		break
	}
}

beta <- drop(beta)
names(beta) <- colnames(X)
list(coefficients=beta,fitted.values=as.vector(mu),deviance=dev,maxit=maxit,iter=iter)
}
meanT <- function(y1,y2) {
#  Mean t-statistic difference between two groups of growth curves
#  Columns are time points, rows are individuals
#  Gordon Smyth
#  14 Feb 2003

	if(is.null(dim(y1)) || is.null(dim(y2))) return(NA)
	y1 <- as.matrix(y1)
	y2 <- as.matrix(y2)
	if(ncol(y1) != ncol(y2)) stop("Number of time points must match")
	m1 <- colMeans(y1,na.rm=TRUE)
	m2 <- colMeans(y2,na.rm=TRUE)
	v1 <- apply(y1,2,var,na.rm=TRUE)
	v2 <- apply(y2,2,var,na.rm=TRUE)
	n1 <- apply(!is.na(y1),2,sum)
	n2 <- apply(!is.na(y2),2,sum)
	s <- ( (n1-1)*v1 + (n2-1)*v2 ) / (n1+n2-2)
	t.stat <- (m1-m2) / sqrt(s*(1/n1+1/n2))
	weighted.mean(t.stat,w=(n1+n2-2)/(n1+n2),na.rm=TRUE)
}

compareTwoGrowthCurves <- function(group,y,nsim=100,fun=meanT) {
#  Permutation test between two groups of growth curves
#  Columns are time points, rows are individuals
#  Gordon Smyth
#  14 Feb 2003

	group <- as.vector(group)
	g <- unique(group)
	if(length(g) != 2) stop("Must be exactly 2 groups")
	stat.obs <- fun(y[group==g[1],,drop=FALSE], y[group==g[2],,drop=FALSE])
	asbig <- 0
	for (i in 1:nsim) {
		pgroup <- sample(group)
		stat <- fun(y[pgroup==g[1],,drop=FALSE], y[pgroup==g[2],,drop=FALSE])
		if(abs(stat) >= abs(stat.obs)) asbig <- asbig+1
	}
	list(stat=stat.obs, p.value=asbig/nsim) 
}

compareGrowthCurves <- function(group,y,levels=NULL,nsim=100,fun=meanT,times=NULL,verbose=TRUE,adjust="holm") {
#  All pairwise permutation tests between groups of growth curves
#  Columns of y are time points, rows are individuals
#  Gordon Smyth
#  14 Feb 2003.  Last modified 17 Nov 2003.

	group <- as.character(group)
	if(is.null(levels)) {
		tab <- table(group)
		tab <- tab[tab >= 2]
		lev <- names(tab)
	} else
		lev <- as.character(levels)
	nlev <- length(lev)
	if(nlev < 2) stop("Less than 2 groups to compare")
	if(is.null(dim(y))) stop("y must be matrix-like")
	y <- as.matrix(y)
	if(!is.null(times)) y <- y[,times,drop=FALSE]

	g1 <- g2 <- rep("",nlev*(nlev-1)/2)
	stat <- pvalue <- rep(0,nlev*(nlev-1)/2)
	pair <- 0
	for (i in 1:(nlev-1)) {
		for (j in (i+1):nlev) {
			if(verbose) cat(lev[i],lev[j])
			pair <- pair+1
			sel <- group %in% c(lev[i],lev[j])
			out <- compareTwoGrowthCurves(group[sel],y[sel,,drop=FALSE],nsim=nsim,fun=fun)
			if(verbose) cat("\ ",round(out$stat,2),"\n")
			g1[pair] <- lev[i]
			g2[pair] <- lev[j]
			stat[pair] <- out$stat
			pvalue[pair] <- out$p.value
		}
	}
	tab <- data.frame(Group1=g1,Group2=g2,Stat=stat,P.Value=pvalue)
	tab$adj.P.Value <- p.adjust(pvalue,method=adjust)
	tab
}
hommel.test <-
#	Multiple testing from Hommel (1988).
#	Similar but very slightly more powerful that Hochberg (1988).
#	Controls Family-Wise Error rate for hypotheses which are independent or
#	which satisfy the free-association condition of Simes (1986).
#	Gordon Smyth, Walter and Eliza Hall Institute, smyth@wehi.edu.au
#	29 Aug 2002

function(p,alpha=0.05) {
	n <- length(p)
	i <- 1:n
	po <- sort(p)
	j <- n
	repeat {
		k <- 1:j
		if(all( po[n - j + k] > k * alpha / j )) break
		j <- j-1
		if(j == 0) break
	}
	p >= alpha/j
}
dinvgauss <- function(x, mu, lambda = 1)
{
#  Density of inverse Gaussian distribution
#  GKS  15 Jan 98
#
	if(any(mu<=0)) stop("mu must be positive")
	if(any(lambda<=0)) stop("lambda must be positive")
	d <- ifelse(x>0,sqrt(lambda/(2*pi*x^3))*exp(-lambda*(x-mu)^2/(2*mu^2*x)),0)
	if(!is.null(Names <- names(x)))
		names(d) <- rep(Names, length = length(d))
	d
}

pinvgauss <- function(q, mu, lambda = 1)
{
#  Inverse Gaussian distribution function
#  GKS  15 Jan 98
#
	if(any(mu<=0)) stop("mu must be positive")
	if(any(lambda<=0)) stop("lambda must be positive")
	n <- length(q)
	if(length(mu)>1 && length(mu)!=n) mu <- rep(mu,length=n)
	if(length(lambda)>1 && length(lambda)!=n) lambda <- rep(lambda,length=n)
	lq <- sqrt(lambda/q)
	qm <- q/mu
	p <- ifelse(q>0,pnorm(lq*(qm-1))+exp(2*lambda/mu)*pnorm(-lq*(qm+1)),0)
	if(!is.null(Names <- names(q)))
		names(p) <- rep(Names, length = length(p))
	p
}

rinvgauss <- function(n, mu, lambda = 1)
{
#  Random variates from inverse Gaussian distribution
#  Reference:
#      Chhikara and Folks, The Inverse Gaussian Distribution,
#      Marcel Dekker, 1989, page 53.
#  GKS  15 Jan 98
#
	if(any(mu<=0)) stop("mu must be positive")
	if(any(lambda<=0)) stop("lambda must be positive")
	if(length(n)>1) n <- length(n)
	if(length(mu)>1 && length(mu)!=n) mu <- rep(mu,length=n)
	if(length(lambda)>1 && length(lambda)!=n) lambda <- rep(lambda,length=n)
	y2 <- rchisq(n,1)
	u <- runif(n)
	r1 <- mu/(2*lambda) * (2*lambda + mu*y2 - sqrt(4*lambda*mu*y2 + mu^2*y2^2))
	r2 <- mu^2/r1
	ifelse(u < mu/(mu+r1), r1, r2)
}

qinvgauss  <- function(p, mu, lambda = 1)
{
#  Quantiles of the inverse Gaussian distribution
#  Dr Paul Bagshaw
#  Centre National d'Etudes des Telecommunications (DIH/DIPS)
#  Technopole Anticipa, France
#  paul.bagshaw@cnet.francetelecom.fr
#  23 Dec 98
#
  if(any(mu <= 0.))
    stop("mu must be positive")
  if(any(lambda <= 0.))
    stop("lambda must be positive")
  n <- length(p)
  if(length(mu) > 1 && length(mu) != n)
    mu <- rep(mu, length = n)
  if(length(lambda) > 1 && length(lambda) != n)
    lambda <- rep(lambda, length = n)
  thi <- lambda / mu
  U <- qnorm (p)
  r1 <- 1 + U / sqrt (thi) + U^2 / (2 * thi) + U^3 / (8 * thi * sqrt(thi))
  x <- r1
  for (i in 1:10) {
    cum <- pinvgauss (x, 1., thi)
    dx <- (cum - p) / dinvgauss (x, 1., thi)
    dx <- ifelse (is.finite(dx), dx, ifelse (p > cum, -1, 1))
    dx[dx < -1] <- -1
    if (all(dx == 0.)) break
    x <- x - dx
  }
  x * mu
}
matvec <- function(M,v) {
#	Multiply the columns of matrix by the elements of a vector,
#	i.e., compute M %*% diag(v)
#	Gordon Smyth
#	5 July 1999
#
	v <- as.vector(v)
	M <- as.matrix(M)
	if(length(v)!=dim(M)[2]) stop("matvec: Dimensions do not match")
	t(v * t(M))
}

vecmat <- function(v,M) {
#	Multiply the rows of matrix by the elements of a vector,
#	i.e., compute diag(v) %*% M
#	Gordon Smyth
#	5 July 1999
#
	v <- as.vector(v)
	M <- as.matrix(M)
	if(length(v)!=dim(M)[1]) stop("vecmat: Dimensions do not match")
	v * M
}
power.fisher.test <- function(p1,p2,n1,n2,alpha=0.05,nsim=100) {
#	Calculation of power for Fisher's exact test for
#	comparing two proportions
#	Gordon smyth
#	3 June 2003

	y1 <- rbinom(nsim,size=n1,prob=p1)
	y2 <- rbinom(nsim,size=n2,prob=p2)
	y <- cbind(y1,n1-y1,y2,n2-y2)
	p.value <- rep(0,nsim)
	for (i in 1:nsim)
		p.value[i] <- fisher.test(matrix(y[i,],2,2))$p.value
	mean(p.value < alpha)
}
## QRES.R

qresiduals <- qresid <- function(glm.obj, dispersion=NULL)
#	Wrapper function for quantile residuals
#	Peter K Dunn
#	28 Sep 2004.  Last modified 5 Oct 2004.
{
glm.family <- glm.obj$family$family
if(substr(glm.family,1,17)=="Negative Binomial") glm.family <- "nbinom"
switch(glm.family,
   binomial = qres.binom( glm.obj),
   poisson = qres.pois(glm.obj),
   Gamma = qres.gamma(glm.obj, dispersion),
   inverse.gaussian = qres.invgauss(glm.obj, dispersion),
   Tweedie = qres.tweedie(glm.obj, dispersion),
   nbinom = qres.nbinom(glm.obj),
   qres.default(glm.obj, dispersion)
)}

qres.binom <- function(glm.obj)
#	Randomized quantile residuals for binomial glm
#	Gordon Smyth
#	20 Oct 96.  Last modified 25 Jan 02.
{
	p <- fitted(glm.obj)
	y <- glm.obj$y
	if(!is.null(glm.obj$prior.weights))
		n <- glm.obj$prior.weights
	else
		n <- rep(1,length(y))
	y <- n * y
	a <- pbinom(y - 1, n, p)
	b <- pbinom(y, n, p)
	u <- runif(n = length(y), min = a, max = b)
	qnorm(u)
}

qres.pois <- function(glm.obj)
#	Quantile residuals for Poisson glm
#	Gordon Smyth
#	28 Dec 96
{
	y <- glm.obj$y
	mu <- fitted(glm.obj)
	a <- ppois(y - 1, mu)
	b <- ppois(y, mu)
	u <- runif(n = length(y), min = a, max = b)
	qnorm(u)
}

qres.gamma <- function(glm.obj, dispersion = NULL)
#	Quantile residuals for gamma glm
#	Gordon Smyth
#	28 Dec 96.  Last modified 10 Jan 97
{
	mu <- fitted(glm.obj)
	y <- glm.obj$y
	df <- glm.obj$df.residual
	w <- glm.obj$prior.weights
	if(is.null(w))
		w <- 1
	if(is.null(dispersion))
		dispersion <- sum(w * ((y - mu)/mu)^2)/df
	u <- pgamma((w * y)/mu/dispersion, w/dispersion)
	qnorm(u)
}

qres.invgauss <- function(glm.obj, dispersion = NULL)
#	Quantile residuals for inverse Gaussian glm
#	Gordon Smyth
#	15 Jan 98
{
	mu <- fitted(glm.obj)
	y <- glm.obj$y
	df <- glm.obj$df.residual
	w <- glm.obj$prior.weights
	if(is.null(w))
		w <- 1
	if(is.null(dispersion))
		dispersion <- sum(w * (y - mu)^2 / (mu^2*y)) / df
	u <- pinvgauss(y, mu, lambda=1/dispersion)
	qnorm(u)
}

qres.nbinom <- function(glm.obj)
{
#	Quantile residuals for Negative Binomial glm
#	Gordon Smyth
#	22 Jun 97.  Last modified 5 Oct 2004.
#
	y <- glm.obj$y
	size <- glm.obj$call$family$theta
	mu <- fitted(glm.obj)
	p <- size/(mu + size)
	a <- ifelse(y > 0, pbeta(p, size, pmax(y, 1)), 0)
	b <- pbeta(p, size, y + 1)
	u <- runif(n = length(y), min = a, max = b)
	qnorm(u)
}

qres.tweedie <- function(glm.obj, dispersion = NULL)
#	Quantile residuals for Tweedie glms
#	Gordon Smyth
#	29 April 98.  Last modified 5 Oct 2004.
{
	require("tweedie")
	mu <- fitted(glm.obj)
	y <- glm.obj$y
	df <- glm.obj$df.residual
	w <- glm.obj$prior.weights
	if(is.null(w))
		w <- 1
	p <- get("p",envir=environment(glm.obj$family$variance))
	if(is.null(dispersion))
		dispersion <- sum((w * (y - mu)^2)/mu^p)/df
	u <- ptweedie(q=y, power=p, mu=fitted(glm.obj), phi=dispersion/w)
	if(p>1&&p<2)
		u[y == 0] <- runif(sum(y == 0), min = 0, max = u[y == 0])
	qnorm(u)
}

qres.default <- function(glm.obj, dispersion=NULL)
#	Quantile residuals for Gaussian and default glms
#	Gordon Smyth
#	5 Oct 2004.
{
	r <- residuals(glm.obj, type="deviance")
	if(is.null(dispersion)) {
		df.r <- glm.obj$df.residual
		if(df.r > 0) {
			if(any(glm.obj$weights==0)) warning("observations with zero weight ", "not used for calculating dispersion")
	        dispersion <- sum(glm.obj$weights*glm.obj$residuals^2)/df.r
	    } else
	    	dispersion <- 1
    }
    r/sqrt(dispersion)
}

#  RANBLOCK.R

randomizedBlock <- function(formula, random, weights=NULL, only.varcomp=FALSE, data=list(), subset=NULL, contrasts=NULL, tol=1e-6, maxit=50, trace=FALSE)
#	REML for mixed linear models
#	Gordon Smyth, Walter and Eliza Hall Institute
#	28 Jan 2003.  Last revised 20 March 2004.
{
#	Extract model from formula
	cl <- match.call()
	mf <- match.call(expand.dots = FALSE)
	mf$only.varcomp <- mf$tol <- mf$tol <- mf$maxit <- NULL
	mf$drop.unused.levels <- TRUE
	mf[[1]] <- as.name("model.frame")
	mf <- eval(mf, parent.frame())
	mt <- attr(mf, "terms")
	xvars <- as.character(attr(mt, "variables"))[-1]
	if((yvar <- attr(mt,"response")) > 0) xvars <- xvars[-yvar]
	xlev <- if(length(xvars) > 0) {
		xlev <- lapply(mf[xvars], levels)
		xlev[!sapply(xlev, is.null)]
	}
	y <- model.response(mf, "numeric")
	w <- model.weights(mf)
	x <- model.matrix(mt, mf, contrasts)
	random <- mf[["(random)"]]

#	Missing values not allowed
	if(any(is.na(y)) || any(is.na(x)) || any(is.na(random))) stop("Missing values not allowed")
	if(!is.null(weights)) if(any(is.na(weights))) stop("Missing values not allowed")

#	Design matrix for random effects
	lev <- unique.default(random)
	z <- 0 + (matrix(random,length(random),length(lev)) == t(matrix(lev,length(lev),length(random))))

	randomizedBlockFit(y,x,z,w=w,only.varcomp=only.varcomp,tol=tol,maxit=maxit,trace=trace)
}

randomizedBlockFit <- function(y,X,Z,w=NULL,only.varcomp=FALSE,tol=1e-6,maxit=50,trace=FALSE) {
#	Restricted maximum likelihood estimation for mixed linear models.
#	Fits the model  Y = X*BETA + Z*U + E  where BETA is fixed
#	and U is random.
#
#	GAMMA holds the variance components.  The errors E and
#	random effects U are assumed to have covariance matrices
#	EYE*GAMMA(1) and EYE*GAMMA(2) respectively.

#	Gordon Smyth, Walter and Eliza Hall Institute
#	Matlab version 19 Feb 94.  Converted to R, 28 Jan 2003.
#	Last revised 20 Mar 2004.

#  Prior weights
if(!is.null(w)) {
	sw <- sqrt(w)
	y <- sw * y
	X <- sw * X
}

#  Find null space Q of X
X <- as.matrix(X)
Z <- as.matrix(Z)
mx <- nrow(X)
nx <- ncol(X)
s <- La.svd(X,nu=mx,nv=0)
if(s$d[1] < 1e-15)
	zero1 <- 1   # X is entirely zero
else {
	zeroeig <- abs(s$d/s$d[1]) < 1e-15
	if(any(zeroeig))
		zero1 <- min((1:nx)[zeroeig])
	else
		zero1 <- nx+1
}
#  Are there any df to estimate error?
if(zero1 > mx) return(list(varcomp=rep(NA,2)))
Q <- s$u[,zero1:mx,drop=FALSE]

#  Apply Q to Z and transform to independent observations
mq <- ncol(Q)
s <- La.svd(crossprod(Q,Z),nu=mq,nv=0)
uqy <- crossprod(s$u,(crossprod(Q,y)))
d <- rep(0,mq)
d[1:length(s$d)] <- s$d^2
dx <- cbind(Residual=1,Block=d)
dy <- uqy^2

#  Try unweighted starting values
dfit <- lm.fit(dx,dy)
varcomp <- dfit$coefficients
dfitted.values <- dfit$fitted.values

#  Main fit
if(mq > 2 && sum(abs(d)>1e-15)>1 && var(d)>1e-15) {
	if(all(dfitted.values >= 0))
		start <- dfit$coefficients
	else
		start <- c(Residual=mean(dy),Block=0)
#	fit gamma glm identity link to dy with dx as covariates
	dfit <- glmgam.fit(dx,dy,start=start,tol=tol,maxit=maxit,trace=trace)
	varcomp <- dfit$coefficients
	dfitted.values <- dfit$fitted.values
}
out <- list(varcomp=dfit$coef)
if(only.varcomp) return(out)

#  Standard errors for variance components
dinfo <- crossprod(dx,vecmat(1/dfitted.values^2,dx))
out$se.varcomp=sqrt(2*diag(chol2inv(chol(dinfo))))

#  fixed effect estimates
s <- La.svd(Z,nu=mx,nv=0)
d <- rep(0,mx)
d[1:length(s$d)] <- s$d^2
v <- drop( cbind(Residual=1,Block=d) %*% varcomp )
mfit <- lm.wfit(crossprod(s$u,X),crossprod(s$u,y),1/v)
out <- c(out,mfit)
out$se.coefficients <- sqrt(diag(chol2inv(mfit$qr$qr)))

out
}

randomizedBlockFit <- function(y,X,Z,w=NULL,only.varcomp=FALSE,tol=1e-6,maxit=50,trace=FALSE) {
#	Restricted maximum likelihood estimation for mixed linear models.
#	Fits the model  Y = X*BETA + Z*U + E  where BETA is fixed
#	and U is random.
#
#	GAMMA holds the variance components.  The errors E and
#	random effects U are assumed to have covariance matrices
#	EYE*GAMMA(1) and EYE*GAMMA(2) respectively.

#	Gordon Smyth, Walter and Eliza Hall Institute
#	Matlab version 19 Feb 94.  Converted to R, 28 Jan 2003.
#	Last revised 24 Mar 2004.

#  Prior weights
if(!is.null(w)) {
	sw <- sqrt(w)
	y <- sw * y
	X <- sw * X
}

#  Find null space Q of X
X <- as.matrix(X)
Z <- as.matrix(Z)
mx <- nrow(X)
nx <- ncol(X)
nz <- ncol(Z)
fit <- lm.fit(X,cbind(Z,y))
r <- fit$rank
QtZ <- fit$effects[(r+1):mx,1:nz]

#  Apply Q to Z and transform to independent observations
mq <- mx-r
if(mq == 0) return(list(varcomp=c(NA,NA)))
s <- La.svd(QtZ,nu=mq,nv=0)
uqy <- crossprod(s$u,fit$effects[(r+1):mx,nz+1])
d <- rep(0,mq)
d[1:length(s$d)] <- s$d^2
dx <- cbind(Residual=1,Block=d)
dy <- uqy^2

#  Try unweighted starting values
dfit <- lm.fit(dx,dy)
varcomp <- dfit$coefficients
dfitted.values <- dfit$fitted.values

#  Main fit
if(mq > 2 && sum(abs(d)>1e-15)>1 && var(d)>1e-15) {
	if(all(dfitted.values >= 0))
		start <- dfit$coefficients
	else
		start <- c(Residual=mean(dy),Block=0)
#	fit gamma glm identity link to dy with dx as covariates
	dfit <- glmgam.fit(dx,dy,start=start,tol=tol,maxit=maxit,trace=trace)
	varcomp <- dfit$coefficients
	dfitted.values <- dfit$fitted.values
}
out <- list(varcomp=dfit$coef)
if(only.varcomp) return(out)

#  Standard errors for variance components
dinfo <- crossprod(dx,vecmat(1/dfitted.values^2,dx))
out$se.varcomp=sqrt(2*diag(chol2inv(chol(dinfo))))

#  fixed effect estimates
s <- La.svd(Z,nu=mx,nv=0)
d <- rep(0,mx)
d[1:length(s$d)] <- s$d^2
v <- drop( cbind(Residual=1,Block=d) %*% varcomp )
mfit <- lm.wfit(crossprod(s$u,X),crossprod(s$u,y),1/v)
out <- c(out,mfit)
out$se.coefficients <- sqrt(diag(chol2inv(mfit$qr$qr)))

out
}

remlscore <- function(y,X,Z,trace=FALSE,tol=1e-5,maxit=40) {
#
#  Mean-variance fit by REML scoring
#  Fit normal(mu,phi) model to y with
#  mu=X%*%beta and log(phi)=Z%*%gam
#
#  Gordon Smyth, Walter and Eliza Hall Institute, smyth@wehi.edu.au
#  11 Sept 2000.  Last modified 10 Dec 2002.

n <- length(y)
p <- dim(X)[2]
q <- dim(Z)[2]
const <- n*log(2*pi)

# initial residuals from unweighted regression
fitm <- lm.fit(X,y)
if(fitm$qr$rank < p) stop("X is of not of full column rank")
Q <- qr.Q(fitm$qr)
h <- as.vector(Q^2 %*% array(1, c(p, 1)))
d <- fitm$residuals^2

# starting values
# use of weights guarantee that regression can be computed even if 1-h = 0
wd <- 1-h
zd <- log( d/(1-h) )+1.27
fitd <- lm.wfit(Z,zd,wd)
gam <- ifelse(is.na(fitd$coef),0,fitd$coef)
g <- fitd$fitted.values
phi <- exp(g)
wm <- 1/phi
fitm <- lm.wfit(X,y,wm)
d <- fitm$residuals^2
dev <- sum(d/phi)+sum(log(phi))+const+2*log(prod(abs(diag(fitm$qr$qr))))

# reml scoring
iter <- 0
if(trace) cat("Iter =",iter,", Dev =",dev," Gamma",gam,"\n")
Q2 <- array(0,c(n,p*(p+1)/2))
repeat {
	iter <- iter+1

	# information matrix and leverages
	Q <- qr.qy(fitm$qr, diag(1, nrow = n, ncol = p))
	j0 <- 0
	for(k in 0:(p-1)) {
		Q2[ ,(j0+1):(j0+p-k)] <- Q[ ,1:(p-k)] * Q[ ,(k+1):p]
		j0 <- j0+p-k
	}
	Q2[ ,(p+1):(p*(p+1)/2)] <- sqrt(2) * Q2[ ,(p+1):(p*(p+1)/2)]
	h <- drop( Q2[ ,1:p] %*% array(1,c(p,1)) )
	Q2Z <- t(Q2) %*% Z
	ZVZ <- ( t(Z) %*% vecmat(1-2*h,Z) + t(Q2Z) %*% Q2Z	)/2
	maxinfo <- max(diag(ZVZ))
	if(iter==1) {
		lambda <- abs(mean(diag(ZVZ)))/q
		I <- diag(q)
	}

	# score vector
	zd <- ( d - (1-h)*phi ) / phi
	dl <- crossprod(Z,zd)/2

	# Levenberg damping
	gamold <- gam
	devold <- dev
	lev <- 0
	repeat {
		lev <- lev+1

		# trial step
		R <- chol(ZVZ + lambda*I)
		dgam <- backsolve(R,backsolve(R,dl,transpose=TRUE))
		gam <- gamold + dgam
		phi <- as.vector(exp( Z %*% gam ))
		wm <- 1/phi
		fitm <- lm.wfit(X,y,wm)
		d <- fitm$residuals^2
		dev <- sum(d/phi)+sum(log(phi))+const+2*log(prod(abs(diag(fitm$qr$qr))))
		if(dev < devold - 1e-15) break

		# exit if too much damping
		if(lambda/maxinfo > 1e15) {
			gam <- gamold
			warning("Too much damping - convergence tolerance not achievable")
			break
		}

		# step not successful so increase damping
		lambda <- 2*lambda
		if(trace) cat("Damping increased to",lambda,"\n")
	}

	# iteration output
	if(trace) cat("Iter =",iter,", Dev =",dev," Gamma",gam,"\n")

	# keep exiting if too much damping
	if(lambda/maxinfo > 1e15) break

	# decrease damping if successful at first try
	if(lev==1) lambda <- lambda/10

	# test for convergence
	if( crossprod(dl,dgam) < tol ) break

	# test for iteration limit
	if(iter > maxit) {
		warning("reml: Max iterations exceeded")
		break
	}
}

# Nominal standard errors
se.gam <- sqrt(diag(chol2inv(chol(ZVZ))))
se.beta <- sqrt(diag(chol2inv(qr.R(fitm$qr))))

list(beta=fitm$coef,se.beta=se.beta,gamma=gam,se.gam=se.gam,mu=fitm$fitted,phi=phi,deviance=dev,h=h)
}
remlscoregamma <- function(y,X,Z,mlink="log",dlink="log",trace=FALSE,tol=1e-5,maxit=40) {
#
#  Mean-dispersion fit by REML scoring for gamma responses
#  Fit ED(mu,phi) model to y with
#  g(mu)=X%*%beta and f(phi)=Z%*%gam
#
#  Gordon Smyth, Walter and Eliza Hall Institute
#  16 Dec 2002.

n <- length(y)
X <- as.matrix(X)
if(is.null(colnames(X))) colnames(X) <- paste("X",as.character(1:ncol(X)),sep="")
Z <- as.matrix(Z)
if(is.null(colnames(Z))) colnames(Z) <- paste("Z",as.character(1:ncol(Z)),sep="")
q <- dim(Z)[2]
const <- 2*sum(log(y))

# Link functions
mli <- make.link(mlink)
dli <- make.link(dlink)

# Mean family
f <- Gamma()
f$linkfun <- mli$linkfun
f$linkinv <- mli$linkinv
f$mu.eta <- mli$mu.eta
f$valideta <- mli$valideta

# initial residuals and leverages assuming constant dispersion
fitm <- glm.fit(X,y,family=f)
mu <- fitted(fitm)
d <- 2*( (y-mu)/mu - log(y/mu) )
p <- fitm$rank

# start from constant dispersion
phi <- -1/canonic.digamma(mean(d))*n/(n-p)
phi <- rep(phi,n)
fitd <- lm.fit(Z,dli$linkfun(phi))
gam <- ifelse(is.na(fitd$coef),0,fitd$coef)
if( mean(abs(fitd$residuals))/phi[1] > 1e-12 ) {
	# intercept is not in span of Z
	phi <- drop(dli$linkinv( Z %*% gam ))
	fitm <- glm.fit(X,y,weights=1/phi,mustart=mu,family=f)
	mu <- fitted(fitm)
	d <- 2*( (y-mu)/mu - log(y/mu) )
} else
	fitm <- glm.fit(X,y,weights=1/phi,mustart=mu,family=f)
dev <- const+sum(2*(lgamma(1/phi)+(1+log(phi))/phi)+d/phi)+const+2*log(prod(abs(diag(fitm$qr$qr)[1:p])))

# reml scoring
iter <- 0
if(trace) cat("Iter =",iter,", Dev =",format(dev,digits=13)," Gamma",gam,"\n")
Q2 <- array(0,c(n,p*(p+1)/2))
repeat {
	iter <- iter+1

	# gradient matrix
	eta <- dli$linkfun(phi)
	phidot <- dli$mu.eta(eta) * Z
	Z2 <- phidot / phi / sqrt(2)

	# information matrix and leverages
	Q <- qr.qy(fitm$qr, diag(1, nrow = n, ncol = p))
	j0 <- 0
	for(k in 0:(p-1)) {
		Q2[ ,(j0+1):(j0+p-k)] <- Q[ ,1:(p-k)] * Q[ ,(k+1):p]
		j0 <- j0+p-k
	}
	if(p>1) Q2[ ,(p+1):(p*(p+1)/2)] <- sqrt(2) * Q2[ ,(p+1):(p*(p+1)/2)]
	h <- drop( Q2[ ,1:p] %*% array(1,c(p,1)) )
	Q2Z <- crossprod(Q2,Z2)
	extradisp <- 2*( trigamma(1/phi) - trigamma(1/phi/h)/h )/phi^2 - (1-h)
	info <- crossprod(Z2,(extradisp+1-2*h)*Z2) + crossprod(Q2Z)

	# score vector
	deltah <- 2*(digamma(1/h/phi)+log(h)-digamma(1/phi))
	dl <- crossprod(phidot, (d - deltah)/(2*phi^2))

	# scoring step
	R <- chol(info)
	dgam <- backsolve(R,backsolve(R,dl,transpose=TRUE))
	gam <- gam + dgam

	# evaluate modified profile likelihood
	phi <- drop(dli$linkinv( Z %*% gam ))
	fitm <- glm.fit(X,y,weights=1/phi,mustart=mu,family=f)
	mu <- fitted(fitm)
	d <- 2*( (y-mu)/mu - log(y/mu) )
	dev <- const+sum(2*(lgamma(1/phi)+(1+log(phi))/phi)+d/phi)+const+2*log(prod(abs(diag(fitm$qr$qr)[1:p])))

	# iteration output
	if(trace) cat("Iter =",iter,", Dev =",format(dev,digits=13)," Gamma",gam,"\n")

	# test for convergence
	if( crossprod(dl,dgam) < tol ) break

	# test for iteration limit
	if(iter > maxit) {
		warning("Max iterations exceeded")
		break
	}
}

# Standard errors
se.gam <- sqrt(diag(chol2inv(chol(info))))
se.beta <- sqrt(diag(chol2inv(qr.R(fitm$qr))))

list(beta=fitm$coef,se.beta=se.beta,gamma=gam,se.gam=se.gam,mu=mu,phi=phi,deviance=dev,h=h)
}
#  SAGE.R

sage.test <- function(x, y, n1=sum(x), n2=sum(y))
#	Binomial probabilities for comparing SAGE libraries
#	Gordon Smyth
#	15 Nov 2003.  Last modified 21 Jan 2004.
{
	if(any(is.na(x)) || any(is.na(y))) stop("missing values not allowed")
	x <- round(x)
	y <- round(y)
	if(any(x<0) || any(y<0)) stop("x and y must be non-negative")
	if(length(x) != length(y)) stop("x and y must have same length")
	n1 <- round(n1)
	n2 <- round(n2)
	if(!missing(n1) && any(x>n1)) stop("x cannot be greater than n1")
	if(!missing(n2) && any(y>n2)) stop("y cannot be greater than n2")
	size <- x+y
	p.value <- rep(1,length(x))
	if(n1==n2) {
		i <- (size>0)
		if(any(i)) {
			x <- pmin(x[i],y[i])
			size <- size[i]
			p.value[i] <- pbinom(x,size=size,prob=0.5)+pbinom(size-x+0.5,size=size,prob=0.5,lower.tail=FALSE)
		}
		return(p.value)
	}
	prob <- n1/(n1+n2)
	if(any(big <- size>10000)) {
		ibig <- (1:length(x))[big]
		for (i in ibig) p.value[i] <- chisq.test(matrix(c(x[i],y[i],n1-x[i],n2-y[i]),2,2))$p.value
	}
	size0 <- size[size>0 & !big]
	if(length(size0)) for (isize in unique(size0)) {
		i <- (size==isize)
		p <- dbinom(0:isize,p=prob,size=isize)
		o <- order(p)
		cumsump <- cumsum(p[o])[order(o)]
		p.value[i] <- cumsump[x[i]+1]
	}
	p.value
}
##  TWEEDIEF.R

tweedie <- function(var.power=0, link.power=1-var.power) {
#	Tweedie generalized linear model family
#	Gordon Smyth
#	22 Oct 2002.  Last modified 20 Sep 2004

	lambda <- link.power
	if(lambda==0) {
		linkfun <- function(mu) log(mu)
		linkinv <- function(eta) pmax(.Machine$double.eps, exp(eta))
		mu.eta <- function(eta) pmax(.Machine$double.eps, exp(eta))
		valideta <- function(eta) TRUE
	} else {
		linkfun <- function(mu) mu^lambda
		linkinv <- function(eta) eta^(1/lambda)
		mu.eta <- function(eta) (1/lambda) * eta^(1/lambda - 1)
		valideta <- function(eta) TRUE
	}
	p <- var.power
	variance <- function(mu) mu^p
	if(p == 0)
		validmu <- function(mu) TRUE
	else if(p > 0)
		validmu <- function(mu) all(mu >= 0)
	else
		validmu <- function(mu) all(mu > 0)
	dev.resids <- function(y, mu, wt) {
		y1 <- y + (y == 0)
		if (p == 1)
			theta <- log(y1/mu)
		else
			theta <- ( y1^(1-p) - mu^(1-p) ) / (1-p)
		if (p == 2)
			kappa <- log(y1/mu)
		else
			kappa <- ( y^(2-p) - mu^(2-p) ) / (2-p)
		2 * wt * (y*theta - kappa)
	}	
	initialize <- expression({
		n <- rep(1, nobs)
		mustart <- y + 0.1 * (y == 0)
	})
	aic <- function(y, n, mu, wt, dev) NA
	structure(list(
		family = "Tweedie", variance = variance, dev.resids = dev.resids, aic = aic,
		link = paste("mu^",as.character(lambda),sep=""), linkfun = linkfun, linkinv = linkinv,
		mu.eta = mu.eta, initialize = initialize, validmu = validmu, valideta = valideta), 
		class = "family")
}

