R/BoxCoxmix.R In boxcoxmix: Box-Cox-Type Transformations for Linear and Logistic Models with Random Effects

Documented in bhatbhatfikfikmikmiknp.bhatnp.bhatnp.boxcoxnp.boxcoxnp.boxcoxmixnp.emnp.emnp.estepnp.estepnp.mstepnp.mstepnp.thetanp.zknp.zkvc.estepvc.estepvc.mstepvc.mstepvc.thetayhatyhatytransytranszk

# Box-Cox Transformation with Random Effects
#'
#' @rdname vc.em
#' @param p ..
#' @param beta ..
#' @param sigma ..
#' @param z ..
#' @export
np.estep <- function(y, x, lambda, p, beta, z, sigma){
n <- length(y)
K <- length(z)
w <-d <- s <- matrix(0, n,K)
f <- fik(y, x, lambda,  beta, z, sigma)
logpf <- t(apply(f, 1, "+", log(p)))
Mi <- apply(logpf, 1, max)
Sik <- logpf - Mi
ifelse(Sik < -760, 0, Sik)
eSik <- exp(Sik)
SeSik <- as.vector(apply(eSik, 1, sum))
w <- eSik/SeSik
return(w) }

##Zk
#'
#' @rdname vc.em
#' @export
np.zk <- function(y, x, w, beta, lambda){
n <- dim(w)[1]
K <- dim(w)[2]
z <- rep(0,K)
if  (all(!is.na(x))){
xbeta <- x%*%beta
} else {xbeta <- 0}
for(k in 1:K){
z[k] <- sum(w[,k]*(ytrans(y, lambda) - xbeta))/sum(w[,k])
}
return(z)
}
## fik
#'
#' @rdname vc.em
#' @export
fik <- function(y, x, lambda,  beta, z, sigma){
n <- length(y)
K <- length(z)
f <- matrix(0,n,K)
theta <- np.theta(y, x, lambda,  beta, z)
for(i in 1:n){
for(k in 1:K){
ytrx <- (1/(2*sigma^2))*theta[i,k]
d <- ((-1/2)*log(2*pi))
j <- log(sigma)
g <- (lambda -1)*log(y[i])
f[i,k] <- d-j-ytrx + g
}
}
return(f)
}
##theta
#' @rdname vc.em
#' @export
np.theta <- function(y, x, lambda,  beta, z){
n <- length(y)
K <- length(z)
theta <- matrix(0, n,K)
for(i in 1:n){
if  (all(!is.na(x))){
xb <- x[i,]%*%beta
} else { xb<- 0 }
for(k in 1:K){
theta[i,k] <- (ytrans(y[i],lambda) - xb - z[k])^2
}
}
return(theta)
}
##yhat
#'
#' @rdname vc.em
#' @param v ..
#' @export
yhat <-  function(v, lambda=1){
if (is.list(v)) { r<- length(v)
for (i in 1:r){
if(lambda == 0){ y[[i]] <- exp(v)
}
else {
y[[i]]<- (lambda*v +1)^(1/lambda)
}
}
} else {
if(lambda == 0){y <- exp(v)
}
else {
y <- (lambda*v + 1)^(1/lambda)
}
}
return(y)
}
#ytrans
#'
#' @rdname vc.em
#' @export
ytrans <-  function(y, lambda=1){
if (is.list(y)) { r<- length(y)
v <- list()
for (i in 1:r){
if(lambda == 0){v[[i]]<-  log(y[[i]])
}
else {
v[[i]] <-  (y[[i]]^lambda - 1)/ lambda
}
}
} else {
if(lambda == 0){v<-  log(y)
}
else {
v <-  (y^lambda - 1)/ lambda
}
}
return(v)
}
## bhat
#' @importFrom stats lm
#' @importFrom stats vcov
#' @rdname vc.em
#' @export
np.bhat <- function(y, x, w, z, lambda){
n <- dim(w)[1]
Ynew <- ytrans(y,lambda)- w%*%z
out <- lm(Ynew ~ -1 + x)
se<- sqrt(diag(vcov(out)))
beta<-  out$coef return(list("beta"=beta,"se"=se)) } ## M-step #' #' @rdname vc.em #' @export np.mstep<- function(y, x, beta, lambda, w){ n <- dim(w)[1] K <- dim(w)[2] p <- apply(w,2,sum)/n z <- rep(0,K) for (S in 1:40){ z <- np.zk(y, x, w, beta, lambda) if (all(!is.na(x))){ bse <- np.bhat(y, x, w, z, lambda) beta<-bse$beta
se<- bse$se } else { beta <- 0 se<-NA } } var <- 0 for(i in 1:n){ theta <- np.theta(y, x, lambda, beta, z) for(k in 1:K){ var<- var+ w[i,k]*theta[i,k] } } sigma <- sqrt(var/n) return(list("p"=p, "z"=z, "beta"=beta, "se"=se, "sigma"=sigma)) } ## EMstep #' @importFrom stats coef #' @importFrom stats sd #' @rdname vc.em #' @export np.em <- function(y, x, K, lambda=1, steps= 500,tol=0.5, start="gq", EMdev.change = 1e-04, plot.opt = 1,verbose = TRUE, ...){ n <- length(y) if (!is.matrix(x)){ x<- matrix(x,n,1) } p <- rep(1/K,K) yt <- ytrans(y, lambda) if (all(!is.na(x))){ a <- lm(formula = yt ~ -1 + x) se <- sqrt(diag(vcov(a))) beta <- coef(a) } else { beta<- 0 se <- NA } sigma <- sd(yt) tol0 <- tol lambda0 <- lambda if (K > 1 && start == "quantile") { z <- mean(yt)+tol*stats::quantile(yt-mean(yt), probs= (1:K)/K-1/(2*K)) } if (K > 1 && start == "gq") { if (K > 60) { K <- 60 cat("K was set equal to 60, since the number of mass points supported are only up to 60 mass points. \n") } if (all(!is.na(x))){ b <- lm(formula = yt ~ x) } else { b <- lm(formula = yt ~ 1) } Z <- b$coef[1] +tol*summary(b)$s*gqz(K, minweight = 1e-50)$location
z<-Z[order(Z)]
}
if (K == 1){
out <- lm(yt ~ x)
s <- 1
z <- out$coef[1] sizes <-"none" se<- sqrt(diag(vcov(out))) beta <- coef(out) w <- matrix(1, n,1) masses <- "none" sigma <- summary(out)$sigma
lik <- 0
for(i in 1:n){
lik <- lik + ((y[i]^(lambda -1))/(2*pi*sigma)^(n/2))*exp(((ytrans(y[i],lambda) - x[i]%*%beta )^2)/2*sigma^2)
}
logLik <- (-n/2)*log(2*pi)-n*log(sigma) - n/2 + (lambda -1)*sum(log(y))
Disp <- -2*logLik
complete.loglik <- "none"
converged <- "none"
Disparities<- "none"
} else {
s<-1
previous_loglik <- -(2^1000)
converged <- FALSE
loglik <- Disparities <- 0
masses<- matrix(0,0,K)
while (s <= steps && !converged){
if (verbose) {
cat(s, "..")
}
w   <- np.estep(y, x, lambda, p, beta, z, sigma)
fit <- np.mstep(y, x, beta, lambda, w)
p  <- fit$p z <- fit$z
beta <- fit$beta sigma <-fit$sigma
se <- fit$se lik <-matrix(0,n, K) f <- fik(y, x, lambda, beta, z, sigma) f <- exp(f) for(k in 1:K){ lik[,k] <- p[k]*f[,k] } loglik[s] <- sum(log(apply(lik,1,sum))) Disparities[s] <- -2*loglik[s] converged <- abs(loglik[s] - previous_loglik)< EMdev.change previous_loglik <- loglik[s] masses <- rbind(masses,z) logLik <- loglik[s] Disp <- Disparities[s] s<-s+1 } if (verbose) { cat("\n") if (converged) { cat("EM algorithm met convergence criteria at iteration # ", s - 1, "\n") } else { cat("EM algorithm failed to meet convergence criteria at iteration # ", s - 1, "\n") } } if (plot.opt==1){ graphics::par(mfrow=c(2,1),cex=0.5,cex.axis=1.5,cex.lab=1.5) } if ((plot.opt == 2 ) && graphics::par("mfrow")[1] > 2) { plot.main <- substitute("tol" == tol0, list(tol0 = tol0)) } else if ( (plot.opt == 3) && graphics::par("mfrow")[1] > 2 ) { plot.main <- substitute("lambda"== lambda0, list(lambda0 = lambda0)) } else { plot.main <- c("") } if (plot.opt == 1 || plot.opt == 2 || plot.opt == 3 ) { graphics::plot(1:(s - 1), Disparities, col = 1, type = "l", xlab = "EM iterations", ylab = "-2logL", main = plot.main) if (verbose) { cat("Disparity trend plotted.\n") } } complete.lik <-matrix(0,n, K) k <- 1:K complete.lik[,k] <- (lik[,k])^w[,k] complete.loglik <- sum(log(complete.lik)) } np.fit<- list("EMiteration"=s-1, "masses"=masses, "kind"=1,"p"=p, "mass.point"=z, "beta"=beta, "se"=se, "sigma"=sigma, "w" =w, "loglik"= logLik, "complete.loglik" = complete.loglik, "Disparities"=Disparities, "disparity"= Disp, "likelihood"= lik, "EMconverged" = converged) class(np.fit)<-"boxcoxmix" return(np.fit) } #np.boxcox #' @importFrom stats model.frame #' @importFrom stats model.matrix #' @importFrom stats model.response #' @rdname vc.em #' @export np.boxcox np.boxcox<- function(formula, groups=1, data, K = 3, tol = 0.5, lambda = 1, steps=500, EMdev.change = 1e-04, plot.opt = 1, verbose = TRUE, start="gq", ...){ call <- match.call() ddim <- dim(data) mf <- match.call(expand.dots = FALSE) m <- match(c("formula", "data"), names(mf), 0) mf <- mf[c(1, m)] mf$drop.unused.levels <- TRUE
mf[[1]] <- as.name("model.frame")
mf <- eval(mf, parent.frame())
mf <- model.frame(formula=formula, data=data)
model <- "mixed"
x <- model.matrix(attr(mf, "terms"), data = mf)[, -1, drop = FALSE]
if (dim(x)[2]== 0) {
x <- model.matrix(attr(mf, "terms"), data = mf)
model <- "pure"
}
y <- model.response(mf, "numeric")
if (any(y <= 0))
stop("response variable must be positive")
n <- NROW(y)
xnames <- dimnames(x)[[2]]
sizes <- groups
mform <- strsplit(as.character(sizes), "\\|")
mform <- gsub(" ", "", mform)
if (length(mform) >= 2) {
stop("Please use function vc.boxcox for two-level models")
}
npfit <- try(np.em(y=y, x=x, sizes= sizes,  K = K, lambda = lambda, steps = steps,tol = tol,
start = start, EMdev.change = EMdev.change,
plot.opt = plot.opt, verbose = verbose))
if (class(npfit) == "try-error") {
cat("Check model specification, change the number of components or specify another value of lambda or tol and try again.")
return()
}
s <- npfit$EMiteration w <- npfit$w
p  <- npfit$p names(p) <- paste('MASS',1:K,sep='') z <- npfit$mass.point
names(z) <- paste('MASS',1:K,sep='')
beta <- npfit$beta sigma <-npfit$sigma
se <- npfit$se disp <- npfit$disparity
Disparities<- npfit$Disparities if (model=='pure'){ aic<- disp+2*(2*K-1) bic<-disp+log(n)*(2*K-1) }else{ aic<- disp+2*(length(beta)+2*K-1) bic<-disp+log(n)*(length(beta)+2*K-1) } loglik <- npfit$loglik
complete.loglik <- npfit$complete.loglik masses <- npfit$masses
EMconverged <- npfit$EMconverged yt <- ytrans(y, lambda) if (K == 1){ylim<- "none"} if (model=='mixed'){ if (K == 1){ predicted.re <-rep(1,n) fitted.transformed <- rep(0,n) x <- model.matrix(attr(mf, "terms"), data=mf) for(i in 1:n){ fitted.transformed[i] <- x[i,]%*%beta } } else{ predicted.re <- rep(0,n) fitted.transformed <- rep(0,n) for(i in 1:n){ predicted.re[i] <- sum(w[i,]%*%z) fitted.transformed[i] <- x[i,]%*%beta + predicted.re[i] } } fitted <- yhat(fitted.transformed, lambda = lambda) residuals <- y- fitted residuals.transformed <- yt - fitted.transformed if (K > 1){ ylim <- c(min(masses[, ]), max(masses[, ])) if (plot.opt == 1) { graphics::plot(1:s, masses[,1], col=1, type = "l", ylim=ylim, ylab='mass points',xlab='EM iterations') for(k in 2:K){ graphics::lines(1:s, masses[,k], type="l", lwd= 1.5,lty=1, col=k)} if (verbose==T){cat("EM Trajectories plotted.\n")} } } npresult<- list("call"=call, "yt"=yt, "aic"=aic, "bic"=bic,"xx"=xnames, "Class"=y,"mform"=mform,"ylim"=ylim, "masses"=masses, "y"=y, "formula"=formula,"data"=data,"EMiteration"= s, "Disparities"=Disparities, "p"=p, "mass.point"=z, "beta"=beta, "se"=se, "sigma"=sigma, "w" =w, "loglik"= loglik, "complete.loglik" = complete.loglik, "disparity"= disp, "EMconverged" = EMconverged, "fitted" = fitted, "fitted.transformed"= fitted.transformed, "predicted.re"= predicted.re, "residuals"=residuals, "residuals.transformed"=residuals.transformed, "kind"=1, "model"=model) class(npresult)<-"boxcoxmix" } else { if (K > 1){ ylim <- c(min(masses[, ]), max(masses[, ])) if (plot.opt == 1) { graphics::plot(1:s, masses[,1], col=1, type = "l", ylim=ylim, ylab='mass points',xlab='EM iterations') for(k in 2:K){ graphics::lines(1:s, masses[,k], type="l", lwd= 1.5,lty=1, col=k)} if (verbose==T){cat("EM Trajectories plotted.\n")} } } fitted <- "none" fitted.transformed <- "none" predicted.re <- "none" npresult<- list("call"=call, "yt"=yt, "aic"=aic, "bic"=bic,"xx"=xnames, "Class"=y,"mform"=mform,"ylim"=ylim, "masses"=masses, "y"=y, "formula"=formula,"data"=data,"EMiteration"= s, "Disparities"=Disparities, "p"=p, "mass.point"=z, "beta"=beta, "se"=se, "sigma"=sigma, "w" =w, "loglik"= loglik, "complete.loglik" = complete.loglik, "disparity"= disp, "EMconverged" = EMconverged,"fitted" = fitted, "fitted.transformed"= fitted.transformed, "predicted.re"= predicted.re, "residuals"=y, "residuals.transformed"=yt, "kind"=1, "model"= model) class(npresult)<-"boxcoxmixpure" } return(npresult) } #Box-Cox Transformation with Variance Component model #' #' @rdname vc.em #' @export vc.estep <- function(Y, X, sizes=1, lambda, p, beta, z, sigma){ K <- length(z) r <- length(sizes) w <-d <- s <- matrix(0, r,K) m <- mik(Y, X, sizes, lambda, beta, z, sigma) logpm <- t(apply(m, 1, "+", log(p))) Mi <- apply(logpm, 1, max) Sik <- logpm - Mi ifelse(Sik < -760, 0, Sik) eSik <- exp(Sik) SeSik <- as.vector(apply(eSik, 1, sum)) w <- eSik/SeSik return(w) } ##Zk #' #' @rdname vc.em #' @export zk <- function(Y, X, sizes, w, beta, lambda){ r <- dim(w)[1] K <- dim(w)[2] z <- rep(0,K) yx <- rep(0,r) for(i in 1:r){ if (all(!is.na(X))){ Xbeta <- X[[i]]%*%beta } else {Xbeta <- 0} yx[i] <- sum(Y[[i]] - Xbeta)} for(k in 1:K){ z[k] <- sum(w[,k]*yx)/sum(sizes*w[,k]) } return(z) } ## bhat #' #' @rdname vc.em #' @param Y .. #' @param X .. #' @param w .. #' @export bhat <- function(Y, X, sizes, w, z, lambda){ r <- dim(w)[1] wz <- w%*%z WZ <- rep(wz, sizes) Xlong<-matrix(0, 0, dim(X[[1]])[2]) for(i in 1:r){ Xlong<-rbind(Xlong, X[[i]])} All<-as.data.frame(cbind(Ynew=unlist(Y)-WZ, Xlong)) out <- lm(Ynew ~ -1 + Xlong, data=All) se<- sqrt(diag(vcov(out))) beta<- out$coef
return(list("beta"=beta,"se"=se))
}

##mik
#'
#' @rdname vc.em
#' @export
mik <- function(Y, X, sizes, lambda,  beta, z, sigma){
K <- length(z)
r <- length(sizes)
m <- matrix(0,r,K)
theta <- vc.theta(Y, X, sizes, lambda,  beta, z)
for(i in 1:r){
for(k in 1:K){
m[i,k] <- (-sizes[i]/2)*log(2*pi)-sizes[i]*log(sigma)-(1/(2*sigma^2))*theta[i,k] + (lambda -1)*sum(log( Y[[i]]))
}
}
return(m)
}
##theta
#' @rdname vc.em
#' @export
vc.theta <- function(Y, X, sizes, lambda,  beta, z){
K <- length(z)
r <- length(sizes)
theta <- matrix(0, r,K)
Ytrans<-ytrans(Y, lambda)
for(i in 1:r){
if  (all(!is.na(X))){
Xb <- X[[i]]%*%beta
} else { Xb<- 0 }
for(k in 1:K){
theta[i,k] <- sum((Ytrans[[i]] - Xb - z[k])^2)
}
}
return(theta)
}
## M-step
#'
#' @rdname vc.em
#' @export
vc.mstep<- function(Y, X, sizes=1, beta, lambda, w){
r <- dim(w)[1]
K <- dim(w)[2]
p <- apply(w,2,sum)/r
z <- rep(0,K)
Ytrans <- ytrans(Y, lambda)
for (S in 1:40){
z <- zk(Ytrans, X, sizes, w, beta, lambda)
if  (all(!is.na(X))){
bse<-  bhat(Ytrans, X, sizes, w, z, lambda)
beta<-bse$beta se<- bse$se
} else {
beta <- 0
se<-NA
}
}
var <- 0
theta <- vc.theta(Y, X, sizes, lambda,  beta, z)
for(i in 1:r){
for(k in 1:K){
var<- var+ w[i,k]*theta[i,k]
}
}
sigma <- sqrt(var/sum(sizes))
return(list("p"=p, "z"=z, "beta"=beta, "se"=se, "sigma"=sigma))
}
## EMstep

#' Internal boxcoxmix functions
#'
#' @title Internal boxcoxmix functions
#' @description auxiliary functions are not intended to be directly called from the user.
#'
#' @rdname vc.em
#' @aliases np.boxcox vc.boxcox np.estep np.zk fik yhat ytrans np.bhat
#' nb.se np.mstep np.em vc.estep bhat vc.se mik vc.mstep vc.em gqz
#' @param y ..
#' @param x ..
#' @param sizes ..
#' @author Amani Almohaimeed and Jochen Einbeck
#' @keywords em
#' @export
vc.em <- function (y, x, sizes = 1, K, lambda, steps = 500, tol = 0.5,
start = "gq", EMdev.change = 1e-04, plot.opt = 1, verbose = TRUE,
...)
{
n <- length(y)
if (length(sizes) == 1) {
if (sizes == 1)
r <- n
else r <- 1
}
else {
r <- length(sizes)
}
if (!is.matrix(x)) {
x <- matrix(x, n, 1)
}
p <- rep(1/K, K)
yt <- ytrans(y, lambda)
if (all(!is.na(x))) {
a <- lm(formula = yt ~ -1 + x)
se <- sqrt(diag(vcov(a)))
beta <- coef(a)
}
else {
beta <- 0
se <- NA
}
sigma <- sd(yt)
tol0 <- tol
lambda0 <- lambda
if (K > 1 && start == "quantile") {
z <- mean(yt) + tol * stats::quantile(yt - mean(yt),
probs = (1:K)/K - 1/(2 * K))
}
if (K > 1 && start == "gq") {
if (K > 60) {
K <- 60
cat("K was set equal to 60, since the number of mass points supported are only up to 60 mass points. \n")
}
if (all(!is.na(x))) {
b <- lm(formula = yt ~ x)
}
else {
b <- lm(formula = yt ~ 1)
}
Z <- b$coef[1] + tol * summary(b)$s * gqz(K, minweight = 1e-50)$location z <- Z[order(Z)] } cumsize <- cumsum(sizes) Y <- Ytrans <- list() Y[[1]] <- y[1:sizes[1]] for(i in 2:r){ Y[[i]] <- y[(cumsize[i - 1] + 1):cumsize[i]]} Ytrans <- ytrans(Y, lambda) X <- list() X[[1]] <- x[1:sizes[1], ] X[[1]] <- matrix(X[[1]], ncol = length(beta)) for(i in 2:r){ X[[i]] <- x[(cumsize[i - 1] + 1):cumsize[i], ] X[[i]] <- matrix(X[[i]], ncol = length(beta)) } if (K == 1) { out <- lm(yt ~ x) s <- 1 z <- coef(out)[1] sizes <- "none" se <- sqrt(diag(vcov(out))) beta <- coef(out) w <- matrix(1, n, 1) masses <- "none" sigma <- summary(out)$sigma
lik <- 0
for(i in 1:n){
lik <- lik + ((y[i]^(lambda - 1))/(2 * pi * sigma)^(n/2)) *
exp(((ytrans(y[i], lambda) - x[i] %*% beta)^2)/2 *
sigma^2)}
logLik <- (-n/2) * log(2 * pi) - n * log(sigma) - n/2 +
(lambda - 1) * sum(log(y))
Disp <- -2 * logLik
complete.loglik <- "none"
converged <- "none"
Disparities<- "none"
}
else {
s <- 1
previous_loglik <- -(2^1000)
converged <- FALSE
loglik <- Disparities <- 0
masses <- matrix(0, 0, K)
while (s <= steps && !converged) {
if (verbose) {
cat(s, "..")
}
w <- vc.estep(Y, X, sizes, lambda, p, beta, z, sigma)
fit <- vc.mstep(Y, X, sizes, beta, lambda, w)
p <- fit$p z <- fit$z
beta <- fit$beta sigma <- fit$sigma
se <- fit$se lik <- matrix(0, r, K) m <- mik(Y, X, sizes, lambda, beta, z, sigma) m <- exp(m) for(k in 1:K){ lik[, k] <- p[k] * m[, k]} loglik[s] <- sum(log(apply(lik, 1, sum))) Disparities[s] <- -2 * loglik[s] converged <- abs(loglik[s] - previous_loglik) < EMdev.change previous_loglik <- loglik[s] masses <- rbind(masses, z) logLik <- loglik[s] Disp <- Disparities[s] s <- s + 1 } if (verbose) { cat("\n") if (converged) { cat("EM algorithm met convergence criteria at iteration # ", s - 1, "\n") } else { cat("EM algorithm failed to meet convergence criteria at iteration # ", s - 1, "\n") } } if (plot.opt == 1) { graphics::par(mfrow = c(2, 1), cex = 0.5, cex.axis = 1.5, cex.lab = 1.5) } if ((plot.opt == 2) && graphics::par("mfrow")[1] > 2) { plot.main <- substitute("tol" == tol0, list(tol0 = tol0)) } else if ((plot.opt == 3) && graphics::par("mfrow")[1] > 2) { plot.main <- substitute("lambda" == lambda0, list(lambda0 = lambda0)) } else { plot.main <- c("") } if (plot.opt == 1 || plot.opt == 2 || plot.opt == 3) { graphics::plot(1:(s - 1), Disparities, col = 1, type = "l", xlab = "EM iterations", ylab = "-2logL", main = plot.main) if (verbose) { cat("Disparity trend plotted.\n") } } complete.lik <- matrix(0, r, K) k <- 1:K complete.lik[, k] <- (lik[, k])^w[, k] complete.loglik <- sum(log(complete.lik)) } fit <- list(EMiteration = s - 1, masses = masses, kind = 1, Disparities = Disparities, p = p, mass.point = z, beta = beta, se = se, sigma = sigma, w = w, loglik = logLik, complete.loglik = complete.loglik, disparity = Disp, likelihood = lik, EMconverged = converged) class(fit) <- "boxcoxmix" return(fit) } # boxcoxmix #' Response Transformations for Random Effect and Variance Component Models #' #' #' #' @rdname np.boxcoxmix #' #' @description The function \code{np.boxcoxmix()} fits an overdispersed generalized linear model #' and variance component models using nonparametric profile maximum likelihood. #' #' #' @details The Box-Cox transformation (Box & Cox, 1964) is applied to the overdispersed #' generalized linear models and variance component models with an unspecified #' mixing distribution. The NPML estimate of the mixing distribution is known #' to be a discrete distribution involving a finite number of mass-points and corresponding #' masses (Aitkin et al., 2009). An Expectation-Maximization (EM) algorithm is #' used for fitting the finite mixture distribution, one needs to specify the #' number of components \code{K} of the finite mixture in advance. To stop the EM-algorithm #' when it reached its convergence point, we need to defined the convergence criteria that is #' the absolute change in the successive log-likelihood function values being less than an arbitrary #' parameter such as \code{EMdev.change} = 0.0001 (Einbeck et at., 2014). This algorithm can be #' implemented using the function \code{np.boxcoxmix()}, which is designed to account for #' overdispersed generalized linear models and variance component models using the non-parametric #' profile maximum likelihood (NPPML) estimation. #' #' #' The ability of the EM algorithm to locate the global maximum in fewer iterations #' can be affected by the choice of initial values, the function \code{np.boxcoxmix()} #' allows us to choose from two different methods to set the initial value of the mass #' points. When option "gq" is set, then Gauss-Hermite masses and mass points are used #' as starting points in the EM algorithm, while setting start= "quantile" uses the #' Quantile-based version to select the starting points. #' #' @aliases np.boxcoxmix #' @param formula a formula describing the transformed response and the fixed #' effect model (e.g. y ~ x). #' @param groups the random effects. To fit overdispersion models , set \code{groups} = 1. #' @param data a data frame containing variables used in the fixed and random #' effect models. #' @param K the number of mass points. #' @param tol a positive scalar (usually, 0< \code{tol} <= 2) #' @param lambda a transformation parameter, setting \code{lambda}=1 means 'no #' transformation'. #' @param steps maximum number of iterations for the EM algorithm. #' @param EMdev.change a small scalar, with default 0.0001, used to determine #' when to stop EM algorithm. #' @param plot.opt Set \code{plot.opt=1}, to plot the disparity against #' iteration number. Use \code{plot.opt=2} for \code{tolfind.boxcox()} and \code{plot.opt=3} #' for \code{optim.boxcox()}. #' @param verbose If set to FALSE, no printed output on progress. #' @param start a description of the initial values to be used in the fitted #' model, Quantile-based version "quantile" or Gaussian Quadrature "gq" can be #' set. #' @param \dots extra arguments will be ignored. #' @return #' \item{mass.point}{the fitted mass points.} \item{p}{the masses corresponding to the mixing proportions.} \item{beta}{the #' vector of coefficients.} \item{sigma}{the standard deviation of the mixing distribution (the square root of the variance).} \item{se}{the standard error of the estimate.} #' \item{w}{a matrix of posterior probabilities that element i comes from cluster k.} #' \item{loglik}{the log-likelihood of the fitted regression model.} #' \item{complete.loglik}{the complete log-likelihood of the fitted regression model.} #' \item{disparity}{the disparity of the fitted regression model.} \item{EMiteration}{provides the number of iterations of the EM algorithm.} #' \item{EMconverged}{TRUE means the EM algorithm converged.} \item{call}{the matched call.} \item{formula}{the formula provided.} #' \item{data}{the data argument.} \item{aic}{the Akaike information criterion of the fitted regression model.} #' \item{bic}{the Bayesian information criterion of the fitted regression model.} #' \item{fitted}{the fitted values for the individual observations.} \item{fitted.transformed}{the fitted values for #' the individual transformed observations.} \item{residuals}{the difference between the observed values and the fitted values.} #' \item{residuals.transformed}{the difference between the transformed observed values and the transformed fitted values.} #' \item{predicted.re}{a vector of predicted residuals.} #' #' The other outcomes are not relevant to users and they are intended for internal use only. #' #' @author Amani Almohaimeed and Jochen Einbeck #' @seealso \code{\link{optim.boxcox}}, \code{\link{tolfind.boxcox}}. #' @references #' Box G. and Cox D. (1964). An analysis of transformations. Journal of #' the Royal Statistical Society. Series B (Methodological), pages 211-252. #' #' Aitkin, M. A., Francis, B., Hinde, J., and Darnell, R. (2009). Statistical #' modelling in R. Oxford University Press Oxford. #' #' Jochen Einbeck, Ross Darnell and John Hinde (2014). npmlreg: #' Nonparametric maximum likelihood estimation for random effect #' models. R package version 0.46-1. #' #' @keywords boxcox random variance #' @examples #' #Pennsylvanian Hospital Stay Data #' data(hosp, package = "npmlreg") #' test1 <- np.boxcoxmix(duration ~ age + wbc1, data = hosp, K = 2, tol = 1, #' start = "quantile", lambda = 1) #' round(summary(test1)$w, digits = 3)
#' # [1,] 1.000 0.000
#'
#' # Refinery yield of gasoline Data
#' data(Gasoline, package = "nlme")
#' test2.vc <- np.boxcoxmix(yield ~ endpoint + vapor, groups = Gasoline$Sample, #' data = Gasoline, K = 3, tol = 1.7, start = "quantile", lambda = 0) #' test2.vc$disparity
#' # [1] 176.9827
#'
#'
#'
#'
#'
#'
#'
#'
# vc.boxcox
#' @export np.boxcoxmix
np.boxcoxmix<-function(formula, groups = 1, data, K = 3, tol = 0.5, lambda = 1,
steps = 500, EMdev.change = 1e-04, plot.opt = 1, verbose = TRUE,
start="gq", ...){
call <- match.call()
mform <- strsplit(as.character(groups), "\\|")
mform <- gsub(" ", "", mform)
if (length(mform) == 1) {
fit <- try(np.boxcox(formula=formula, groups= groups, data=data, K=K, lambda=lambda,steps= steps,
tol = tol, start=start, EMdev.change = EMdev.change, plot.opt = plot.opt, verbose = verbose))
if (class(fit) == "try-error") {
cat("Check model specification, change the number of components or specify another value of lambda or tol and try again.")
return()
}
}
else {
fit <- try(vc.boxcox(formula=formula, groups= groups, data=data, K=K, lambda=lambda, steps= steps,
tol = tol, start=start, EMdev.change = EMdev.change, plot.opt = plot.opt, verbose = verbose))
if (class(fit) == "try-error") {
cat("Check model specification, change the number of components or specify another value of lambda or tol and try again.")
return()
}
}
W <- fit$w P <- fit$p
se <- fit$se iter <- fit$EMiteration
names(P) <- paste('MASS',1:K,sep='')
Z <- fit$mass.point names(Z) <- paste('MASS',1:K,sep='') Beta <- fit$beta
Sigma<- fit$sigma aic<- fit$aic
bic<- fit$bic y <- fit$y
yt <- fit$yt fitted <- fit$fitted
fitted.transformed <- fit$fitted.transformed masses<- fit$masses
ylim<- fit$ylim residuals<- fit$residuals
residuals.transformed<- fit$residuals.transformed predicted.re <- fit$predicted.re
Class<-fit$Class xx <- fit$xx
model<- fit$model Disp <- fit$disparity
Disparities <- fit$Disparities Loglik <- fit$loglik
complete.loglik <- fit$complete.loglik kind <- fit$kind
EMconverged <- fit$EMconverged model <- fit$model
if (model == "mixed") {
result<- list("call"=call,  "p"=P, "mass.point"=Z, "beta"=Beta, "sigma"=Sigma, "se"=se, "w" =W, "Disparities" = Disparities,
"formula" = formula, "data" = data, "loglik"= Loglik, "aic"=aic, "bic"=bic,"masses"=masses, "y"=y, "yt"=yt,
"complete.loglik" = complete.loglik,
"disparity"= Disp, "EMconverged" = EMconverged, "mform"=length(mform),"ylim"=ylim,
"fitted" = fitted, "Class"= Class, "fitted.transformed"= fitted.transformed, "predicted.re"= predicted.re,
"residuals"=residuals, "residuals.transformed"=residuals.transformed,"kind"=kind,
"EMiteration"= iter, "xx" = xx, model= model)
class(result)<-"boxcoxmix"
}else{
result<- list("call"=call,  "p"=P, "mass.point"=Z, "beta"=Z, "sigma"=Sigma, "se"=se, "w" =W, "Disparities" = Disparities,
"formula" = formula, "data" = data, "loglik"= Loglik, "aic"=aic, "bic"=bic, "masses"=masses, "y"=y, "yt"=yt,
"complete.loglik" = complete.loglik,
"disparity"= Disp, "EMconverged" = EMconverged, "mform"=length(mform),"ylim"=ylim,
"fitted" = fitted, "Class"= Class, "fitted.transformed"= fitted.transformed, "predicted.re"= predicted.re,
"residuals"=residuals, "residuals.transformed"=residuals.transformed,"kind"=kind,
"EMiteration"= iter, "xx" = xx, model= model)
class(result)<-"boxcoxmixpure"
}
return(result)
}
#' @rdname vc.em
#' @param formula a formula describing the transformed response and the fixed
#' effect model (e.g. y ~ x).
#' @param groups the random effects. To fit overdispersion models , set \code{groups} = 1.
#' @param data a data frame containing variables used in the fixed and random
#' effect models.
#' @param K the number of mass points.
#' @param tol a positive scalar (usually, 0< \code{tol} <= 2)
#' @param lambda a transformation parameter, setting \code{lambda}=1 means 'no
#' transformation'.
#' @param steps maximum number of iterations for the EM algorithm.
#' @param EMdev.change a small scalar, with default 0.0001, used to determine
#' when to stop EM algorithm.
#' @param plot.opt Set plot.opt=1, to plot the disparity against
#' iteration number. Use \code{plot.opt=2} for \code{tolfind.boxcox} and \code{plot.opt=3}
#' for \code{optim.boxcox}.
#' @param verbose If set to FALSE, no printed output on progress.
#' @param start a description of the initial values to be used in the fitted
#' model, Quantile-based version "quantile" or Gaussian Quadrature "gq" can be
#' set.
#' @param \dots extra arguments will be ignored.
#' @export
vc.boxcox<- function (formula, groups = 1, data, K = 3, tol = 0.5, lambda = 1,
steps = 500, EMdev.change = 1e-04, plot.opt = 1, verbose = TRUE,
start="gq", ...)
{
call <- match.call()
data <- as.data.frame(data)
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "data"), names(mf), 0)
mf <- mf[c(1, m)]
mf$drop.unused.levels <- TRUE mf[[1]] <- as.name("model.frame") mf <- eval(mf, parent.frame()) datay <- data groupy <- groups ordered <- data[with(data, order(groups)), ] mf <- model.frame(formula = formula, data = ordered) model <- "mixed" x <- model.matrix(attr(mf, "terms"), data = mf)[, -1, drop = FALSE] if (dim(x)[2]== 0) { x <- model.matrix(attr(mf, "terms"), data = mf) model <- "pure" } y <- model.response(mf) if (any(y <= 0)) stop("response variable must be positive") N <- NROW(y) data <- if (is.matrix(y)) data[dimnames(y)[[1]], ] else data[names(y), ] xnames <- dimnames(x)[[2]] groups <- groups[match(rownames(ordered), rownames(datay))] sizes <- table(groups) mform <- strsplit(as.character(sizes), "\\|") mform <- gsub(" ", "", mform) if (length(mform) == 1) { stop("Please use function np.boxcox for overdispersion models") } vfit <- try(vc.em(y = y, x = x, sizes = sizes, K = K, lambda = lambda, steps = steps, tol = tol, start = start, EMdev.change = EMdev.change, plot.opt = plot.opt, verbose = verbose)) if (class(vfit) == "try-error") { cat("Check model specification, change the number of components or specify another value of lambda or tol and try again.") return() } EMiter <- vfit$EMiteration
if (K == 1) {
w <- vfit$w } else { w <- vfit$w
w <- w[match(unique(groupy), unique(groups)), ]
rownames(w) <- unique(groupy)
}
p <- vfit$p names(p) <- paste("MASS", 1:K, sep = "") z <- vfit$mass.point
names(z) <- paste("MASS", 1:K, sep = "")
beta <- vfit$beta sigma <- vfit$sigma
se <- vfit$se disp <- vfit$disparity
Disparities <- vfit$Disparities loglik <- vfit$loglik
EMconverged <- vfit$EMconverged complete.loglik <- vfit$complete.loglik
masses <- vfit$masses if (model=='pure'){ aic<- disp+2*(2*K-1) bic<-disp+log(N)*(2*K-1) }else{ aic<- disp+2*(length(beta)+2*K-1) bic<-disp+log(N)*(length(beta)+2*K-1) } yt <- ytrans(y, lambda) if (K == 1) { ylim <- "none" } if (model == "mixed") { if (K == 1) { predicted.re <- rep(1, N) } else { r <- dim(w)[1] predicted <- rep(0, r) predicted[1] <- sum(w[1, ] %*% z) i <- 2:r predicted[i] <- sum(w[i, ] %*% z) names(predicted) <- unique(groupy) predicted.re <- predicted } if (K == 1) { fitted.transformed <- rep(0, N) x <- model.matrix(attr(mf, "terms"), data = mf) for(i in 1:N){ fitted.transformed[i] <- x[i, ] %*% beta } } else { predicted <- predicted[match(unique(groups), unique(groupy))] fitted.transformed <- matrix(0, N, 1) cumsize <- cumsum(sizes) for(j in 1:cumsize[1]){ fitted.transformed[j, ] <- x[j, ] %*% beta + predicted[1] } for(i in 2:r){ for(j in (cumsize[i - 1] + 1):cumsize[i]){ fitted.transformed[j, ] <- x[j, ] %*% beta + predicted[i] } } rownames(fitted.transformed) <- rownames(ordered) fitted.transformed <- fitted.transformed[match(row.names(datay), row.names(ordered)), 1, drop = FALSE] } fitted <- yhat(fitted.transformed, lambda = lambda) residuals <- y - fitted residuals.transformed <- yt - fitted.transformed if (K > 1) { ylim <- c(min(masses[, ]), max(masses[, ])) if (plot.opt == 1) { graphics::plot(1:EMiter, masses[, 1], col = 1, type = "l", ylim = ylim, ylab = "mass points", xlab = "EM iterations") for(k in 2:K){ graphics::lines(1:EMiter, masses[, k], type = "l", lwd = 1.5, lty = 1, col = k)} if (verbose == T) { cat("EM Trajectories plotted.\n") } } } result <- list(call = call, aic = aic, bic = bic, xx = xnames, yt = yt, Disparities = Disparities, Class = sizes, mform = mform, ylim = ylim, masses = masses, y = y, formula = formula, data = data, EMiteration = EMiter, p = p, mass.point = z, beta = beta, se = se, sigma = sigma, w = w, loglik = loglik, complete.loglik = complete.loglik, disparity = disp, EMconverged = EMconverged, fitted = fitted, fitted.transformed = fitted.transformed, predicted.re = predicted.re, residuals = residuals, residuals.transformed = residuals.transformed, kind = 1,model=model) class(result) <- "boxcoxmix" } else { if (K > 1) { ylim <- c(min(masses[, ]), max(masses[, ])) if (plot.opt == 1) { graphics::plot(1:EMiter, masses[, 1], col = 1, type = "l", ylim = ylim, ylab = "mass points", xlab = "EM iterations") for(k in 2:K){ graphics::lines(1:EMiter, masses[, k], type = "l", lwd = 1.5, lty = 1, col = k)} if (verbose == T) { cat("EM Trajectories plotted.\n") } } } fitted <- "none" fitted.transformed <- "none" predicted.re <- "none" result <- list(call = call, aic = aic, bic = bic, xx = xnames, yt = yt, Disparities = Disparities, Class = sizes, mform = mform, ylim = ylim, masses = masses, y = y, formula = formula, data = data, EMiteration = EMiter, p = p, mass.point = z, beta = beta, se = se, sigma = sigma, w = w, loglik = loglik, complete.loglik = complete.loglik, fitted = fitted, fitted.transformed = fitted.transformed, predicted.re = predicted.re, disparity = disp, EMconverged = EMconverged, residuals = y, residuals.transformed = yt, kind = 1, model= model) class(result) <- "boxcoxmixpure" } return(result) } ##gqz #' @import "statmod" #' @rdname vc.em #' @param numnodes .. #' @param minweight .. #' @export gqz <- function (numnodes = 20, minweight = 1e-06) { out <- statmod::gauss.quad(numnodes, "hermite") h <- rbind(out$nodes * sqrt(2), out$weights/sum(out$weights))
ord <- order(h[1, ], decreasing = TRUE)
h <- h[, ord]
h <- cbind(h[1, ], h[2, ])
h <- subset(as.data.frame(h), (h[, 2] >= minweight))
names(h) <- c("location", "weight")
h
}


Try the boxcoxmix package in your browser

Any scripts or data that you put into this service are public.

boxcoxmix documentation built on May 2, 2019, 5:42 a.m.