Nothing
R/pois.envMU.R
In Renvlp: Computing Envelope Estimators
Defines functions
pois.envMU
pois.envMU <- function(X, Y, u, initial = NULL){
X <- as.matrix(X)
Y <- as.matrix(Y)
a <- dim(Y)
n <- a[1]
r <- a[2]
p <- ncol(X)
if (a[1] != nrow(X))
stop("X and Y should have the same number of observations.")
if (u > p || u < 0)
stop("u must be an integer between 0 and r.")
if (sum(duplicated(cbind(X, Y), MARGIN = 2)) > 0)
stop("Some responses also appear in the predictors, or there maybe duplicated columns in X or Y.")
options(warn = -1)
fit <- stats::glm(Y ~ X, family = stats::poisson)
mu <- as.vector(fit$coefficients[1])
beta <- as.vector(fit$coefficients[2 : (p + 1)])
theta <- mu + X %*% beta
c.theta <- - exp(theta)
c.theta.mean <- mean(c.theta)
weight <- c.theta / c.theta.mean
V <- theta + ((Y - c.theta) / weight)
wx <- sweep(X, 1, weight, '*')
mean.wx <- colMeans(wx)
wxx <- sweep(X, 2, mean.wx, '-')
sigwx <- crossprod(wxx, sweep(wxx, 1, weight, '*')) / n
wv <- weight * V
mean.wv <- mean(wv)
wvv <- as.matrix(V - mean.wv)
sigwxv <- crossprod(wxx, sweep(wvv, 1, weight, '*')) / n
inv.sigwx <- chol2inv(chol(sigwx))
M.init <- inv.sigwx / (- c.theta.mean)
U.init <- tcrossprod(beta)
tmp1 <- envMU(M.init, U.init, u, initial = initial)
gamma.init <- tmp1$Gammahat
gamma0.init <- tmp1$Gamma0hat
sigX <- stats::cov(X) * (n - 1) / n
invsigX <- chol2inv(chol(sigX))
if (u == 0) {
Gammahat <- NULL
Gamma0hat <- diag(p)
tmp.M <- eigen(sigX)
mu <- log(mean(Y))
muh <- matrix(1, n, 1) %*% mu
Cn1 <- t(Y) %*% muh - colSums(exp(muh))
eta <- NULL
var <- NULL
weight <- rep(1, n)
V <- Y
objfun <- Cn1 - n/2 * sum(log(tmp.M$values))
} else if (u == p) {
Gammahat <- diag(p)
Gamma0hat <- NULL
tmp.M <- eigen(sigX)
eta <- beta
var <- inv.sigwx / (- c.theta.mean)
Cn1 <- t(Y) %*% theta - colSums(exp(theta))
objfun <- Cn1 - n/2 * sum(log(tmp.M$values))
} else if (u == p - 1) {
maxiter = 1000
ftol = 0.001
M <- sigX
MU <- sigX
tmp.MU <- eigen(MU)
invMU <- sweep(tmp.MU$vectors, MARGIN = 2, 1/tmp.MU$values, "*") %*% t(tmp.MU$vectors)
Cn1 <- t(Y) %*% theta - colSums(exp(theta))
e1 <- eigen(crossprod(gamma.init, M) %*% gamma.init)
e2 <- eigen(crossprod(gamma.init, invMU) %*% gamma.init)
e3 <- eigen(sigX)
temp1 <- sum(log(e1$values)) + sum(log(e2$values)) + sum(log(e3$values))
obj1 <- Cn1 - (n / 2) * temp1
init <- gamma.init
GEidx <- GE(init)
Ginit <- init %*% solve(init[GEidx[1:u], ])
GiX <- X %*% Ginit
fit <- stats::glm(Y ~ GiX, family = stats::poisson)
eta <- as.matrix(fit$coefficients[2 : (u + 1)])
j <- GEidx[p]
g <- as.matrix(Ginit[j, ])
g1 <- Ginit[-j, ]
t2 <- crossprod(g1, as.matrix(M[-j, j]))/M[j, j]
t3 <- crossprod(g1, as.matrix(invMU[-j, j]))/invMU[j, j]
GUGt2 <- g + t2
GUG <- crossprod(Ginit, (M %*% Ginit)) - tcrossprod(GUGt2, GUGt2) * M[j, j]
GVGt2 <- g + t3
GVG <- crossprod(Ginit, (invMU %*% Ginit)) - tcrossprod(GVGt2, GVGt2) * invMU[j, j]
invC1 <- chol2inv(chol(GUG))
invC2 <- chol2inv(chol(GVG))
X1 <- X[ , j]
X2 <- X[ , -j]
t5 <- X1 %*% t(g) %*% eta
t6 <- matrix(1, n, 1) %*% mu + X2 %*% g1 %*% eta
t7 <- t6 + t5
t8 <- t(Y) %*% t7
et6 <- exp(t7)
t9 <- colSums(et6)
fobj <- function(x) {
tmp2 <- x + t2
tmp3 <- x + t3
Ginit[j, ] <- x
GiX <- X %*% Ginit
fit <- stats::glm(Y ~ GiX, family = stats::poisson)
eta <- as.matrix(fit$coefficients[2 : (u + 1)])
tmp4 <- X1 %*% t(x) %*% eta + t6
T2 <- invC1 %*% tmp2
T3 <- invC2 %*% tmp3
T4 <- exp(tmp4)
n /2 * (- 2 * log(1 + sum(x^2)) + log(1 + M[j, j] *
crossprod(tmp2,T2)) + log(1 + invMU[j, j] *
crossprod(tmp3,T3))) - t(Y) %*% tmp4 + colSums(T4)
}
gobj <- function(x) {
tmp2 <- x + t2
tmp3 <- x + t3
Ginit[j, ] <- x
GiX <- X %*% Ginit
fit <- stats::glm(Y ~ GiX, family = stats::poisson)
eta <- as.matrix(fit$coefficients[2 : (u + 1)])
tmp4 <- X1 %*% t(x) %*% eta + t6
Cn <- Y - exp(tmp4)
T2 <- invC1 %*% tmp2
T3 <- invC2 %*% tmp3
A1 <- crossprod(g1, sigwx[-j, -j]) %*% g1 + as.matrix(x) %*% sigwx[j, -j] %*% g1 + crossprod(g1, sigwx[-j, j]) %*% t(x) + tcrossprod(x) * sigwx[j , j]
invC3 <- chol2inv(chol(A1))
T5 <- crossprod(X, Cn) %*% t(eta)
tmp5 <- X1 %*% t(x) + X2 %*% g1
T6 <- tcrossprod(sigwxv, Cn) %*% tmp5 %*% invC3
tmp6 <- as.matrix(sigwx[ , -j]) %*% g1 + tcrossprod(sigwx[ , j], x)
T7 <- tmp6 %*% invC3 %*% crossprod(tmp5, Cn) %*% t(eta)
T8 <- tmp6 %*% eta %*% crossprod(Cn, tmp5) %*% invC3
r1 <- rep(0, p)
r1[j] <- 1
T9 <- t(r1) %*% (T5 + T6 - T7 - T8)
- t(T9) + n * (- 4 * x %*% solve(1 + sum(x^2)) + 2 * T2/as.numeric(1/M[j, j] +
crossprod(tmp2, T2)) + 2 * T3/as.numeric(1/invMU[j, j] + crossprod(tmp3, T3)))
}
i <- 1
while (i < maxiter) {
res <- stats::optim(Ginit[j, ], fobj, gobj, method = "BFGS")
Ginit[j, ] <- res$par
GX <- X %*% Ginit
fit <- stats::glm(Y ~ GX, family = stats::poisson)
mu <- as.vector(fit$coefficients[1])
eta <- as.matrix(fit$coefficients[2 : (u + 1)])
beta <- Ginit %*% eta
theta <- mu + X %*% beta
c.theta <- - exp(theta)
c.theta.mean <- mean(c.theta)
weight <- c.theta / c.theta.mean
V <- theta + ((Y - c.theta) / weight)
wx <- sweep(X, 1, weight, '*')
mean.wx <- colMeans(wx)
wxx <- sweep(X, 2, mean.wx, '-')
sigwx <- crossprod(wxx, sweep(wxx, 1, weight, '*')) / n
wv <- weight * V
mean.wv <- mean(wv)
wvv <- as.matrix(V - mean.wv)
sigwxv <- crossprod(wxx, sweep(wvv, 1, weight, '*')) / n
e1 <- eigen(t(Ginit) %*% M %*% Ginit)
e2 <- eigen(t(Ginit) %*% invMU %*% Ginit)
e3 <- eigen(crossprod(Ginit))
e4 <- mu + X %*% beta
e5 <- crossprod(Y, e4) - colSums(exp(e4))
obj2 <- - n/2 * (sum(log(e1$values)) + sum(log(e2$values))) + sum(log(e3$values)) + e5
if (abs(obj1 - obj2) < ftol) {
break
}
else {
obj1 <- obj2
i <- i + 1
}
}
a <- qr(Ginit)
Gammahat <- qr.Q(a)
GX <- X %*% Gammahat
fit <- stats::glm(Y ~ GX, family = stats::poisson)
mu <- as.vector(fit$coefficients[1])
eta <- matrix(fit$coefficients[2 : (u + 1)])
beta <- Gammahat %*% eta
theta <- mu + X %*% beta
c.theta <- - exp(theta)
c.theta.mean <- mean(c.theta)
weight <- c.theta / c.theta.mean
wx <- sweep(X, 1, weight, '*')
mean.wx <- colMeans(wx)
wxx <- sweep(X, 2, mean.wx, '-')
sigwx <- crossprod(wxx, sweep(wxx, 1, weight, '*')) / n
inv.sigwx <- chol2inv(chol(sigwx))
var <- inv.sigwx / (- c.theta.mean)
e1 <- eigen(t(Gammahat) %*% M %*% Gammahat)
e2 <- eigen(t(Gammahat) %*% invMU %*% Gammahat)
e3 <- eigen(M)
e4 <- mu + X %*% Gammahat %*% eta
e5 <- crossprod(Y, e4) - colSums(exp(e4))
Gamma0hat <- qr.Q(a, complete = TRUE)[, (u + 1):r, drop = FALSE]
objfun <- - n/2 * (sum(log(e1$values)) + sum(log(e2$values))) + sum(log(e3$values)) + e5
Gammahat <- as.matrix(Gammahat)
Gamma0hat <- as.matrix(Gamma0hat)
} else {
maxiter = 100
ftol = 0.001
M <- sigX
MU <- sigX
tmp.MU <- eigen(MU)
invMU <- sweep(tmp.MU$vectors, MARGIN = 2, 1/tmp.MU$values, "*") %*% t(tmp.MU$vectors)
Cn1 <- t(Y) %*% theta - colSums(exp(theta))
e1 <- eigen(crossprod(gamma.init, M) %*% gamma.init)
e2 <- eigen(crossprod(gamma.init, invMU) %*% gamma.init)
e3 <- eigen(sigX)
temp1 <- sum(log(e1$values)) + sum(log(e2$values)) + sum(log(e3$values))
obj1 <- Cn1 - (n / 2) * temp1
init <- gamma.init
GEidx <- GE(init)
Ginit <- init %*% solve(init[GEidx[1:u], ])
GiX <- X %*% Ginit
fit <- stats::glm(Y ~ GiX, family = stats::poisson)
eta <- as.matrix(fit$coefficients[2 : (u + 1)])
GUG <- crossprod(Ginit, (M %*% Ginit))
GVG <- crossprod(Ginit, (invMU %*% Ginit))
t4 <- crossprod(Ginit[GEidx[(u + 1):p], ], Ginit[GEidx[(u + 1):p], ]) + diag(u)
i <- 1
while (i < maxiter) {
for (j in GEidx[(u + 1):p]) {
g <- as.matrix(Ginit[j, ])
g1 <- as.matrix(Ginit[-j, ])
t2 <- crossprod(g1, as.matrix(M[-j, j]))/M[j, j]
t3 <- crossprod(g1, as.matrix(invMU[-j, j]))/invMU[j, j]
GUGt2 <- g + t2
GUG <- GUG - tcrossprod(GUGt2, GUGt2) * M[j, j]
GVGt2 <- g + t3
GVG <- GVG - tcrossprod(GVGt2, GVGt2) * invMU[j, j]
t4 <- t4 - tcrossprod(g, g)
invC1 <- chol2inv(chol(GUG))
invC2 <- chol2inv(chol(GVG))
invt4 <- chol2inv(chol(t4))
X1 <- X[ , j]
X2 <- X[ , -j]
t5 <- X1 %*% t(g) %*% eta
t6 <- mu + X2 %*% g1 %*% eta
t7 <- t5 + t6
t8 <- t(Y) %*% t7
et7 <- exp(t7)
t9 <- colSums(et7)
fobj <- function(x) {
tmp2 <- x + t2
tmp3 <- x + t3
Ginit[j, ] <- x
GiX <- X %*% Ginit
fit <- stats::glm(Y ~ GiX, family = stats::poisson)
eta <- as.matrix(fit$coefficients[2 : (u + 1)])
tmp4 <- X1 %*% t(x) %*% eta + t6
T1 <- invt4 %*% x
T2 <- invC1 %*% tmp2
T3 <- invC2 %*% tmp3
T4 <- exp(tmp4)
n /2 * (-2 * log(1 + x %*% T1) + log(1 + M[j, j] *
crossprod(tmp2, T2)) + log(1 + invMU[j, j] *
crossprod(tmp3, T3))) - t(Y) %*% tmp4 + colSums(T4)
}
gobj <- function(x) {
tmp2 <- x + t2
tmp3 <- x + t3
Ginit[j, ] <- x
GiX <- X %*% Ginit
fit <- stats::glm(Y ~ GiX, family = stats::poisson)
eta <- as.matrix(fit$coefficients[2 : (u + 1)])
tmp4 <- X1 %*% t(x) %*% eta + t6
Cn <- Y - exp(tmp4)
T1 <- invt4 %*% x
T2 <- invC1 %*% tmp2
T3 <- invC2 %*% tmp3
A1 <- crossprod(g1, sigwx[-j, -j]) %*% g1 + as.matrix(x) %*% sigwx[j, -j] %*% g1 + crossprod(g1, sigwx[-j, j]) %*% t(x) + tcrossprod(x) * sigwx[j , j]
invC3 <- chol2inv(chol(A1))
T5 <- crossprod(X, Cn) %*% t(eta)
tmp5 <- X1 %*% t(x) + X2 %*% g1
T6 <- tcrossprod(sigwxv, Cn) %*% tmp5 %*% invC3
tmp6 <- as.matrix(sigwx[ , -j]) %*% g1 + tcrossprod(sigwx[ , j], x)
T7 <- tmp6 %*% invC3 %*% crossprod(tmp5, Cn) %*% t(eta)
T8 <- tmp6 %*% eta %*% crossprod(Cn, tmp5) %*% invC3
r1 <- rep(0, p)
r1[j] <- 1
T9 <- t(r1) %*% (T5 + T6 - T7 - T8)
- t(T9) + n * (- 4 * T1/as.numeric(1 + x %*% T1) + 2 * T2/as.numeric(1/M[j, j] +
crossprod(tmp2, T2)) + 2 * T3/as.numeric(1/invMU[j, j] + crossprod(tmp3, T3)))
}
res <- stats::optim(Ginit[j, ], fobj, gobj, method = "BFGS")
Ginit[j, ] <- res$par
g <- as.matrix(Ginit[j, ])
t4 <- t4 + tcrossprod(g, g)
GUGt2 <- g + t2
GUG <- GUG + tcrossprod(GUGt2, GUGt2) * M[j, j]
GVGt2 <- g + t3
GVG <- GVG + tcrossprod(GVGt2, GVGt2) * invMU[j, j]
}
GX <- X %*% Ginit
fit <- stats::glm(Y ~ GX, family = stats::poisson)
mu <- as.vector(fit$coefficients[1])
eta <- as.matrix(fit$coefficients[2 : (u + 1)])
beta <- Ginit %*% eta
theta <- mu + X %*% beta
c.theta <- - exp(theta)
c.theta.mean <- mean(c.theta)
weight <- c.theta / c.theta.mean
V <- theta + ((Y - c.theta) / weight)
wx <- sweep(X, 1, weight, '*')
mean.wx <- colMeans(wx)
wxx <- sweep(X, 2, mean.wx, '-')
sigwx <- crossprod(wxx, sweep(wxx, 1, weight, '*')) / n
wv <- weight * V
mean.wv <- mean(wv)
wvv <- as.matrix(V - mean.wv)
sigwxv <- crossprod(wxx, sweep(wvv, 1, weight, '*')) / n
e1 <- eigen(t(Ginit) %*% M %*% Ginit)
e2 <- eigen(t(Ginit) %*% invMU %*% Ginit)
e3 <- eigen(crossprod(Ginit))
e4 <- matrix(1, n, 1) %*% mu + X %*% beta
e5 <- crossprod(Y, e4) - colSums(exp(e4))
obj2 <- - n/2 * (sum(log(e1$values)) + sum(log(e2$values))) + sum(log(e3$values)) + e5
if (abs(obj1 - obj2) < ftol) {
break
}
else {
obj1 <- obj2
i <- i + 1
}
}
a <- qr(Ginit)
Gammahat <- qr.Q(a)
GX <- X %*% Gammahat
fit <- stats::glm(Y ~ GX, family = stats::poisson)
mu <- as.vector(fit$coefficients[1])
eta <- matrix(fit$coefficients[2 : (u + 1)])
beta <- Gammahat %*% eta
theta <- matrix(1, n, 1) %*% mu + X %*% beta
c.theta <- - exp(theta)
c.theta.mean <- mean(c.theta)
weight <- c.theta / c.theta.mean
wx <- sweep(X, 1, weight, '*')
mean.wx <- colMeans(wx)
wxx <- sweep(X, 2, mean.wx, '-')
sigwx <- crossprod(wxx, sweep(wxx, 1, weight, '*')) / n
inv.sigwx <- chol2inv(chol(sigwx))
var <- inv.sigwx / (- c.theta.mean)
e1 <- eigen(t(Gammahat) %*% M %*% Gammahat)
e2 <- eigen(t(Gammahat) %*% invMU %*% Gammahat)
e3 <- eigen(M)
e4 <- mu + X %*% Gammahat %*% eta
e5 <- crossprod(Y, e4) - colSums(exp(e4))
Gamma0hat <- qr.Q(a, complete = TRUE)[, (u + 1):p]
objfun <- - n/2 * (sum(log(e1$values)) + sum(log(e2$values))) + sum(log(e3$values)) + e5
Gammahat <- as.matrix(Gammahat)
Gamma0hat <- as.matrix(Gamma0hat)
}
return(list(Gammahat = Gammahat, Gamma0hat = Gamma0hat, muhat = mu,
etahat = eta, weighthat = weight, Vhat = V,
avar = var, objfun = objfun))
}
Try the Renvlp package in your browser
Any scripts or data that you put into this service are public.
Renvlp documentation built on Oct. 11, 2023, 1:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions pois.envMU
pois.envMU <- function(X, Y, u, initial = NULL){
X <- as.matrix(X)
Y <- as.matrix(Y)
a <- dim(Y)
n <- a[1]
r <- a[2]
p <- ncol(X)
if (a[1] != nrow(X))
stop("X and Y should have the same number of observations.")
if (u > p || u < 0)
stop("u must be an integer between 0 and r.")
if (sum(duplicated(cbind(X, Y), MARGIN = 2)) > 0)
stop("Some responses also appear in the predictors, or there maybe duplicated columns in X or Y.")
options(warn = -1)
fit <- stats::glm(Y ~ X, family = stats::poisson)
mu <- as.vector(fit$coefficients[1])
beta <- as.vector(fit$coefficients[2 : (p + 1)])
theta <- mu + X %*% beta
c.theta <- - exp(theta)
c.theta.mean <- mean(c.theta)
weight <- c.theta / c.theta.mean
V <- theta + ((Y - c.theta) / weight)
wx <- sweep(X, 1, weight, '*')
mean.wx <- colMeans(wx)
wxx <- sweep(X, 2, mean.wx, '-')
sigwx <- crossprod(wxx, sweep(wxx, 1, weight, '*')) / n
wv <- weight * V
mean.wv <- mean(wv)
wvv <- as.matrix(V - mean.wv)
sigwxv <- crossprod(wxx, sweep(wvv, 1, weight, '*')) / n
inv.sigwx <- chol2inv(chol(sigwx))
M.init <- inv.sigwx / (- c.theta.mean)
U.init <- tcrossprod(beta)
tmp1 <- envMU(M.init, U.init, u, initial = initial)
gamma.init <- tmp1$Gammahat
gamma0.init <- tmp1$Gamma0hat
sigX <- stats::cov(X) * (n - 1) / n
invsigX <- chol2inv(chol(sigX))
if (u == 0) {
Gammahat <- NULL
Gamma0hat <- diag(p)
tmp.M <- eigen(sigX)
mu <- log(mean(Y))
muh <- matrix(1, n, 1) %*% mu
Cn1 <- t(Y) %*% muh - colSums(exp(muh))
eta <- NULL
var <- NULL
weight <- rep(1, n)
V <- Y
objfun <- Cn1 - n/2 * sum(log(tmp.M$values))
} else if (u == p) {
Gammahat <- diag(p)
Gamma0hat <- NULL
tmp.M <- eigen(sigX)
eta <- beta
var <- inv.sigwx / (- c.theta.mean)
Cn1 <- t(Y) %*% theta - colSums(exp(theta))
objfun <- Cn1 - n/2 * sum(log(tmp.M$values))
} else if (u == p - 1) {
maxiter = 1000
ftol = 0.001
M <- sigX
MU <- sigX
tmp.MU <- eigen(MU)
invMU <- sweep(tmp.MU$vectors, MARGIN = 2, 1/tmp.MU$values, "*") %*% t(tmp.MU$vectors)
Cn1 <- t(Y) %*% theta - colSums(exp(theta))
e1 <- eigen(crossprod(gamma.init, M) %*% gamma.init)
e2 <- eigen(crossprod(gamma.init, invMU) %*% gamma.init)
e3 <- eigen(sigX)
temp1 <- sum(log(e1$values)) + sum(log(e2$values)) + sum(log(e3$values))
obj1 <- Cn1 - (n / 2) * temp1
init <- gamma.init
GEidx <- GE(init)
Ginit <- init %*% solve(init[GEidx[1:u], ])
GiX <- X %*% Ginit
fit <- stats::glm(Y ~ GiX, family = stats::poisson)
eta <- as.matrix(fit$coefficients[2 : (u + 1)])
j <- GEidx[p]
g <- as.matrix(Ginit[j, ])
g1 <- Ginit[-j, ]
t2 <- crossprod(g1, as.matrix(M[-j, j]))/M[j, j]
t3 <- crossprod(g1, as.matrix(invMU[-j, j]))/invMU[j, j]
GUGt2 <- g + t2
GUG <- crossprod(Ginit, (M %*% Ginit)) - tcrossprod(GUGt2, GUGt2) * M[j, j]
GVGt2 <- g + t3
GVG <- crossprod(Ginit, (invMU %*% Ginit)) - tcrossprod(GVGt2, GVGt2) * invMU[j, j]
invC1 <- chol2inv(chol(GUG))
invC2 <- chol2inv(chol(GVG))
X1 <- X[ , j]
X2 <- X[ , -j]
t5 <- X1 %*% t(g) %*% eta
t6 <- matrix(1, n, 1) %*% mu + X2 %*% g1 %*% eta
t7 <- t6 + t5
t8 <- t(Y) %*% t7
et6 <- exp(t7)
t9 <- colSums(et6)
fobj <- function(x) {
tmp2 <- x + t2
tmp3 <- x + t3
Ginit[j, ] <- x
GiX <- X %*% Ginit
fit <- stats::glm(Y ~ GiX, family = stats::poisson)
eta <- as.matrix(fit$coefficients[2 : (u + 1)])
tmp4 <- X1 %*% t(x) %*% eta + t6
T2 <- invC1 %*% tmp2
T3 <- invC2 %*% tmp3
T4 <- exp(tmp4)
n /2 * (- 2 * log(1 + sum(x^2)) + log(1 + M[j, j] *
crossprod(tmp2,T2)) + log(1 + invMU[j, j] *
crossprod(tmp3,T3))) - t(Y) %*% tmp4 + colSums(T4)
}
gobj <- function(x) {
tmp2 <- x + t2
tmp3 <- x + t3
Ginit[j, ] <- x
GiX <- X %*% Ginit
fit <- stats::glm(Y ~ GiX, family = stats::poisson)
eta <- as.matrix(fit$coefficients[2 : (u + 1)])
tmp4 <- X1 %*% t(x) %*% eta + t6
Cn <- Y - exp(tmp4)
T2 <- invC1 %*% tmp2
T3 <- invC2 %*% tmp3
A1 <- crossprod(g1, sigwx[-j, -j]) %*% g1 + as.matrix(x) %*% sigwx[j, -j] %*% g1 + crossprod(g1, sigwx[-j, j]) %*% t(x) + tcrossprod(x) * sigwx[j , j]
invC3 <- chol2inv(chol(A1))
T5 <- crossprod(X, Cn) %*% t(eta)
tmp5 <- X1 %*% t(x) + X2 %*% g1
T6 <- tcrossprod(sigwxv, Cn) %*% tmp5 %*% invC3
tmp6 <- as.matrix(sigwx[ , -j]) %*% g1 + tcrossprod(sigwx[ , j], x)
T7 <- tmp6 %*% invC3 %*% crossprod(tmp5, Cn) %*% t(eta)
T8 <- tmp6 %*% eta %*% crossprod(Cn, tmp5) %*% invC3
r1 <- rep(0, p)
r1[j] <- 1
T9 <- t(r1) %*% (T5 + T6 - T7 - T8)
- t(T9) + n * (- 4 * x %*% solve(1 + sum(x^2)) + 2 * T2/as.numeric(1/M[j, j] +
crossprod(tmp2, T2)) + 2 * T3/as.numeric(1/invMU[j, j] + crossprod(tmp3, T3)))
}
i <- 1
while (i < maxiter) {
res <- stats::optim(Ginit[j, ], fobj, gobj, method = "BFGS")
Ginit[j, ] <- res$par
GX <- X %*% Ginit
fit <- stats::glm(Y ~ GX, family = stats::poisson)
mu <- as.vector(fit$coefficients[1])
eta <- as.matrix(fit$coefficients[2 : (u + 1)])
beta <- Ginit %*% eta
theta <- mu + X %*% beta
c.theta <- - exp(theta)
c.theta.mean <- mean(c.theta)
weight <- c.theta / c.theta.mean
V <- theta + ((Y - c.theta) / weight)
wx <- sweep(X, 1, weight, '*')
mean.wx <- colMeans(wx)
wxx <- sweep(X, 2, mean.wx, '-')
sigwx <- crossprod(wxx, sweep(wxx, 1, weight, '*')) / n
wv <- weight * V
mean.wv <- mean(wv)
wvv <- as.matrix(V - mean.wv)
sigwxv <- crossprod(wxx, sweep(wvv, 1, weight, '*')) / n
e1 <- eigen(t(Ginit) %*% M %*% Ginit)
e2 <- eigen(t(Ginit) %*% invMU %*% Ginit)
e3 <- eigen(crossprod(Ginit))
e4 <- mu + X %*% beta
e5 <- crossprod(Y, e4) - colSums(exp(e4))
obj2 <- - n/2 * (sum(log(e1$values)) + sum(log(e2$values))) + sum(log(e3$values)) + e5
if (abs(obj1 - obj2) < ftol) {
break
}
else {
obj1 <- obj2
i <- i + 1
}
}
a <- qr(Ginit)
Gammahat <- qr.Q(a)
GX <- X %*% Gammahat
fit <- stats::glm(Y ~ GX, family = stats::poisson)
mu <- as.vector(fit$coefficients[1])
eta <- matrix(fit$coefficients[2 : (u + 1)])
beta <- Gammahat %*% eta
theta <- mu + X %*% beta
c.theta <- - exp(theta)
c.theta.mean <- mean(c.theta)
weight <- c.theta / c.theta.mean
wx <- sweep(X, 1, weight, '*')
mean.wx <- colMeans(wx)
wxx <- sweep(X, 2, mean.wx, '-')
sigwx <- crossprod(wxx, sweep(wxx, 1, weight, '*')) / n
inv.sigwx <- chol2inv(chol(sigwx))
var <- inv.sigwx / (- c.theta.mean)
e1 <- eigen(t(Gammahat) %*% M %*% Gammahat)
e2 <- eigen(t(Gammahat) %*% invMU %*% Gammahat)
e3 <- eigen(M)
e4 <- mu + X %*% Gammahat %*% eta
e5 <- crossprod(Y, e4) - colSums(exp(e4))
Gamma0hat <- qr.Q(a, complete = TRUE)[, (u + 1):r, drop = FALSE]
objfun <- - n/2 * (sum(log(e1$values)) + sum(log(e2$values))) + sum(log(e3$values)) + e5
Gammahat <- as.matrix(Gammahat)
Gamma0hat <- as.matrix(Gamma0hat)
} else {
maxiter = 100
ftol = 0.001
M <- sigX
MU <- sigX
tmp.MU <- eigen(MU)
invMU <- sweep(tmp.MU$vectors, MARGIN = 2, 1/tmp.MU$values, "*") %*% t(tmp.MU$vectors)
Cn1 <- t(Y) %*% theta - colSums(exp(theta))
e1 <- eigen(crossprod(gamma.init, M) %*% gamma.init)
e2 <- eigen(crossprod(gamma.init, invMU) %*% gamma.init)
e3 <- eigen(sigX)
temp1 <- sum(log(e1$values)) + sum(log(e2$values)) + sum(log(e3$values))
obj1 <- Cn1 - (n / 2) * temp1
init <- gamma.init
GEidx <- GE(init)
Ginit <- init %*% solve(init[GEidx[1:u], ])
GiX <- X %*% Ginit
fit <- stats::glm(Y ~ GiX, family = stats::poisson)
eta <- as.matrix(fit$coefficients[2 : (u + 1)])
GUG <- crossprod(Ginit, (M %*% Ginit))
GVG <- crossprod(Ginit, (invMU %*% Ginit))
t4 <- crossprod(Ginit[GEidx[(u + 1):p], ], Ginit[GEidx[(u + 1):p], ]) + diag(u)
i <- 1
while (i < maxiter) {
for (j in GEidx[(u + 1):p]) {
g <- as.matrix(Ginit[j, ])
g1 <- as.matrix(Ginit[-j, ])
t2 <- crossprod(g1, as.matrix(M[-j, j]))/M[j, j]
t3 <- crossprod(g1, as.matrix(invMU[-j, j]))/invMU[j, j]
GUGt2 <- g + t2
GUG <- GUG - tcrossprod(GUGt2, GUGt2) * M[j, j]
GVGt2 <- g + t3
GVG <- GVG - tcrossprod(GVGt2, GVGt2) * invMU[j, j]
t4 <- t4 - tcrossprod(g, g)
invC1 <- chol2inv(chol(GUG))
invC2 <- chol2inv(chol(GVG))
invt4 <- chol2inv(chol(t4))
X1 <- X[ , j]
X2 <- X[ , -j]
t5 <- X1 %*% t(g) %*% eta
t6 <- mu + X2 %*% g1 %*% eta
t7 <- t5 + t6
t8 <- t(Y) %*% t7
et7 <- exp(t7)
t9 <- colSums(et7)
fobj <- function(x) {
tmp2 <- x + t2
tmp3 <- x + t3
Ginit[j, ] <- x
GiX <- X %*% Ginit
fit <- stats::glm(Y ~ GiX, family = stats::poisson)
eta <- as.matrix(fit$coefficients[2 : (u + 1)])
tmp4 <- X1 %*% t(x) %*% eta + t6
T1 <- invt4 %*% x
T2 <- invC1 %*% tmp2
T3 <- invC2 %*% tmp3
T4 <- exp(tmp4)
n /2 * (-2 * log(1 + x %*% T1) + log(1 + M[j, j] *
crossprod(tmp2, T2)) + log(1 + invMU[j, j] *
crossprod(tmp3, T3))) - t(Y) %*% tmp4 + colSums(T4)
}
gobj <- function(x) {
tmp2 <- x + t2
tmp3 <- x + t3
Ginit[j, ] <- x
GiX <- X %*% Ginit
fit <- stats::glm(Y ~ GiX, family = stats::poisson)
eta <- as.matrix(fit$coefficients[2 : (u + 1)])
tmp4 <- X1 %*% t(x) %*% eta + t6
Cn <- Y - exp(tmp4)
T1 <- invt4 %*% x
T2 <- invC1 %*% tmp2
T3 <- invC2 %*% tmp3
A1 <- crossprod(g1, sigwx[-j, -j]) %*% g1 + as.matrix(x) %*% sigwx[j, -j] %*% g1 + crossprod(g1, sigwx[-j, j]) %*% t(x) + tcrossprod(x) * sigwx[j , j]
invC3 <- chol2inv(chol(A1))
T5 <- crossprod(X, Cn) %*% t(eta)
tmp5 <- X1 %*% t(x) + X2 %*% g1
T6 <- tcrossprod(sigwxv, Cn) %*% tmp5 %*% invC3
tmp6 <- as.matrix(sigwx[ , -j]) %*% g1 + tcrossprod(sigwx[ , j], x)
T7 <- tmp6 %*% invC3 %*% crossprod(tmp5, Cn) %*% t(eta)
T8 <- tmp6 %*% eta %*% crossprod(Cn, tmp5) %*% invC3
r1 <- rep(0, p)
r1[j] <- 1
T9 <- t(r1) %*% (T5 + T6 - T7 - T8)
- t(T9) + n * (- 4 * T1/as.numeric(1 + x %*% T1) + 2 * T2/as.numeric(1/M[j, j] +
crossprod(tmp2, T2)) + 2 * T3/as.numeric(1/invMU[j, j] + crossprod(tmp3, T3)))
}
res <- stats::optim(Ginit[j, ], fobj, gobj, method = "BFGS")
Ginit[j, ] <- res$par
g <- as.matrix(Ginit[j, ])
t4 <- t4 + tcrossprod(g, g)
GUGt2 <- g + t2
GUG <- GUG + tcrossprod(GUGt2, GUGt2) * M[j, j]
GVGt2 <- g + t3
GVG <- GVG + tcrossprod(GVGt2, GVGt2) * invMU[j, j]
}
GX <- X %*% Ginit
fit <- stats::glm(Y ~ GX, family = stats::poisson)
mu <- as.vector(fit$coefficients[1])
eta <- as.matrix(fit$coefficients[2 : (u + 1)])
beta <- Ginit %*% eta
theta <- mu + X %*% beta
c.theta <- - exp(theta)
c.theta.mean <- mean(c.theta)
weight <- c.theta / c.theta.mean
V <- theta + ((Y - c.theta) / weight)
wx <- sweep(X, 1, weight, '*')
mean.wx <- colMeans(wx)
wxx <- sweep(X, 2, mean.wx, '-')
sigwx <- crossprod(wxx, sweep(wxx, 1, weight, '*')) / n
wv <- weight * V
mean.wv <- mean(wv)
wvv <- as.matrix(V - mean.wv)
sigwxv <- crossprod(wxx, sweep(wvv, 1, weight, '*')) / n
e1 <- eigen(t(Ginit) %*% M %*% Ginit)
e2 <- eigen(t(Ginit) %*% invMU %*% Ginit)
e3 <- eigen(crossprod(Ginit))
e4 <- matrix(1, n, 1) %*% mu + X %*% beta
e5 <- crossprod(Y, e4) - colSums(exp(e4))
obj2 <- - n/2 * (sum(log(e1$values)) + sum(log(e2$values))) + sum(log(e3$values)) + e5
if (abs(obj1 - obj2) < ftol) {
break
}
else {
obj1 <- obj2
i <- i + 1
}
}
a <- qr(Ginit)
Gammahat <- qr.Q(a)
GX <- X %*% Gammahat
fit <- stats::glm(Y ~ GX, family = stats::poisson)
mu <- as.vector(fit$coefficients[1])
eta <- matrix(fit$coefficients[2 : (u + 1)])
beta <- Gammahat %*% eta
theta <- matrix(1, n, 1) %*% mu + X %*% beta
c.theta <- - exp(theta)
c.theta.mean <- mean(c.theta)
weight <- c.theta / c.theta.mean
wx <- sweep(X, 1, weight, '*')
mean.wx <- colMeans(wx)
wxx <- sweep(X, 2, mean.wx, '-')
sigwx <- crossprod(wxx, sweep(wxx, 1, weight, '*')) / n
inv.sigwx <- chol2inv(chol(sigwx))
var <- inv.sigwx / (- c.theta.mean)
e1 <- eigen(t(Gammahat) %*% M %*% Gammahat)
e2 <- eigen(t(Gammahat) %*% invMU %*% Gammahat)
e3 <- eigen(M)
e4 <- mu + X %*% Gammahat %*% eta
e5 <- crossprod(Y, e4) - colSums(exp(e4))
Gamma0hat <- qr.Q(a, complete = TRUE)[, (u + 1):p]
objfun <- - n/2 * (sum(log(e1$values)) + sum(log(e2$values))) + sum(log(e3$values)) + e5
Gammahat <- as.matrix(Gammahat)
Gamma0hat <- as.matrix(Gamma0hat)
}
return(list(Gammahat = Gammahat, Gamma0hat = Gamma0hat, muhat = mu,
etahat = eta, weighthat = weight, Vhat = V,
avar = var, objfun = objfun))
}
Try the Renvlp package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.