Nothing
## Objective value of multiple regression
##' @title obj.multiRegression
##' @param Xs eQTL matrices
##' @param Ys gene expression matrices
##' @param fit regression fit result object
##' @param trans if rows for sample, trans = TRUE, otherwise, trans = FALSE. Default FALSE
##' @return error squared norm of ||(I-B)Y - FX||_2^2
obj.multiRegression = function(Xs, Ys, fit, trans = F) {
m = length(Ys)
if (is.matrix(Xs)) {
Xs = lapply(1:m, function(i) {
Xs
})
}
if (!trans) {
Xs = lapply(1:m, function(i) {
t(Xs[[i]])
})
Ys = lapply(1:m, function(i) {
t(Ys[[i]])
})
}
error = .Call("ObjMultiReg", Xs, Ys, fit, PACKAGE = "fssemR")
error
}
## center Xs and Ys
##' @title proc.centerFSSEM
##' @param Xs eQTL matrices
##' @param Ys list of gene expression matrices
##' @return centered Xs and Ys and mean vectors
proc.centerFSSEM = function(Xs, Ys) {
mX = rowMeans(Xs)
mY = lapply(Ys, rowMeans)
center = function(m) {
apply(m, 1, function(x){ x - mean(x) })
}
Xs = center(Xs)
Ys = lapply(Ys, center)
list(Xs = t(Xs),
Ys = lapply(Ys, t),
mX = mX, mY = mY)
}
## logliklihood of FSSEM
##' @title logLikFSSEM
##' @param Bs Network matrices
##' @param Wl Weights for lasso term
##' @param Wf Weights for fused term
##' @param lambda Hyperparameter of lasso term
##' @param rho Hyperparameter of fused lasso term
##' @param sigma2 noise variance
##' @param Dets determinants of I-B matrices
##' @param n number of observations
##' @param p number of genes
##' @return objective value of FSSEM with specified hyper-paramters
logLikFSSEM = function(Bs, Wl, Wf, lambda, rho, sigma2, Dets, n, p) {
K = length(Bs)
loglik = 0
l1 = 0
lf = 0
rho = min(rho, 1e12)
lf = rho * Wf * abs(Bs[[2]] - Bs[[1]])
for (k in 1:K) {
l1 = l1 + lambda * Wl[[k]] * abs(Bs[[k]])
loglik = loglik - n[k] / 2 * log(Dets[k]**2)
}
diag(l1) = 0
diag(lf) = 0
loglik + p * sum(n) / 2 * log(sigma2) + sum(l1) + sum(lf)
}
##' function generator function
##' @title cwiseGradient4FSSEM
##' @param n number of observations
##' @param c cofactor vector
##' @param Y Matrix of gene expression
##' @param R Residual matrix
##' @param Y2norm Column of YtY
##' @param sigma2 noise variance
##' @return function whose argument is column vector bi
cwiseGradient4FSSEM = function(n, c, Y, R, Y2norm, sigma2) {
function(x) {
n * t(c) + (Y2norm * x - tcrossprod(R, Y)) / sigma2
}
}
## lipschitz computation
# cwiseLipschitzFSSEMv0 = function(n, o, g, s, c2, Deti, Y2norm, sigma2, p) {
# oxg = tcrossprod(o, g)
# Imo = diag(p - 1) - g %*% oxg
# sog = s %*% oxg
# c = 1 - s %*% tcrossprod(o, s)
# zta = 1e-12
# x = -1 * tcrossprod(chol2inv(Imo + diag(p - 1) * zta), sog)
# Li = n * c2 / (crossprod(x, Imo %*% x) + 2 * sog %*% x + c) / (Deti + 1e-16) + Y2norm / sigma2
# while(Li < 0) {
# zta = zta * 10
# x = -1 * tcrossprod(chol2inv(Imo + diag(p - 1) * zta), sog)
# Li = n * c2 / (crossprod(x, Imo %*% x) + 2 * sog %*% x + c) / (Deti + 1e-16) + Y2norm / sigma2
# }
# Li
#}
## speed-up lipschitz computation
cwiseLipschitz4FSSEM = function(n, z, c2, b2, Y2norm, sigma2) {
n * c2 / (abs(z) - sqrt(c2 * b2))**2 + Y2norm / sigma2
}
## speed-up lipschitz computation
## abs(z - c %*% b) >= abs(z) - abs(c %*% b) >= abs(z) - sqrt(c2 * b2)
cwiseLipschitz4FSSEM2B = function(n, z, c2, b2, Y2norm, sigma2, ImB, i) {
if (abs(z) - sqrt(c2 * b2) > 0) {
n * c2 / (abs(z) - sqrt(c2 * b2))**2 + Y2norm / sigma2
} else {
## n * c2 / (abs(z) - sqrt(c2 * b2))**2 + Y2norm / sigma2
cwiseLipschitzFSSEMv0(n, c2, Y2norm, sigma2, ImB, i)[1]
## Inf
}
}
cwiseLipschitzFSSEMv0 = function(n, c2, Y2norm, sigma2, ImB, i) {
## x^tWx + 2R^tx + C
p = nrow(ImB)
e = 1e-6
I = chol2inv(chol(crossprod(ImB[, -i]) + diag(e, p - 1)))
D = det(crossprod(ImB[, -i]) + diag(e, p - 1))
W = diag(p - 1) - ImB[-i, -i] %*% tcrossprod(I, ImB[-i, -i])
R = ImB[-i, -i] %*% tcrossprod(I, ImB[i, -i, drop = FALSE])
C = 1 - ImB[i, -i, drop = FALSE] %*% tcrossprod(I, ImB[i, -i, drop = FALSE])
v = solve(W + diag(e, p - 1), -1 * R)
L = n * c2 / (crossprod(v, I %*% v) + 2 * crossprod(R, v) + C) / (D + e) + Y2norm / sigma2
if (L < 0) {
L = Inf
}
L
}
## Proximal operator for FLSA
proxFLSA = function(u, w, r, lambda, rho, c) {
D = u[[1]] - u[[2]]
eQ = (abs(D) <= 2 * rho * r / c)
Df = 1 - eQ
rho = min(rho, 1e16)
x = list(
(u[[1]] + u[[2]]) / 2 * eQ + (u[[1]] - sign(D) * rho * r / c) * Df,
(u[[1]] + u[[2]]) / 2 * eQ + (u[[2]] + sign(D) * rho * r / c) * Df
)
lapply(1:2, function(k) {
xe = pmax(x[[k]] - lambda * (w[[1]] + w[[2]]) / 2 / c, 0) + pmin(x[[k]] + lambda * (w[[1]] + w[[2]]) / 2 / c, 0)
xd = pmax(x[[k]] - lambda * w[[k]] / c, 0) + pmin(x[[k]] + lambda * w[[k]] / c, 0)
xe * eQ + xd * Df
})
}
## Proximal operator for PNSA
proxPNSA = function(u, w, r, lambda, rho, c) {
p = length(u[[1]])
D = cbind(diag(p), -diag(p))
U = matrix(c(u[[1]], u[[2]]), nrow = 1)
v = tcrossprod(D, U)[,1]
eQ = (norm(v, "2")[1] <= 2 * rho * r / c)
Df = 1 - eQ
rho = min(rho, 1e16)
x = if (eQ) {
list((u[[1]] + u[[2]]) / 2,
(u[[1]] + u[[2]]) / 2)
} else {
v = rho * r / c * v / norm(v, "2")[1]
list(u[[1]] - v, u[[2]] + v)
}
lapply(1:2, function(k) {
xe = pmax(x[[k]] - lambda * (w[[1]] + w[[2]]) / 2 / c, 0) + pmin(x[[k]] + lambda * (w[[1]] + w[[2]]) / 2 / c, 0)
xd = pmax(x[[k]] - lambda * w[[k]] / c, 0) + pmin(x[[k]] + lambda * w[[k]] / c, 0)
xe * eQ + xd * Df
})
}
sigma2FSSEM = function(Xs, Ys, Bs, Fs, n, p, m) {
error = 0
for(k in 1:m) {
error = error + (norm(Ys[[k]] - Bs[[k]] %*% Ys[[k]] - Fs[[k]] %*% Xs, "f"))**2
}
error / (sum(n) * p)
}
inert_opt = function(opts = c("continuous", "linear"), init = 0) {
switch(
opts,
"continuous" = function(k) {
init
},
"linear" = function(k) {
(k - 1) / (k + 2)
}
)
}
##' @title inverseB
##' @description inverse matrices of B network for adaptive FSSEM
##' @param Bs list of network matrices
##' @return list of inversed B matrices
##' @export
inverseB = function(Bs) {
lapply(Bs, function(B) { 1 / abs(B) })
}
##' @title invoneB
##' @description if you do not want to get inversed B matrces, invoneB gives you a matrix with constant 1 instead in FSSEM
##' @param Bs list of network matrices
##' @return list of invone B matrices
##' @export
invoneB = function(Bs) {
lapply(Bs, function(B) { 1 / abs(sign(B)) })
}
##' @title flinvB
##' @description inversed difference of two B matrices. For adaptive fused lasso penalty
##' @param Bs list of network matrices
##' @return inversed difference matrices
##' @export
flinvB = function(Bs) {
1 / abs(Bs[[2]] - Bs[[1]])
}
##' @title floneB
##' @description if you do not want adaptive fused lasso penalty, floneB replace flinvB
##' @param Bs list of network matrices
##' @return matrix whose entries are all 1
##' @export
floneB = function(Bs) {
1 / abs(sign(Bs[[2]] - Bs[[1]]))
}
##' @title pninvB
##' @description inversed column l2 norm for perturbed group lasso penalty
##' @param Bs list of network matrices
##' @return inversed l2 norm of B2 - B1
##' @export
pninvB = function(Bs) {
apply(Bs[[1]] - Bs[[2]], 2, function(x) {
1 / norm(x, "2")
})
}
##' @title pnoneB
##' @description if you do not want adaptive group lasso penalty, pnoneB replace pninvB
##' @param Bs list of network matrices
##' @return inversed l2 norm of B2 - B1 with all entries is 1
##' @export
pnoneB = function(Bs) {
apply(Bs[[1]] - Bs[[2]], 2, function(x) {
1
})
}
logLiklihood = function(Xs, Ys, Bs, Fs, mu, Dets, sigma2, p) {
err = 0
logl = 0
m = length(Ys)
n = sapply(Ys, ncol)
df = 0
for(i in 1:m) {
err = err + norm(Ys[[i]] - Bs[[i]] %*% Ys[[i]] - Fs[[i]] %*% Xs - tcrossprod(mu[[i]], rep(1, n[i])), "f")**2
logl = logl - n[i] / 2 * log(Dets[i]**2)
df = df + sum(Bs[[2]] - Bs[[1]] != 0 & Bs[[i]] != 0) + p
}
df = df + sum(Bs[[2]] - Bs[[1]] == 0 & Bs[[1]] != 0) + 1
sigma2 = err / (sum(n) * p)
2 * (logl + sum(n) * p / 2 * log(2 * pi * sigma2)) + df * (log(n[1]) + log(n[2]))
}
bayesianInfocriterion = function(Xs, Ys, Bs, Fs, mu, Dets, sigma2, p) {
logl = 0
m = length(Ys)
n = sapply(Ys, ncol)
df = 0
for(i in 1:m) {
err = norm(Ys[[i]] - Bs[[i]] %*% Ys[[i]] - Fs[[i]] %*% Xs - tcrossprod(mu[[i]], rep(1, n[i])), "f")**2
logl = logl - n[i] / 2 * log(Dets[i]**2) + n[i] * p / 2 * log(sigma2) + err / (2 * sigma2)
df = df + sum(Bs[[2]] - Bs[[1]] != 0 & Bs[[i]] != 0) + p
}
df = df + sum(Bs[[2]] - Bs[[1]] == 0 & Bs[[1]] != 0) + 1
2 * logl + df * (log(n[1]) + log(n[2]))
}
L2lamax = function(Xs, Ys, Sk, n, p, k) {
m = length(Ys)
Xs = lapply(1:m, function(i) {
t(Xs[[i]])
})
Ys = lapply(1:m, function(i) {
t(Ys[[i]])
})
lambda = .Call("L2lamax", Xs, Ys, Sk, n, p, k, PACKAGE = "fssemR")
lambda
}
##' @title TPR
##' @description Power of detection for network prediction
##' @param X list of predicted network matrices
##' @param B list of true network matrices
##' @param PREC precision threshold for FDR test. Default 0.
##' @export
TPR = function(X, B, PREC = 0) {
X = as.matrix(X)
B = as.matrix(B)
sum(abs(X[B != 0]) > PREC) / sum(B != 0)
}
##' @title FDR
##' @description False discovery rate for network prediction
##' @param X list of predicted network matrices
##' @param B list of true network matrices
##' @param PREC precision threshold for FDR test. Default 0.
##' @export
FDR = function(X, B, PREC = 0) {
X = as.matrix(X)
B = as.matrix(B)
# PREC = max(PREC, 0.05)
if(sum(X != 0) != 0) {
sum(abs(X[B == 0]) > PREC) / sum(X != 0)
} else {
0
}
}
## center Xs and Ys
##' @title proc.centerFSSEM2
##' @param Xs list of eQTL matrices
##' @param Ys list of gene expression matrices
##' @return centered Xs and Ys and mean vectors
proc.centerFSSEM2 = function(Xs, Ys) {
mX = lapply(Xs, rowMeans)
mY = lapply(Ys, rowMeans)
center = function(m) {
apply(m, 1, function(x){ x - mean(x) })
}
Xs = lapply(Xs, center)
Ys = lapply(Ys, center)
list(Xs = lapply(Xs, t),
Ys = lapply(Ys, t),
mX = mX, mY = mY)
}
sigma2FSSEM2 = function(Xs, Ys, Bs, Fs, n, p, m) {
error = 0
for(k in 1:m) {
error = error + (norm(Ys[[k]] - Bs[[k]] %*% Ys[[k]] - Fs[[k]] %*% Xs[[k]], "f"))**2
}
error / (sum(n) * p)
}
logLiklihood2 = function(Xs, Ys, Bs, Fs, mu, Dets, sigma2, p) {
err = 0
logl = 0
m = length(Ys)
n = sapply(Ys, ncol)
df = 0
for(i in 1:m) {
err = err + norm(Ys[[i]] - Bs[[i]] %*% Ys[[i]] - Fs[[i]] %*% Xs[[i]] - tcrossprod(mu[[i]], rep(1, n[i])), "f")**2
logl = logl - n[i] / 2 * log(Dets[i]**2)
df = df + sum(Bs[[2]] - Bs[[1]] != 0 & Bs[[i]] != 0) + p
}
df = df + sum(Bs[[2]] - Bs[[1]] == 0 & Bs[[1]] != 0) + 1
sigma2 = err / (sum(n) * p)
2 * (logl + sum(n) * p / 2 * log(2 * pi * sigma2)) + df * (log(n[1]) + log(n[2]))
}
bayesianInfocriterion2 = function(Xs, Ys, Bs, Fs, mu, Dets, sigma2, p) {
logl = 0
m = length(Ys)
n = sapply(Ys, ncol)
df = 0
for(i in 1:m) {
err = norm(Ys[[i]] - Bs[[i]] %*% Ys[[i]] - Fs[[i]] %*% Xs[[i]] - tcrossprod(mu[[i]], rep(1, n[i])), "f")**2
logl = logl - n[i] / 2 * log(Dets[i]**2) + n[i] * p / 2 * log(sigma2) + err / (2 * sigma2)
df = df + sum(Bs[[2]] - Bs[[1]] != 0 & Bs[[i]] != 0) + p
}
df = df + sum(Bs[[2]] - Bs[[1]] == 0 & Bs[[1]] != 0) + 1
2 * logl + df * (log(n[1]) + log(n[2]))
}
## BIC for node-based penalty
bayesianInfocriterion3 = function(Xs, Ys, Bs, Fs, mu, Dets, sigma2, p) {
logl = 0
m = length(Ys)
n = sapply(Ys, ncol)
df = 0
for(i in 1:m) {
err = norm(Ys[[i]] - Bs[[i]] %*% Ys[[i]] - Fs[[i]] %*% Xs[[i]] - tcrossprod(mu[[i]], rep(1, n[i])), "f")**2
logl = logl - n[i] / 2 * log(Dets[i]**2) + n[i] * p / 2 * log(sigma2) + err / (2 * sigma2)
df = df + sum(colSums((Bs[[2]] - Bs[[1]])**2) != 0 & colSums(Bs[[i]]**2) != 0) * p + p
}
df = df + sum(colSums((Bs[[2]] - Bs[[1]])**2) == 0 & colSums(Bs[[1]]**2) != 0) * p + 1
2 * logl + df * (log(n[1]) + log(n[2]))
}
## logliklihood of NFSSEM
##' @title logLikNFSSEM
##' @param Bs Network matrices
##' @param Wl Weights for lasso term
##' @param Wf Weights for group perturb lasso term
##' @param lambda Hyperparameter of lasso term
##' @param rho Hyperparameter of group fused lasso term
##' @param sigma2 noise variance
##' @param Dets determinants of I-B matrices
##' @param n number of observations
##' @param p number of genes
##' @return objective value of NFSSEM with specified hyper-paramters
logLikNFSSEM = function(Bs, Wl, Wf, lambda, rho, sigma2, Dets, n, p) {
K = length(Bs)
loglik = 0
l1 = 0
lf = 0
rho = min(rho, 1e12)
lf = rho * Wf * sqrt(colSums((Bs[[2]] - Bs[[1]])**2))
for (k in 1:K) {
l1 = l1 + lambda * Wl[[k]] * abs(Bs[[k]])
loglik = loglik - n[k] / 2 * log(Dets[k]**2)
}
diag(l1) = 0
loglik + p * sum(n) / 2 * log(sigma2) + sum(l1) + sum(lf)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.