R/DMD.MN.R In MGLM: Multivariate Response Generalized Linear Models

Defines functions DMD.MN.fitDMD.MN.regDMD.MN.predictDMD.MN.loss

```##============================================================##
## Function: fit Multinomial
##============================================================##
DMD.MN.fit <- function(data, weight) {

N <- nrow(data)
d <- ncol(data)

t <- sum(data)
m <- rowSums(data)
if (missing(weight)) weight <- rep(1, N)
estimate <-apply(data/m, 2, mean)
SE <- sqrt(estimate * (1-estimate)/N)
p_MN <- t(matrix(estimate, d, N))
logL <- sum(weight * data * log(p_MN)) + sum(lgamma(m + 1)) - sum(lgamma(data + 1))
BIC <- -2 * logL + log(N) * (d - 1)
AIC <- -2 * logL + 2 * (d-1)

##----------------------------------------##
## Clean up the results
##----------------------------------------##
estimate = apply(data/m, 2, mean)
if (!is.null(colnames(data))) {
names(estimate) <- paste("alpha", colnames(data), sep = "_")
} else {
names(estimate) <- paste("alpha", 1:d, sep = "_")
}

##----------------------------------------##
list(estimate = estimate, SE = SE, logL = logL, BIC = BIC,
AIC = AIC, distribution = "Multinomial", LRT = numeric(),
LRTpvalue = numeric(), iter = numeric(),
vcov = matrix(), gradient = matrix(), fitted = vector())
}

## ============================================================##
## Function: fit Multinomial Reg The dimension of Y is n by d (d>=3)
## ============================================================##
DMD.MN.reg <- function(Y, init, X, weight, epsilon, maxiters, display, parallel,
cores, cl, sys) {

## ----------------------------------------##
## Function to calculate log-likelihood
## ----------------------------------------##
dmultn_L <- function(B1) {
B <- cbind(B1, 0)
alpha <- exp(X %*% B)
return(sum(Y * (X %*% B) * weight) - sum(Y * log(rowSums(alpha)) * weight) +
sum(lgamma(m + 1) * weight) - sum(lgamma(Y + 1) * weight))
}

##----------------------------------------##
##----------------------------------------##
dl.MN_fctn <- function(i, beta) {
x <- (weight * X)[i, ]
y <- Y[i, ]
prob <- exp(colSums(x * beta))/(sum(exp(colSums(x * beta))) + 1)
dl <- y[-ncol(Y)] - sum(y) * prob
return(matrix(c(outer(x, dl)), (length(y) - 1) * length(x), 1))
}

## ----------------------------------------##
## Function to calculate hessian matrix
## ----------------------------------------##
H.MN_fctn <- function(i, beta) {
x <- (weight * X)[i, ]
y <- Y[i, ]
prob <- exp(colSums(x * beta))/(sum(exp(colSums(x * beta))) + 1)
if (d > 2)
dprob <- diag(prob) else if (d == 2)
dprob <- prob
d2l <- sum(y) * (outer(prob, prob) - dprob)
o.x <- outer(x, x)
H <- d2l %x% o.x
return(H)
}

## ----------------------------------------##
## Keep some original values
## ----------------------------------------##
pred <- matrix(NA, nrow(Y), ncol(Y))
emptyRow <- rowSums(Y) == 0
Xf <- X
ow <- getOption("warn")

weight <- weight[rowSums(Y) != 0]
X <- as.matrix(X[rowSums(Y) != 0, ])
Y <- as.matrix(Y[rowSums(Y) != 0, colSums(Y) != 0])
d <- ncol(Y)
m <- rowSums(Y)
p <- ncol(X)
N <- nrow(Y)
beta <- init
lliter <- rep(0, maxiters)
lliter[1] <- dmultn_L(init)
niter <- 1
options(warn = -1)

## ----------------------------------------##
## Begin the main loop
## ----------------------------------------##
while (((niter <= 2) || ((ll2 - ll1)/(abs(ll1) + 1) > epsilon)) & (niter < maxiters)) {

niter <- niter + 1
ll1 <- lliter[niter - 1]

## ----------------------------------------##
## Newton update
## ----------------------------------------##
P <- exp(X %*% beta)
P <- P/(rowSums(P) + 1)
dl <- Y[, -d] - P * m
if (d > 2) {
dl <- colSums(kr(dl, X, weight))
} else if (d == 2) {
dl <- colSums(dl * weight * X)
}
# dl <- Reduce('+', lapply(1:N, dl.MN_fctn, beta) )
H <- Reduce("+", lapply(1:N, H.MN_fctn, beta))
update <- NULL
try(update <- solve(H, dl), silent = TRUE)
if (is.null(update)) {
ll.Newton <- NA
} else if (is.numeric(update)) {
beta_Newton <- beta - matrix(update, p, d - 1)
ll.Newton <- dmultn_L(beta_Newton)
## ----------------------------------------##
## Half stepping
## ----------------------------------------##
if (is.nan(ll.Newton) || ll.Newton >= 0) {
ll.Newton <- NA
} else if (!is.na(ll.Newton) & ll1 >= ll.Newton) {
for (step in 1:20) {
beta_N <- beta - matrix(update/(2^step), p, d - 1)
llnew <- dmultn_L(beta_N)
if (is.nan(llnew) || llnew > 0) {
next
} else if (llnew > ll.Newton) {
ll.Newton <- llnew
beta_Newton <- beta_N
}
if (is.na(llnew) | llnew > ll1) {
break
}
}
}
if (is.nan(ll.Newton) | ll.Newton >= 0)
ll.Newton <- NA
} else {
ll.Newton <- NA
}

if (is.na(ll.Newton) || ll.Newton < ll1) {
## ----------------------------------------## MM update
denominator <- 1 + rowSums(exp(X %*% beta))
a <- m/denominator
beta_MM <- matrix(0, p, (d - 1))
if (!parallel) {
for (i in 1:(d - 1)) {
beta_MM[, i] <- glm.fit(X, Y[, i]/a, weights = weight * a, family = poisson(link = "log"))\$coefficients
}
} else {
y.list <- split(Y/a, rep(1:d, each = nrow(Y)))
y.list[[d]] <- NULL
if (sys == "Windows") {
fit.list <- clusterMap(cl, glm.private, y.list, .scheduling = "dynamic",
MoreArgs = list(weights = weight * a, X = X, family = poisson(link = "log")))
} else if (sys != "Windows") {
fit.list <- mcmapply(cl, glm.private, y.list, MoreArgs = list(weights = weight *
a, X = X, family = poisson(link = "log")), mc.cores = cores,
mc.preschedule = FALSE)
}
beta_MM <- do.call("cbind", fit.list)
}
ll.MM <- dmultn_L(beta_MM)
## ----------------------------------------## Choose the update
if (is.na(ll.Newton) | (ll.MM < 0 & ll.MM > ll1)) {
if (display)
print(paste("Iteration ", niter, " MM update", ll.MM, sep = ""))
beta <- beta_MM
lliter[niter] <- ll.MM
ll2 <- ll.MM
}
} else {
if (display)
print(paste("Iteration ", niter, " Newton's update", ll.Newton, sep = ""))
beta <- beta_Newton
lliter[niter] <- ll.Newton
ll2 <- ll.Newton
}
}
## ----------------------------------------##
## End of the main loop
## ----------------------------------------##
options(warn = ow)

## ----------------------------------------##
## Compute output statistics
## ----------------------------------------##
BIC <- -2 * lliter[niter] + log(N) * p * (d - 1)
AIC <- 2 * p * (d - 1) - 2 * lliter[niter]
P <- exp(X %*% beta)
P <- cbind(P, 1)
P <- P/rowSums(P)
if (d > 2) {
H <- kr(P[, 1:(d - 1)], X)
} else if (d == 2) {
H <- P[, 1] * X
}
H <- t(H) %*% (weight * m * H)
for (i in 1:(d - 1)) {
id <- (i - 1) * p + (1:p)
H[id, id] = H[id, id] - t(X) %*% (weight * m * P[, i] * X)
}

SE <- matrix(NA, p, (d - 1))
wald <- rep(NA, p)
wald.p <- rep(NA, p)

## ----------------------------------------##
## Check the estimate
## ----------------------------------------##
dl <- Reduce("+", lapply(1:N, dl.MN_fctn, beta))
if (mean(dl^2) > 0.001) {
"\n", sep = " "))
}
eig <- eigen(H)\$values
if (any(is.complex(eig)) || any(eig > 0)) {
warning("The estimate is a saddle point.")
} else if (any(eig == 0)) {
warning(paste("The Hessian matrix is almost singular.", "\n", sep = ""))
} else {
## ----------------------------------------##
## If Hessian is negative definite, estimate SE and Wald
## ----------------------------------------##
Hinv <- chol2inv(chol(-H))
SE <- matrix(sqrt(diag(Hinv)), p, (d - 1))
wald <- rep(0, p)
wald.p <- rep(0, p)
for (j in 1:p) {
id <- c(0:(d - 2)) * p + j
wald[j] <- beta[j, ] %*% solve(Hinv[id, id]) %*% beta[j, ]
wald.p[j] <- pchisq(wald[j], (d - 1), lower.tail = FALSE)
}
}

## ----------------------------------------##
## Fitted value
## ----------------------------------------##
tmppred <- exp(Xf %*% beta)
pred <- tmppred / (rowSums(tmppred) + 1)
pred <- cbind(pred, (1 - rowSums(pred)))
pred[emptyRow, ] <- 0

## ----------------------------------------##
## Clean up the results
## ----------------------------------------##
colnames(beta) <- colnames(Y)[1:(ncol(Y) - 1)]
colnames(SE) <- colnames(Y)[1:(ncol(Y) - 1)]

## ----------------------------------------##
list(coefficients = beta, SE = SE, Hessian = H, wald.value = wald, wald.p = wald.p,
DoF = p * (d - 1), logL = lliter[niter], BIC = BIC, AIC = AIC, iter = niter,
gradient = dl, fitted = pred, cl = cl, distribution = "Multinomial")

}

##============================================================##
## Function: predict
##============================================================##

DMD.MN.predict <- function(beta, d, newdata) {
pred <- exp(newdata %*% beta)/(rowSums(exp(newdata %*% beta)) + 1)
pred <- cbind(pred, (1 - rowSums(pred)))
return(pred)
}

##============================================================##
## Function: loss ; called by MGLMsparsereg
##============================================================##
DMD.MN.loss <- function(Y, X, beta, dist, weight){

N <- nrow(Y)
d <- ncol(Y)
p <- ncol(X)
m <- rowSums(Y)

P <- matrix(NA, N, d)
P[, d] <- rep(1, N)
P[, 1:(d - 1)] <- exp(X %*% beta)
P <- P/rowSums(P)
loss <- -sum(weight * dmn(Y, P))
kr1 <- Y[, -d] - P[, 1:(d - 1)] * m
if (d > 2) {
lossD1 <- -colSums(kr(kr1, X, weight))
} else if (d == 2) {
lossD1 <- -colSums(kr1 * weight * X)
}
lossD1 <- matrix(lossD1, p, (d - 1))

return(list(loss, lossD1))
}
```

Try the MGLM package in your browser

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

MGLM documentation built on May 2, 2019, 1:38 p.m.