Package 'lglasso'

Title: Longitudinal Graphical Lasso
Description: For high-dimensional correlated observations, this package carries out the L_1 penalized maximum likelihood estimation of the precision matrix (network) and the correlation parameters. The correlated data can be longitudinal data (may be irregularly spaced) with dampening correlation or clustered data with uniform correlation. For the details of the algorithms, please see the paper Jie Zhou et al. Identifying Microbial Interaction Networks Based on Irregularly Spaced Longitudinal 16S rRNA sequence data <doi:10.1101/2021.11.26.470159>.
Authors: Jie Zhou [aut, cre, cph], Jiang Gui [aut], Weston Viles [aut], Anne Hoen [aut]
Maintainer: Jie Zhou <[email protected]>
License: GPL (>= 3)
Version: 1.0.0
Built: 2026-05-28 07:17:16 UTC
Source: https://github.com/jiezhou-2/lglasso

Help Index

Title

Description

Title

Usage

AA(
  B,
  data,
  type = c("expFixed"),
  expFix,
  maxit = 30,
  tol = 10^(-4),
  lower = c(0.01, 0.1),
  upper = c(10, 5),
  ...
)

Arguments

B

a p by p given precision matrix
data

a (p+2)-by-n data frame
type

specify how to model the covariance matrix
expFix

the parameter in variance function when the data are longitudinal
maxit

the maximum iteration number
tol

the minimum difference for the algorithm to be called converged
lower

vector of length 1 or 2 which specifies the lower bounds for alpha_1 (and alpha_2) in the correlation matrix
upper

vector of length 1 or 2 which specifies the upper bounds for alpha_1 (and alpha_2) in the correlation matrix
...

other unspecifed parameterss

Value

a list of matrices

Estimate the phimatrix in heterogeneous model

Description

Estimate the phimatrix in heterogeneous model

Usage

AAheter(data, wi, alpha, group, l)

Arguments

data

longitudinal data set
wi

given precision matrix
alpha

expoential distribution with rate alpha
group

specify how data is grouped
l

number of random samples in importance sampling

Value

a list for estimates of tau and AA

Title

Description

Title

Usage

BB(A, data, lambda, type = c("expFixed"), diagonal = TRUE)

Arguments

A

a list representing the phimatrices
data

a (p+2)-by-n data frame
lambda

given tuning parameter(s)
type

specify model type

Value

a list with the same length as A

density function in EM algorithm

Description

density function in EM algorithm

Usage

conDensityTau(tau, expFixed, datai, wi, alpha, groupi)

Arguments

tau

the dampening rate
datai

the data set for subject i
wi

the given precision matrices
alpha

the rate in exponential distribution
groupi

the data point indices

Compute the cross validation error

Description

Compute the cross validation error

Usage

cvError(data.train, data.valid, B, group.train = NULL, group.valid = NULL)

Arguments

data.train

trainng data
data.valid

testing data
B

given network (or network list)
group.train

group in training data
group.valid

group in testing data

Value

a matrix

Cross validation for lglasso

Description

The function computes the cross validation errors for one of the three network models in lglasso command.

Usage

CVlglasso(
  type = c("expFixed"),
  data,
  group = NULL,
  random = FALSE,
  lambda = NULL,
  nlam = 10,
  lam.min.ratio = 0.01,
  K,
  expFix = 1,
  trace = FALSE,
  NN = 500
)

Arguments

type

underlying model type, either general, longihomo or longiheter.
data

raw data
group

group variable
lambda

tuning parameter
nlam

number of tuning parameter
lam.min.ratio

ratio of largest lambda vs smallest lambda
K

cv folds
expFix

given parameter
trace

whether show the process
cores

parallel computing

Value

list of which the first component is the cross validation errors and the second component is the corresponding tuning parameters

start = match.arg(start)

Sminus = S diag(Sminus) = 0 if (is.null(lambda)) if (!((lam.min.ratio <= 1) && (lam.min.ratio > 0))) cat("\nlam.min.ratio must be in (0, 1]... setting to 1e-2!") lam.min.ratio = 0.01

if (!((nlam > 0) && (nlam%%1 == 0))) cat("\nnlam must be a positive integer... setting to 10!") nlam = 10

lam.max = max(abs(Sminus)) lam.min = lam.min.ratio * lam.max lambda = 10^seq(log10(lam.min), log10(lam.max), length = nlam) if (!is.null(group)) lambda=expand.grid(lambda,lambda)

else if (is.null(group)) lambda = sort(lambda)

nnlambda=ifelse(is.null(group),length(lambda),nrow(lambda)) cv_error=matrix(0,nrow=nnlambda,ncol=K)

k=NULL k=1 CV = foreach(k = 1:K, .packages = "lglasso", .combine = "cbind", .inorder = FALSE) %dopar% 

             if (trace) {
               progress = utils::txtProgressBar(max = K, style = 3)
             }

             leave.out =subjects[ind[(1 + floor((k - 1) * n/K)):floor(k *
                                                                        n/K)]]
             indexValid=which(data[,1] %in% leave.out)
             data.train = data[-indexValid, , drop = FALSE]
             data_bar = apply(data.train[,-c(1,2)], 2, mean)
             data.train[,-c(1,2)] = scale(data.train[,-c(1,2)], center = data_bar, scale = FALSE)
             data.valid = data[indexValid,, drop = FALSE]
             data.valid[,-c(1,2)] = scale(data.valid[,-c(1,2)], center = data_bar, scale = FALSE)
             group.train=group[-indexValid]
             group.valid=group[indexValid]
             #S.train = crossprod(data.train[,-c(1,2)])/(dim(data.train)[1])
             #S.valid = crossprod(data.valid[,-c(1,2)])/(dim(data.valid)[1])


             if (is.null(group)){

               aa= sapply(lambda, function(x) {lglasso(data=data.train,lambda=x)$wi})
               cc=lapply(aa, function(B){
                 M=ifelse(abs(B)<=10^(-2), 0,1)
                 diag(M)=2
                 M
               })
               bb=lapply(cc, function(B) {cvError(data.train=data.train,data.valid=data.valid,B=B)})
               bb=do.call(c,bb)
             }
             if (!is.null(group)){
               aa= apply(lambda,1, function(x) {lglasso(data=data.train,lambda=x, expFix = expFix, group=group.train)$wi})
               cc=lapply(aa, function(B){
                 lapply(B, function(Z){
                   M=ifelse(abs(Z)<=10^(-1), 0,1)
                   diag(M)=2
                   M
                 }
                 )
               }
               )

               bb=lapply(cc, function(BB){cvError(data.train=data.train,data.valid=data.valid,B=BB,
                                                  group.valid=group.valid, group.train = group.train)})
               bb=do.call(c,bb)
             }

             cv_error=bb

             if (trace) {
               utils::setTxtProgressBar(progress,  k)
             }

             return(cv_error=cv_error)
           }

stopCluster(cluster) output=structure(list(cv_error=CV, lambda=lambda),class="cvlglasso") return(output)

Cross validation for lglasso

Description

Cross validation for lglasso

Usage

cvlglassofull(
  data,
  group = NULL,
  lambda = NULL,
  random = FALSE,
  nlam = 10,
  lam.min.ratio = 0.01,
  K,
  expFix = 1,
  trace = FALSE,
  NN = NN
)

Arguments

data

raw data
group

group variable
lambda

tuning parameter
nlam

number of tuning parameter
lam.min.ratio

ratio of largest lambda vs smallest lambda
K

cv folds
expFix

given parameter
trace

whether show the process
type

model type
cores

parallel computing

Value

list

Function for computing estimates in EM algorithm

Description

Function for computing estimates in EM algorithm

Usage

importanceEstimates(importancesSample, datai, groupi)

Arguments

importancesSample

the random samples
datai

the data for subject i
groupi

specify how datai is grouped

Value

a list of estimated

Function for generating the samples from posteior distribution in EM algorithm

Description

Function for generating the samples from posteior distribution in EM algorithm

Usage

importanceSample(n, datai, wi, alpha, groupi)

Arguments

n

the number of random samples
datai

the data for subject i
wi

the given precision matrix
alpha

the exponential distribution with rate alpha
groupi

specify how datai is grouped

Value

a data frame for samples and their weights

Longitudinal graphical lasso

Description

This is the main function of the package, which identifies the underlying network model from clustered data. Here clustered data include longitudinal data, or spatially correlated data, e.g, metabolites in different tissues of a same subject. . or more broadly, clustered data for given tuning parameters.

Usage

lglasso(
  data,
  lambda,
  group = NULL,
  random = FALSE,
  expFix = 1,
  N = 100,
  maxit = 30,
  tol = 10^(-1),
  lower = c(0.01, 0.1),
  upper = c(10, 5),
  start = c("cold", "warm"),
  w.init = NULL,
  wi.init = NULL,
  trace = FALSE,
  type = c("expFixed"),
  ...
)

Arguments

data

a n by (p+2) data frame in which the first column is subject ID, the second column is the time point for longitudinal data or tissue ID.
lambda

vector of length 1 or 2, which is the tuning parameter for the identification of the networks. For details, see the explanations in the below.
group

vector of length n if supplied which specify which data points need to be grouped together to infer the heterogeneous networks for, e.g, pre/post vaccination.
expFix

numeric number used in the model specification
maxit

the maximum iterations for the estimation.
tol

the minimum value for convergence criterion
lower

vector of length 1 or 2 which specifies the lower bounds for alpha_1 (and alpha_2) in the correlation matrix
upper

vector of length 1 or 2 which specifies the upper bounds for alpha_1 (and alpha_2) in the correlation matrix
start

how to start the initial values for lglasso algorithm
w.init

initial value for covariance matrix
wi.init

inital value for precision matrix
trace

whether or not show the progress of the computation
type

a string specifying which model need to be fitted. There are three models available, which are referred as general, expFixed respectively. Please see the details in the below for the meaning of each.
...

other inputs

Details

This function implements three statistical models for network inference, according to how the correlations is specified between time points (or tissues or contents in some clinical studies). These three models are referred as general, expFixed.Let's say we have two time points,t_i,t_j, then in model general, the correlation is tau_ij, while in model expFixed, we have tau=exp(-alpha_1|t_1-t_2|^(-alpha_2)) with alpha_2 need to be pre-specified (default is alpha_2=1). In model twoPara, both alpha_1 and alpha_2 is unknown and need to be inferred from the data. For longitudinal data, model expFixed is recommended while for omics data from different tissues or contents, model general should be adopted.

Value

list which include following components:

w the general covariance matrix estimate;

wList list representing the individual covariance matrix estimate;

wi the general precision matrix estimate;

wiList list representing the individual precision matrix estimate;

v the correlation matrix between specified classes;

vList list representing the individual correlation matrix;

tauhat the correlation parameters for longitudinal data

Title

Description

Title

Usage

lglassoHeter(
  data,
  lambda,
  expFix,
  group,
  maxit,
  tol = 10^(-3),
  trace = FALSE,
  start = c("warm", "cold"),
  w.init = NULL,
  wi.init = NULL,
  N
)

Arguments

trace

Title

Description

Title

Usage

phifunction(t, tau)

Arguments

t

a vector specify the time points corresponding to the data
tau

the damping rate parameter with length 1 or 2

Value

a square matrix used to constrct the likelihood

Plot function for CVlglasso

Description

Plot function for CVlglasso

Usage

## S3 method for class 'cvlglasso'
plot(x, xvar = c("lambda", "step"), ...)

Arguments

x

CVlglasso object
xvar

character which specify the x axis of the plot
...

other plot arguments

Value

If group is NULL in CVlglasso, then a line plot will produced; otherwise, a heatmap will be produced.