Title: | Multivariate Response Generalized Linear Models |
---|---|
Description: | Provides functions that (1) fit multivariate discrete distributions, (2) generate random numbers from multivariate discrete distributions, and (3) run regression and penalized regression on the multivariate categorical response data. Implemented models include: multinomial logit model, Dirichlet multinomial model, generalized Dirichlet multinomial model, and negative multinomial model. Making the best of the minorization-maximization (MM) algorithm and Newton-Raphson method, we derive and implement stable and efficient algorithms to find the maximum likelihood estimates. On a multi-core machine, multi-threading is supported. |
Authors: | Yiwen Zhang <[email protected]> and Hua Zhou <[email protected]> |
Maintainer: | Juhyun Kim <[email protected]> |
License: | GPL (>= 2) |
Version: | 0.2.1 |
Built: | 2024-08-22 06:28:35 UTC |
Source: | CRAN |
The package provides functions that (1) fit multivariate discrete distributions, (2) generate random numbers from multivariate discrete distributions, and (3) run regression and penalized regression on the multivariate categorical response data. Implemented models include: multinomial logit model, Dirichlet multinomial model, generalized Dirichlet multinomial model, and negative multinomial model. Making the best of the minorization-maximization (MM) algorithm and Newton-Raphson method, we derive and implement stable and efficient algorithms to find the maximum likelihood estimates. On a multi-core machine, multi-threading is supported.
Package: | MGLM |
Type: | Package |
Version: | 0.0.9 |
Date: | 2017-12-14 |
License: | GPL (>= 2) |
Depends: | R (>= 3.0.0), methods, stats, parallel |
Yiwen Zhang and Hua Zhou
Calculates the Akaike's information criterion (AIC) for a fitted model object.
## S4 method for signature 'MGLMfit' AIC(object) ## S4 method for signature 'MGLMreg' AIC(object) ## S4 method for signature 'MGLMsparsereg' AIC(object) ## S4 method for signature 'MGLMtune' AIC(object)
## S4 method for signature 'MGLMfit' AIC(object) ## S4 method for signature 'MGLMreg' AIC(object) ## S4 method for signature 'MGLMsparsereg' AIC(object) ## S4 method for signature 'MGLMtune' AIC(object)
object |
MGLM object. |
Returns a numeric value with the corresponding AIC.
For the class "MGLMtune"
, the function returns AIC
based on the optimal tuning parameter.
set.seed(124) n <- 200 d <- 4 alpha <- rep(1, d-1) beta <- rep(1, d-1) m <- 50 Y <- rgdirmn(n, m, alpha, beta) gdmFit <- MGLMfit(Y, dist="GDM") AIC(gdmFit)
set.seed(124) n <- 200 d <- 4 alpha <- rep(1, d-1) beta <- rep(1, d-1) m <- 50 Y <- rgdirmn(n, m, alpha, beta) gdmFit <- MGLMfit(Y, dist="GDM") AIC(gdmFit)
Calculates the Bayesian information criterion (BIC) for a fitted model object.
## S4 method for signature 'MGLMfit' BIC(object) ## S4 method for signature 'MGLMreg' BIC(object) ## S4 method for signature 'MGLMsparsereg' BIC(object) ## S4 method for signature 'MGLMtune' BIC(object)
## S4 method for signature 'MGLMfit' BIC(object) ## S4 method for signature 'MGLMreg' BIC(object) ## S4 method for signature 'MGLMsparsereg' BIC(object) ## S4 method for signature 'MGLMtune' BIC(object)
object |
MGLM object. |
Returns a numeric value with the corresponding BIC.
For the class "MGLMtune"
, the function returns BIC
based on the optimal tuning parameter.
set.seed(124) n <- 200 d <- 4 alpha <- rep(1, d-1) beta <- rep(1, d-1) m <- 50 Y <- rgdirmn(n, m, alpha, beta) gdmFit <- MGLMfit(Y, dist="GDM") BIC(gdmFit)
set.seed(124) n <- 200 d <- 4 alpha <- rep(1, d-1) beta <- rep(1, d-1) m <- 50 Y <- rgdirmn(n, m, alpha, beta) gdmFit <- MGLMfit(Y, dist="GDM") BIC(gdmFit)
coef
extracts estimated model coefficients of class. coefficients
is an alias for it.
## S4 method for signature 'MGLMfit' coef(object) ## S4 method for signature 'MGLMreg' coef(object) ## S4 method for signature 'MGLMsparsereg' coef(object) ## S4 method for signature 'MGLMtune' coef(object)
## S4 method for signature 'MGLMfit' coef(object) ## S4 method for signature 'MGLMreg' coef(object) ## S4 method for signature 'MGLMsparsereg' coef(object) ## S4 method for signature 'MGLMtune' coef(object)
object |
an object for which the extraction of model coefficients is meaningful.
One of the following classes |
Method coef.
Coefficients extracted from the model object object
.
For the class "MGLMtune"
, the function returns model coefficients
based on the optimal tuning parameter.
library("MGLM") data("rnaseq") data <- rnaseq[, 1:6] mnreg <- MGLMreg(formula = cbind(X1, X2, X3, X4, X5, X6) ~ log(totalReads) + treatment + age + gender, data = rnaseq, dist = "MN") coef(mnreg)
library("MGLM") data("rnaseq") data <- rnaseq[, 1:6] mnreg <- MGLMreg(formula = cbind(X1, X2, X3, X4, X5, X6) ~ log(totalReads) + treatment + age + gender, data = rnaseq, dist = "MN") coef(mnreg)
An object that specifies the distribution to be fitted by the MGLMfit
function, or the regression model to be fitted by the MGLMreg
or MGLMsparsereg
functions.
Can be chosen from "MN"
, "DM"
, "NegMN"
, or "GDM"
.
A multinomial distribution models the counts of possible outcomes.
The counts of categories are negatively correlated.
The density of a
category count vector
with parameter
is
where ,
, and
.
Here,
, often read as "
choose
", refers the number of
combinations from a set of
elements.
The MGLMreg
function with dist="MN"
calculates the MLE of regression coefficients of the multinomial logit model, which has link function
,
. The
MGLMsparsereg
function with dist="MN"
fits regularized multinomial logit model.
When the multivariate count data exhibits over-dispersion, the traditional
multinomial model is insufficient. Dirichlet multinomial distribution models the
probabilities of the categories by a Dirichlet distribution.
The density of a category count vector
, with
parameter
,
, is
where . Here,
, often read as "
choose
",
refers the number of
combinations from a set of
elements.
The MGLMfit
function with dist="DM"
calculates the maximum likelihood estimate (MLE) of . The
MGLMreg
function with dist="DM"
calculates the MLE of regression coefficients of the Dirichlet multinomial regression model, which has link function
,
. The
MGLMsparsereg
function with dist="DM"
fits regularized Dirichlet multinomial regression model.
The more flexible Generalized Dirichlet multinomial model can be used when the counts of categories have both positive and negative correlations.
The probability mass of a count vector over
trials with parameter
,
, is
where and
. Here,
, often read as "
choose
",
#' refers the number of
combinations from a set of
elements.
The MGLMfit
with dist="GDM"
calculates the MLE of . The
MGLMreg
function with dist="GDM"
calculates the MLE of regression coefficients of the generalized Dirichlet multinomial regression model, which has link functions
and
,
. The
MGLMsparsereg
function with dist="GDM"
fits regularized generalized Dirichlet multinomial regression model.
Both the multinomial distribution and Dirichlet multinomial distribution are good for
negatively correlated counts. When the counts of categories are positively
correlated, the negative multinomial distribution is preferred.
The probability mass function of a category count vector
with parameter
,
,
,
, is
where . Here,
, often read as "
choose
", refers the number of
combinations from a set of
elements.
The MGLMfit
function with dist="NegMN"
calculates the MLE of . The
MGLMreg
function with dist="NegMN"
and regBeta=FALSE
calculates the MLE of regression coefficients of the negative multinomial regression model, which has link function
,
,
. When
dist="NegMN"
and regBeta=TRUE
, the overdispersion parameter is linked to covariates via , and the
function
MGLMreg
outputs an estimated matrix of
. The
MGLMsparsereg
function with dist="NegMN"
fits regularized negative multinomial regression model.
Yiwen Zhang and Hua Zhou
MGLMfit
, MGLMreg
, MGLMsparsereg
,
dmn
, ddirmn
, dgdirmn
, dnegmn
Fit the specified multivariate discrete distribution.
DMD.DM.fit( data, init, weight, epsilon = 1e-08, maxiters = 150, display = FALSE ) DMD.GDM.fit( data, init, weight, epsilon = 1e-08, maxiters = 150, display = FALSE ) DMD.NegMN.fit( data, init, weight, epsilon = 1e-08, maxiters = 150, display = FALSE ) MGLMfit( data, dist, init, weight, epsilon = 1e-08, maxiters = 150, display = FALSE )
DMD.DM.fit( data, init, weight, epsilon = 1e-08, maxiters = 150, display = FALSE ) DMD.GDM.fit( data, init, weight, epsilon = 1e-08, maxiters = 150, display = FALSE ) DMD.NegMN.fit( data, init, weight, epsilon = 1e-08, maxiters = 150, display = FALSE ) MGLMfit( data, dist, init, weight, epsilon = 1e-08, maxiters = 150, display = FALSE )
data |
a data frame or matrix containing the count data. Rows of the matrix represent observations and columns are the categories. Rows and columns of all zeros are automatically removed. |
init |
an optional vector of initial value of the parameter estimates. Should have the same dimension as the estimated parameters. See |
weight |
an optional vector of weights assigned to each row of the data. Should be Null or a numeric vector with the length equal to the number of rows of |
epsilon |
an optional numeric controlling the stopping criterion. The algorithm terminates when the relative change in the log-likelihoods of two successive iterates is less than |
maxiters |
an optional number controlling the maximum number of iterations. The default value is |
display |
an optional logical variable controlling the display of iterations. The default value is FALSE. |
dist |
a description of the distribution to fit. Choose from |
See dist
for details about model parameterization.
Returns an object of S4 class "MGLMfit"
. An object of class "MGLMfit"
is a list containing at least the following components:
estimate
the vector of the distribution prameter estimates.
SE
the vector of standard errors of the estimates.
vcov
the variance-covariance matrix of the estimates.
logL
the loglikelihood value.
iter
the number of iterations used.
BIC
Bayesian information criterion.
AIC
Akaike information criterion.
distribution
the distribution fitted.
LRT
when dist="DM"
or "GDM"
, it is the likelihood ratio test statistic for comparing the current model to the multinomial model. No LRT provided when dist="NegMN"
.
LRTpvalue
the likelihood ratio test P value.
gradient
the gradient at the estimated parameter values.
DoF
the degrees of freedom of the model.
Yiwen Zhang and Hua Zhou
data(rnaseq) Y <- as.matrix(rnaseq[, 1:6]) fit <- MGLMfit(data=Y, dist="GDM")
data(rnaseq) Y <- as.matrix(rnaseq[, 1:6]) fit <- MGLMfit(data=Y, dist="GDM")
dof
extracts the degrees of freedom of the estimated parameter
from the object of class MGLMsparsereg
.
## S4 method for signature 'MGLMsparsereg' dof(object)
## S4 method for signature 'MGLMsparsereg' dof(object)
object |
an object of class |
Returns degrees of freedom of object
.
library("MGLM") dist <- "DM" n <- 100 p <- 10 d <- 5 set.seed(118) m <- rbinom(n, 200, 0.8) X <- matrix(rnorm(n * p), n, p) alpha <- matrix(0, p, d) alpha[c(1, 3, 5), ] <- 1 Alpha <- exp(X %*% alpha) Y <- rdirmn(size = m, alpha = Alpha) pen <- "group" ngridpt <- 30 spmodelfit <- MGLMsparsereg(formula = Y ~ 0 + X, dist = dist, lambda = Inf, penalty = pen) df <- dof(spmodelfit)
library("MGLM") dist <- "DM" n <- 100 p <- 10 d <- 5 set.seed(118) m <- rbinom(n, 200, 0.8) X <- matrix(rnorm(n * p), n, p) alpha <- matrix(0, p, d) alpha[c(1, 3, 5), ] <- 1 Alpha <- exp(X %*% alpha) Y <- rdirmn(size = m, alpha = Alpha) pen <- "group" ngridpt <- 30 spmodelfit <- MGLMsparsereg(formula = Y ~ 0 + X, dist = dist, lambda = Inf, penalty = pen) df <- dof(spmodelfit)
Return the Khatri-Rao product of two matrices, which is a column-wise Kronecker product.
kr(A, B, w, byrow = TRUE)
kr(A, B, w, byrow = TRUE)
A , B
|
matrices. The two matrices |
w |
the weights vector. The length of the vector should match with the dimension of the matrices. If performing column-wise Kronecker product, the length of w should be the same as the column number of A and B. If performing row-wise Kronecker prodoct, the length of w should be the same as the row number of A and B. The default is a vector of 1 if no value provided. |
byrow |
a logical variable controlling whether to perform row/column-wise Kronecker product.
The default is |
The column/row-wise Kronecker product.
A matrix of the Khatri-Rao product.
Yiwen Zhang and Hua Zhou
X <- matrix(rnorm(30), 10, 3) Y <- matrix(runif(50), 10, 5) C <- kr(X, Y)
X <- matrix(rnorm(30), 10, 3) Y <- matrix(runif(50), 10, 5) C <- kr(X, Y)
logLik
extracts log-likelihood for classes "MGLMfit"
,
"MGLMreg"
, "MGLMsparsereg"
.
## S4 method for signature 'MGLMfit' logLik(object) ## S4 method for signature 'MGLMreg' logLik(object) ## S4 method for signature 'MGLMsparsereg' logLik(object)
## S4 method for signature 'MGLMfit' logLik(object) ## S4 method for signature 'MGLMreg' logLik(object) ## S4 method for signature 'MGLMsparsereg' logLik(object)
object |
an object from which a log-likelihood value can be extracted. |
Returns a log-likelihood value of object
.
library("MGLM") data("rnaseq") data <- rnaseq[, 1:6] dmFit <- MGLMfit(data, dist = "DM") logLik(dmFit)
library("MGLM") data("rnaseq") data <- rnaseq[, 1:6] dmFit <- MGLMfit(data, dist = "DM") logLik(dmFit)
maxlambda
extracts the maximum tuning parameter that ensures
the estimated regression coefficients are not all zero for the object of class MGLMsparsereg
.
## S4 method for signature 'MGLMsparsereg' maxlambda(object)
## S4 method for signature 'MGLMsparsereg' maxlambda(object)
object |
an object of class |
Returns a maximum lambda value of object
.
library("MGLM") dist <- "DM" n <- 100 p <- 10 d <- 5 set.seed(118) m <- rbinom(n, 200, 0.8) X <- matrix(rnorm(n * p), n, p) alpha <- matrix(0, p, d) alpha[c(1, 3, 5), ] <- 1 Alpha <- exp(X %*% alpha) Y <- rdirmn(size = m, alpha = Alpha) pen <- "group" ngridpt <- 30 spmodelfit <- MGLMsparsereg(formula = Y ~ 0 + X, dist = dist, lambda = Inf, penalty = pen) maxlambda <- maxlambda(spmodelfit)
library("MGLM") dist <- "DM" n <- 100 p <- 10 d <- 5 set.seed(118) m <- rbinom(n, 200, 0.8) X <- matrix(rnorm(n * p), n, p) alpha <- matrix(0, p, d) alpha[c(1, 3, 5), ] <- 1 Alpha <- exp(X %*% alpha) Y <- rdirmn(size = m, alpha = Alpha) pen <- "group" ngridpt <- 30 spmodelfit <- MGLMsparsereg(formula = Y ~ 0 + X, dist = dist, lambda = Inf, penalty = pen) maxlambda <- maxlambda(spmodelfit)
These functions are provided for compatibility with older version of the yourPackageName package. They may eventually be completely removed.
ddirm(...) rdirm(...) dgdirm(...) rgdirm(...) dneg(Y, alpha, beta)
ddirm(...) rdirm(...) dgdirm(...) rgdirm(...) dneg(Y, alpha, beta)
... |
parameters to be passed to the modern version of the function |
Y , alpha , beta
|
for functions |
ddirm |
now a synonym for ddirmn
|
dgdirm |
now a synonym for dgdirmn
|
dneg |
now a synonym for dnegmn
|
rdirm |
now a synonym for rdirmn
|
rgdirm |
now a synonym for rgdirmn
|
The function dneg
has been deprecated. Use dnegmn
instead.
Note the change in argument order:
dneg(Y, prob, beta)
and dnegmn(Y, alpha, beta)
from MGLM_0.0.8 have been deprecated;
use dnegmn(Y, beta, prob = alpha/(rowSums(alpha)+1), alpha=NULL)
instead.
"MGLMfit"
A class containing the model fitting results from the MGLMfit
.
estimate
object of class "vector"
, containing the parameter estimates.
SE
object of class "vector"
,
containing the standard errors of the estimates.
vcov
object of class "matrix"
,
the variance covariance matrix of the parameter estimates.
logL
object of class "numeric"
,
the fitted log likelihood.
BIC
object of class "numeric"
,
Bayesian information criterion.
AIC
object of class "numeric"
,
Akaike information criterion.
LRTpvalue
object of class "numeric"
,
likelihood ratio test p value.
gradient
object of class "numeric"
or "matrix"
,
containing the gradient.
iter
object of class "numeric"
,
number of iteration used.
distribution
object of class "character"
,
the distribution fitted.
fitted
object of class "vector"
,
the fitted mean of each category.
LRT
object of class "numeric"
,
the likelihood ratio test statistic.
Yiwen Zhang and Hua Zhou
showClass("MGLMfit")
showClass("MGLMfit")
MGLMreg
fits multivariate response generalized linear models, specified by a symbolic description of the linear predictor and a description of the error distribution.
MGLMreg( formula, data, dist, init = NULL, weight = NULL, epsilon = 1e-08, maxiters = 150, display = FALSE, LRT = FALSE, parallel = FALSE, cores = NULL, cl = NULL, sys = NULL, regBeta = FALSE ) MGLMreg.fit( Y, init = NULL, X, dist, weight = NULL, epsilon = 1e-08, maxiters = 150, display = FALSE, LRT = FALSE, parallel = FALSE, cores = NULL, cl = NULL, sys = NULL, regBeta = FALSE )
MGLMreg( formula, data, dist, init = NULL, weight = NULL, epsilon = 1e-08, maxiters = 150, display = FALSE, LRT = FALSE, parallel = FALSE, cores = NULL, cl = NULL, sys = NULL, regBeta = FALSE ) MGLMreg.fit( Y, init = NULL, X, dist, weight = NULL, epsilon = 1e-08, maxiters = 150, display = FALSE, LRT = FALSE, parallel = FALSE, cores = NULL, cl = NULL, sys = NULL, regBeta = FALSE )
formula |
an object of class |
data |
an optional data frame, list or environment (or object coercible by as.data.frame to a data frame) containing the variables in the model. If not found in |
dist |
a description of the error distribution to fit. See |
init |
an optional matrix of initial value of the parameter estimates. Should have the compatible dimension with |
weight |
an optional vector of weights assigned to each row of the data. Should be |
epsilon |
an optional numeric controlling the stopping criterion. The algorithm terminates when the relative change in the loglikelihoods of two successive iterates is less than |
maxiters |
an optional numeric controlling the maximum number of iterations. The default value is |
display |
an optional logical variable controlling the display of iterations. The default value is |
LRT |
an optional logical variable controlling whether to perform likelihood ratio test on each predictor. The default value is |
parallel |
an optional logical variable controlling whether to perform parallel computing. On a multi-core Windows machine, a cluster is created based on socket; on a multi-core Linux/Mac machine, a cluster is created based on forking. The default value is |
cores |
an optional value specifying the number of cores to use. Default value is half of the logical cores. |
cl |
a cluster object, created by the package parallel or by package snow. If |
sys |
the operating system. Will be used when choosing parallel type. |
regBeta |
an optional logical variable. When |
Y , X
|
for |
The formula should be in the form responses ~ covariates where the responses are the multivariate count matrix or a few columns from a data frame which is specified by data
. The covariates are either matrices or from the data frame. The covariates can be numeric or character or factor.
See dist
for details about distributions.
Instead of using the formula, the user can directly input the design matrix and the response vector using MGLMreg.fit
function.
Returns an object of class "MGLMreg"
. An object of class "MGLMreg"
is a list containing the following components:
coefficients
the estimated regression coefficients.
SE
the standard errors of the estimates.
Hessian
the Hessian at the estimated parameter values.
gradient
the gradient at the estimated parameter values.
wald.value
the Wald statistics.
wald.p
the p values of Wald test.
test
test statistic and the corresponding p-value. If LRT=FALSE
, only returns test resultsfrom Wald test; if LRT=TRUE
, returns the test results from both Wald test and likelihood ratio test.
logL
the final loglikelihood.
BIC
Bayesian information criterion.
AIC
Akaike information criterion.
fitted
the fitted values from the regression model
iter
the number of iterations used.
call
the matched call.
distribution
the distribution fitted.
data
the data used to fit the model.
Dof
degrees of freedom.
Yiwen Zhang and Hua Zhou
See also MGLMfit
for distribution fitting.
##----------------------------------------## ## Generate data n <- 2000 p <- 5 d <- 4 m <- rep(20, n) set.seed(1234) X <- 0.1* matrix(rnorm(n*p),n, p) alpha <- matrix(1, p, d-1) beta <- matrix(1, p, d-1) Alpha <- exp(X %*% alpha) Beta <- exp(X %*% beta) gdm.Y <- rgdirmn(n, m, Alpha, Beta) ##----------------------------------------## ## Regression gdm.reg <- MGLMreg(gdm.Y~X, dist="GDM", LRT=FALSE)
##----------------------------------------## ## Generate data n <- 2000 p <- 5 d <- 4 m <- rep(20, n) set.seed(1234) X <- 0.1* matrix(rnorm(n*p),n, p) alpha <- matrix(1, p, d-1) beta <- matrix(1, p, d-1) Alpha <- exp(X %*% alpha) Beta <- exp(X %*% beta) gdm.Y <- rgdirmn(n, m, Alpha, Beta) ##----------------------------------------## ## Regression gdm.reg <- MGLMreg(gdm.Y~X, dist="GDM", LRT=FALSE)
"MGLMreg"
Objects can be created by calls of the form new("MGLMreg", ...)
.
call
object of class "call"
.
data
object of class "list"
,
consists of both the predictor matrix and the response matrix.
coefficients
object of class "list"
or "matrix"
,
the estimated parameters.
SE
object of class "list"
or "matrix"
,
the standard errors of the parameters.
test
object of class "matrix"
,
the test statistics and p-values.
Hessian
object of class "matrix"
,
the Hessian matrix.
logL
object of class "numeric"
,
the loglikelihood.
BIC
object of class "numeric"
,
AIC
object of class "numeric"
,
Akaike information criterion.
iter
object of class "numeric"
,
the number of iteration used.
distribution
object of class "character"
,
the distribution fitted.
fitted
object of class "vector"
,
the fitted value.
gradient
object of class "numeric"
or "matrix"
,
the gradient at the estimated parameter values.
wald.value
object of class "numeric"
or "logical"
,
the Wald statistics.
wald.p
object of class "numeric"
or "logical"
,
the p values of Wald test.
Dof
object of class "numeric"
,
the degrees of freedom.
Yiwen Zhang and Hua Zhou
showClass("MGLMreg")
showClass("MGLMreg")
Fit sparse regression in multivariate generalized linear models.
MGLMsparsereg( formula, data, dist, lambda, penalty, weight, init, penidx, maxiters = 150, ridgedelta, epsilon = 1e-05, regBeta = FALSE, overdisp ) MGLMsparsereg.fit( Y, X, dist, lambda, penalty, weight, init, penidx, maxiters = 150, ridgedelta, epsilon = 1e-05, regBeta = FALSE, overdisp )
MGLMsparsereg( formula, data, dist, lambda, penalty, weight, init, penidx, maxiters = 150, ridgedelta, epsilon = 1e-05, regBeta = FALSE, overdisp ) MGLMsparsereg.fit( Y, X, dist, lambda, penalty, weight, init, penidx, maxiters = 150, ridgedelta, epsilon = 1e-05, regBeta = FALSE, overdisp )
formula |
an object of class |
data |
an optional data frame, list or environment (or object coercible by
|
dist |
a description of the error distribution to fit. See |
lambda |
penalty parameter. |
penalty |
penalty type for the regularization term. Can be chosen from |
weight |
an optional vector of weights assigned to each row of the data.
Should be |
init |
an optional matrix of initial value of the parameter estimates.
Should have the compatible dimension with the data. See |
penidx |
a logical vector indicating the variables to be penalized. The default value is |
maxiters |
an optional numeric controlling the maximum number of iterations. The default value is maxiters=150. |
ridgedelta |
an optional numeric controlling the behavior of the Nesterov's accelerated proximal gradient method. The default value is |
epsilon |
an optional numeric controlling the stopping criterion. The algorithm terminates when the relative change in the objective values of two successive iterates is less then |
regBeta |
an optional logical variable used when running negative multinomial regression ( |
overdisp |
an optional numerical variable used only when fitting sparse negative multinomial
model |
Y |
a matrix containing the multivariate categorical response data.
Rows of the matrix represent observations, while columns are the different
categories. Rows and columns of all zeros are automatically removed when
|
X |
design matrix (including intercept).
Number of rows of the matrix should match that of |
In general, we consider regularization problem
where is the loglikelihood function and
is the
regularization function.
Sparsity in the individual elements of the parameter matrix is achieved
by the lasso penalty (
dist="sweep"
)
Sparsity in the rows of the regression parameter matrix is achieved
by the group penalty (
dist="group"
)
where is the
norm of a vector
. In other words,
is the
norm of the
-th row of the
parameter matrix
.
Sparsity in the rank of the parameter matrix is achieved by the nuclear norm penalty (
dist="nuclear"
)
where are the singular values of the parameter matrix
.
The nuclear norm
is a convex relaxation of
.
See dist
for details about distributions.
Returns an object of class "MGLMsparsereg"
. An object of class "MGLMsparsereg"
is a list containing at least the following components:
coefficients
the estimated matrix of regression coefficients.
logL
the final loglikelihood value.
AIC
Akaike information criterion.
BIC
Bayesian information criterion.
Dof
degrees of freedom of the estimated parameter.
iter
number of iterations used.
maxlambda
the maxmum tuning parameter such that
the estimated coefficients are not all zero. This value is returned only
when the tuning parameter lambda
given to the function is large enough
such that all the parameter estimates are zero; otherwise, maxlambda
is not computed.
call
a matched call.
data
the data used to fit the model: a list of the predictor matrix
and the response matrix.
penalty
the penalty chosen when running the penalized regression.
Yiwen Zhang and Hua Zhou
## Generate Dirichlet Multinomial data dist <- "DM" n <- 100 p <- 15 d <- 5 m <- runif(n, min=0, max=25) + 25 set.seed(134) X <- matrix(rnorm(n*p),n, p) alpha <- matrix(0, p, d) alpha[c(1,3, 5), ] <- 1 Alpha <- exp(X%*%alpha) Y <- rdirmn(size=m, alpha=Alpha) ## Tuning ngridpt <- 10 p <- ncol(X) d <- ncol(Y) pen <- 'nuclear' spfit <- MGLMsparsereg(formula=Y~0+X, dist=dist, lambda=Inf, penalty=pen)
## Generate Dirichlet Multinomial data dist <- "DM" n <- 100 p <- 15 d <- 5 m <- runif(n, min=0, max=25) + 25 set.seed(134) X <- matrix(rnorm(n*p),n, p) alpha <- matrix(0, p, d) alpha[c(1,3, 5), ] <- 1 Alpha <- exp(X%*%alpha) Y <- rdirmn(size=m, alpha=Alpha) ## Tuning ngridpt <- 10 p <- ncol(X) d <- ncol(Y) pen <- 'nuclear' spfit <- MGLMsparsereg(formula=Y~0+X, dist=dist, lambda=Inf, penalty=pen)
"MGLMsparsereg"
A class containing the results from the MGLMsparsereg
.
call
object of class "call"
.
data
object of class "list"
,
consists of both the predictor matrix and the response matrix.
coefficients
object of class "matrix"
,
the estimated parameters.
logL
object of class "numeric"
,
the loglikelihood.
BIC
object of class "numeric"
,
AIC
object of class "numeric"
,
Akaike information criterion.
Dof
object of class "numeric"
,
the degrees of freedom.
iter
object of class "numeric"
,
the number of iteration used.
maxlambda
object of class "numeric"
,
the maximum tuning parameter that ensures the estimated regression coefficients are not all zero.
lambda
object of class "numeric"
,
the tuning parameter used.
distribution
object of class "character"
,
the distribution fitted.
penalty
Object of class "character"
,
the chosen penalty when running penalized regression.
Yiwen Zhang and Hua Zhou
showClass("MGLMsparsereg")
showClass("MGLMsparsereg")
Finds the tuning parameter value that yields the smallest BIC.
MGLMtune( formula, data, dist, penalty, lambdas, ngridpt, warm.start = TRUE, keep.path = FALSE, display = FALSE, init, weight, penidx, ridgedelta, maxiters = 150, epsilon = 1e-05, regBeta = FALSE, overdisp )
MGLMtune( formula, data, dist, penalty, lambdas, ngridpt, warm.start = TRUE, keep.path = FALSE, display = FALSE, init, weight, penidx, ridgedelta, maxiters = 150, epsilon = 1e-05, regBeta = FALSE, overdisp )
formula |
an object of class |
data |
an optional data frame, list or environment (or object coercible by |
dist |
a description of the distribution to fit. See |
penalty |
penalty type for the regularization term. Can be chosen from |
lambdas |
an optional vector of the penalty values to tune. If missing, the vector of penalty values will be set inside the function. |
ngridpt |
an optional numeric variable specifying the number of grid points to tune. If |
warm.start |
an optional logical variable to specify whether to give warm start at each tuning grid point. If |
keep.path |
an optional logical variable controling whether to output the whole solution path. The default is |
display |
an optional logical variable to specify whether to show each tuning step. |
init |
an optional matrix of initial value of the parameter estimates. Should have the compatible dimension with the data. See |
weight |
an optional vector of weights assigned to each row of the data. Should be |
penidx |
a logical vector indicating the variables to be penalized. The default value is |
ridgedelta |
an optional numeric controlling the behavior of the Nesterov's accelerated proximal gradient method. The default value is |
maxiters |
an optional numeric controlling the maximum number of iterations. The default value is |
epsilon |
an optional numeric controlling the stopping criterion. The algorithm terminates when the relative change in the objective values of two successive iterates is less then |
regBeta |
an optional logical variable used when running negative multinomial regression ( |
overdisp |
an optional numerical variable used only when fitting sparse negative multinomial model and |
select
the final sparse regression result, using the optimal tuning parameter.
path
a data frame with degrees of freedom and BICs at each lambda.
Yiwen Zhang and Hua Zhou
set.seed(118) n <- 50 p <- 10 d <- 5 m <- rbinom(n, 100, 0.8) X <- matrix(rnorm(n * p), n, p) alpha <- matrix(0, p, d) alpha[c(1, 3, 5), ] <- 1 Alpha <- exp(X %*% alpha) Y <- rdirmn(size=m, alpha=Alpha) sweep <- MGLMtune(Y ~ 0 + X, dist="DM", penalty="sweep", ngridpt=10) show(sweep)
set.seed(118) n <- 50 p <- 10 d <- 5 m <- rbinom(n, 100, 0.8) X <- matrix(rnorm(n * p), n, p) alpha <- matrix(0, p, d) alpha[c(1, 3, 5), ] <- 1 Alpha <- exp(X %*% alpha) Y <- rdirmn(size=m, alpha=Alpha) sweep <- MGLMtune(Y ~ 0 + X, dist="DM", penalty="sweep", ngridpt=10) show(sweep)
"MGLMtune"
A class containing the results from the MGLMtune
.
call
object of class "call"
.
select
object of class "MGLMsparsereg"
,
regularized regression results given by the optimal tuning parameter.
path
object of class "data.frame"
,
the BIC, AIC, log-likelihood and degrees of freedom given each tuning parameter.
select.list
object of class "list"
,
the regularized regression results at each tuning grid point.
Yiwen Zhang and Hua Zhou
showClass("MGLMtune")
showClass("MGLMtune")
path
extracts from object of the class MGLMtune
the path of
BIC, AIC, log-likelihood and degrees of freedom given each tuning parameter.
## S4 method for signature 'MGLMtune' path(object)
## S4 method for signature 'MGLMtune' path(object)
object |
an object of class |
Returns a path of object
.
library("MGLM") dist <- "DM" n <- 100 p <- 10 d <- 5 set.seed(118) m <- rbinom(n, 200, 0.8) X <- matrix(rnorm(n * p), n, p) alpha <- matrix(0, p, d) alpha[c(1, 3, 5), ] <- 1 Alpha <- exp(X %*% alpha) Y <- rdirmn(size = m, alpha = Alpha) select <- MGLMtune(Y ~ 0 + X, dist = "DM", penalty = "nuclear", ngridpt = 10, display = FALSE) select_path <- path(select)
library("MGLM") dist <- "DM" n <- 100 p <- 10 d <- 5 set.seed(118) m <- rbinom(n, 200, 0.8) X <- matrix(rnorm(n * p), n, p) alpha <- matrix(0, p, d) alpha[c(1, 3, 5), ] <- 1 Alpha <- exp(X %*% alpha) Y <- rdirmn(size = m, alpha = Alpha) select <- MGLMtune(Y ~ 0 + X, dist = "DM", penalty = "nuclear", ngridpt = 10, display = FALSE) select_path <- path(select)
Predict using the fitted model from MGLMreg
when given a new set of covariates.
## S4 method for signature 'MGLMreg' predict(object, newdata)
## S4 method for signature 'MGLMreg' predict(object, newdata)
object |
model object. |
newdata |
new covariates data matrix. |
Outputs the probabilities of each category.
This helps answer questions such as whether certain features increase the probability of observing category j.
n <- 200 p <- 5 d <- 4 X <- matrix(runif(p * n), n, p) alpha <- matrix(c(0.6, 0.8, 1), p, d - 1, byrow=TRUE) alpha[c(1, 2),] <- 0 Alpha <- exp(X %*% alpha) beta <- matrix(c(1.2, 1, 0.6), p, d - 1, byrow=TRUE) beta[c(1, 2),] <- 0 Beta <- exp(X %*% beta) m <- runif(n, min=0, max=25) + 25 Y <- rgdirmn(n, m, Alpha, Beta) gdmReg <- MGLMreg(Y~0+X, dist="GDM") newX <- matrix(runif(1*p), 1, p) pred <- predict(gdmReg, newX)
n <- 200 p <- 5 d <- 4 X <- matrix(runif(p * n), n, p) alpha <- matrix(c(0.6, 0.8, 1), p, d - 1, byrow=TRUE) alpha[c(1, 2),] <- 0 Alpha <- exp(X %*% alpha) beta <- matrix(c(1.2, 1, 0.6), p, d - 1, byrow=TRUE) beta[c(1, 2),] <- 0 Beta <- exp(X %*% beta) m <- runif(n, min=0, max=25) + 25 Y <- rgdirmn(n, m, Alpha, Beta) gdmReg <- MGLMreg(Y~0+X, dist="GDM") newX <- matrix(runif(1*p), 1, p) pred <- predict(gdmReg, newX)
ddirmn
computes the log of the Dirichlet multinomial probability mass function.
rdirmn
generates Dirichlet multinomially distributed random number vectors.
rdirmn(n, size, alpha) ddirmn(Y, alpha)
rdirmn(n, size, alpha) ddirmn(Y, alpha)
n |
number of random vectors to generate. When |
size |
a number or vector specifying the total number of objects that are put into d categories in the Dirichlet multinomial distribution. |
alpha |
the parameter of the Dirichlet multinomial distribution. Can be a numerical positive vector or matrix.
For For |
Y |
The multivariate count matrix with dimensions |
When the multivariate count data exhibits over-dispersion, the traditional
multinomial model is insufficient. Dirichlet multinomial distribution models the
probabilities of the categories by a Dirichlet distribution.
Given the parameter vector ,
the probability mass of
-category count vector
,
under Dirichlet multinomial distribution is
where . Here,
, often read as "
choose
",
refers the number of
combinations from a set of
elements.
The parameter can be a vector of length
,
such as the results from the distribution fitting.
can also be a matrix with
rows, such as the inverse link
calculated from the regression parameter estimate
.
For each count vector and each corresponding parameter vector
, the function
ddirmn
returns the value .
When
Y
is a matrix of rows,
ddirmn
returns a vector of length .
rdirmn
returns a matrix of the generated random observations.
m <- 20 alpha <- c(0.1, 0.2) dm.Y <- rdirmn(n=10, m, alpha) pdfln <- ddirmn(dm.Y, alpha)
m <- 20 alpha <- c(0.1, 0.2) dm.Y <- rdirmn(n=10, m, alpha) pdfln <- ddirmn(dm.Y, alpha)
rgdirmn
generates random observations from the generalized Dirichlet multinomial distribution.
dgdirmn
computes the log of the generalized Dirichlet multinomial probability mass function.
rgdirmn(n, size, alpha, beta) dgdirmn(Y, alpha, beta)
rgdirmn(n, size, alpha, beta) dgdirmn(Y, alpha, beta)
n |
the number of random vectors to generate. When |
size |
a number or vector specifying the total number of objects that are put into d categories in the generalized Dirichlet multinomial distribution. |
alpha |
the parameter of the generalized Dirichlet multinomial distribution.
For For |
beta |
the parameter of the generalized Dirichlet multinomial distribution. For |
Y |
the multivariate count matrix with dimensions |
are the
category count vectors. Given the parameter vector
, and
,
the generalized Dirichlet multinomial probability mass function is
where and
.
Here,
, often read as "
choose
",
refers the number of
combinations from a set of
elements.
The and
parameters can be vectors, like the results from the
distribution
fitting function, or they can be matrices with
rows,
like the estimate
from the regression function multiplied by the covariate matrix
and
dgdirmn
returns the value of
.
When
Y
is a matrix of rows, the function
dgdirmn
returns a vector of length .
rgdirmn
returns a matrix of the generated random observations.
# example 1 m <- 20 alpha <- c(0.2, 0.5) beta <- c(0.7, 0.4) Y <- rgdirmn(10, m, alpha, beta) dgdirmn(Y, alpha, beta) # example 2 set.seed(100) alpha <- matrix(abs(rnorm(40)), 10, 4) beta <- matrix(abs(rnorm(40)), 10, 4) size <- rbinom(10, 10, 0.5) GDM.rdm <- rgdirmn(size=size, alpha=alpha, beta=beta) GDM.rdm1 <- rgdirmn(n=20, size=10, alpha=abs(rnorm(4)), beta=abs(rnorm(4)))
# example 1 m <- 20 alpha <- c(0.2, 0.5) beta <- c(0.7, 0.4) Y <- rgdirmn(10, m, alpha, beta) dgdirmn(Y, alpha, beta) # example 2 set.seed(100) alpha <- matrix(abs(rnorm(40)), 10, 4) beta <- matrix(abs(rnorm(40)), 10, 4) size <- rbinom(10, 10, 0.5) GDM.rdm <- rgdirmn(size=size, alpha=alpha, beta=beta) GDM.rdm1 <- rgdirmn(n=20, size=10, alpha=abs(rnorm(4)), beta=abs(rnorm(4)))
rmn
generates random number vectors given alpha
.
The function rmn(n, size, alpha)
calls rmultinom(n, size, prob)
after converting alpha
to probability.
dmn
computes the log of multinomial probability mass function.
rmn(n, size, alpha) dmn(Y, prob)
rmn(n, size, alpha) dmn(Y, prob)
n |
number of random vectors to generate. |
size |
a scalar or a vector. |
alpha |
a vector or a matrix. |
Y |
the multivariate count matrix with dimension |
prob |
the probability parameter of the multinomial distribution. |
A multinomial distribution models the counts of possible outcomes.
The counts of categories are negatively correlated.
is a
category count vector.
Given the parameter vector
,
,
, the function calculates the log of the multinomial pmf
where . Here,
, often read as "
choose
",
refers the number of
combinations from a set of
elements.
The parameter can be one vector, like the result from the distribution
fitting function; or,
can be a matrix with
rows, like the estimate
from the regression function,
where
and
.
The
-th column of the coefficient matrix
is set to
to avoid the identifiability issue.
The function dmn
returns the value of .
When
Y
is a matrix of rows, the function returns a
vector of length
.
The function rmn
returns multinomially distributed random number vectors
Yiwen Zhang and Hua Zhou
m <- 20 prob <- c(0.1, 0.2) dm.Y <- rdirmn(n=10, m, prob) pdfln <- dmn(dm.Y, prob)
m <- 20 prob <- c(0.1, 0.2) dm.Y <- rdirmn(n=10, m, prob) pdfln <- dmn(dm.Y, prob)
RNA-seq data simulated following the standard procedures (provided by Dr. Wei Sun, [email protected]).
rnaseq
rnaseq
A data frame containing 10 columns and 100 rows. The first 6 columns are the expression counts of 6 exons of a gene; the last four columns are the covariates: age, gender, treatment, and total number of reads.
Dr. Sun Wei, [email protected]
dnegmn
calculates the log of the negative multinomial probability mass function.
rnegmn
generates random observations from the negative multinomial distribution.
rnegmn(n, beta, prob) dnegmn(Y, beta, prob = alpha/(rowSums(alpha) + 1), alpha = NULL)
rnegmn(n, beta, prob) dnegmn(Y, beta, prob = alpha/(rowSums(alpha) + 1), alpha = NULL)
n |
number of random vectors to generate. When |
beta |
the over dispersion parameter of the negative multinomial distribution. |
prob |
the probability parameter of the negative multinomial distribution. Should be a numerical non-negative vector or matrix. For For |
Y |
the multivariate response matrix of dimension |
alpha |
an alternative way to specify the probability. Default value is |
is a
category vector. Given the parameter vector
,
,
and
,
, the negative multinomial probability mass function is
where . Here,
, often read as "
choose
",
refers the number of
combinations from a set of
elements.
alpha
is an alternative way to specify the probability:
for and
.
The parameter prob
can be a vector and beta
is a scalar; prob
can also
be a matrix with rows, and
beta
is a vector of length
like the estimate from the regression function
multiplied by the covariate matrix.
dnegmn
returns the value of
. When
Y
is a matrix of rows, the function
returns a vector of length
.
rnegmn
returns a matrix of the generated random observations.
Yiwen Zhang and Hua Zhou
###-----------------------### set.seed(128) n <- 100 d <- 4 p <- 5 a <- -matrix(1,p,d) X <- matrix(runif(n*p), n, p ) alpha <- exp(X%*%a) prob <- alpha/(rowSums(alpha)+1) beta <- exp(X%*%matrix(1,p)) Y <- rnegmn(n, beta, prob) ###-----------------------### m <- 20 n <- 10 p <- 5 d <- 6 a <- -matrix(1,p,d) X <- matrix(runif(n*p), n, p ) alpha <- exp(X%*%a) prob <- alpha/(rowSums(alpha)+1) b <- exp(X%*%rep(0.3,p)) Y <- rnegmn(prob=prob, beta=rep(10, n)) dnegmn(Y, b, prob)
###-----------------------### set.seed(128) n <- 100 d <- 4 p <- 5 a <- -matrix(1,p,d) X <- matrix(runif(n*p), n, p ) alpha <- exp(X%*%a) prob <- alpha/(rowSums(alpha)+1) beta <- exp(X%*%matrix(1,p)) Y <- rnegmn(n, beta, prob) ###-----------------------### m <- 20 n <- 10 p <- 5 d <- 6 a <- -matrix(1,p,d) X <- matrix(runif(n*p), n, p ) alpha <- exp(X%*%a) prob <- alpha/(rowSums(alpha)+1) b <- exp(X%*%rep(0.3,p)) Y <- rnegmn(prob=prob, beta=rep(10, n)) dnegmn(Y, b, prob)
Display the object by printing its class.
## S4 method for signature 'MGLMfit' show(object) ## S4 method for signature 'MGLMreg' show(object) ## S4 method for signature 'MGLMsparsereg' show(object) ## S4 method for signature 'MGLMtune' show(object)
## S4 method for signature 'MGLMfit' show(object) ## S4 method for signature 'MGLMreg' show(object) ## S4 method for signature 'MGLMsparsereg' show(object) ## S4 method for signature 'MGLMtune' show(object)
object |
an object to be printed. Should be of class |
library("MGLM") data("rnaseq") data <- rnaseq[, 1:6] gdmFit <- MGLMfit(data, dist = "GDM") show(gdmFit)
library("MGLM") data("rnaseq") data <- rnaseq[, 1:6] gdmFit <- MGLMfit(data, dist = "GDM") show(gdmFit)