Identify networks from longitudinal data: homogeneous vs heterogeneous model

# number of nodes
p=10
# number of edge in general network
m1=20
# the difference between the number of edges in individual networks and general network
m2=2
# number of subjects
n=20

One-stage model

Estimate the network based on homogeneous one-stage model

Generate data


set.seed(1)
dd=Simulate(type="longihomo",n=n,p=p,m1=m1,m2=m2,tau=2,tt=10)
ddata=dd$data
subjects=fct_inorder(ddata$subject)
group=rep(0,nrow(ddata))
homoData=list(data=ddata)
saveRDS(homoData,file="homoOneStageData.rds")

Estimation

set.seed(1)
aa=lglasso(data=ddata,lambda = 0.01,trace=T)
#> [1] "iteration 1 precision difference: 2.919 /correlation tau difference: 0.734"
#> [1] "iteration 2 precision difference: 2.036 /correlation tau difference: 1.672"
#> [1] "iteration 3 precision difference: 1.671 /correlation tau difference: 1.55"
#> [1] "iteration 4 precision difference: 1.36 /correlation tau difference: 1.398"
#> [1] "iteration 5 precision difference: 1.318 /correlation tau difference: 1.192"
#> [1] "iteration 6 precision difference: 0.896 /correlation tau difference: 1.153"
#> [1] "iteration 7 precision difference: 0.884 /correlation tau difference: 0.829"
#> [1] "iteration 8 precision difference: 0.466 /correlation tau difference: 0.818"
#> [1] "iteration 9 precision difference: 0.462 /correlation tau difference: 0.443"
#> [1] "iteration 10 precision difference: 0.189 /correlation tau difference: 0.44"
#> [1] "iteration 11 precision difference: 0.188 /correlation tau difference: 0.178"
#> [1] "iteration 12 precision difference: 0.064 /correlation tau difference: 0.177"
#> [1] "iteration 13 precision difference: 0.064 /correlation tau difference: 0.06"
saveRDS(aa,file="homoOneStageResult.rds")
## estimated network 
estimates=lapply(aa$wi,function(ll){ifelse(abs(ll)>10^(-3),1,0)})
estimates
#> [[1]]
#>     V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
#> V1   1  1  1  1  0  0  1  1  0   1
#> V2   1  1  1  0  0  1  1  1  1   1
#> V3   1  1  1  1  0  0  0  1  1   1
#> V4   1  0  1  1  1  1  0  1  1   1
#> V5   0  0  0  1  1  1  1  1  1   1
#> V6   0  1  0  1  1  1  1  0  1   1
#> V7   1  1  0  0  1  1  1  1  1   0
#> V8   1  1  1  1  1  0  1  1  0   1
#> V9   0  1  1  1  1  1  1  0  1   0
#> V10  1  1  1  1  1  1  0  1  0   1
## true network 
dd$network
#>     V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
#> V1   1  1  1  0  0  0  1  0  1   1
#> V2   1  1  1  0  0  1  0  0  0   0
#> V3   1  1  1  1  1  0  0  1  1   1
#> V4   0  0  1  1  1  1  0  0  1   0
#> V5   0  0  1  1  1  1  0  0  0   1
#> V6   0  1  0  1  1  1  0  0  0   1
#> V7   1  0  0  0  0  0  1  1  0   0
#> V8   0  0  1  0  0  0  1  1  0   1
#> V9   1  0  1  1  0  0  0  0  1   0
#> V10  1  0  1  0  1  1  0  1  0   1

## correlation parameter
aa$tau
#> [1] 1.672447

Estimate the network based on heterogeneous one-stage model

Generate data

set.seed(1)
dd=Simulate(type="longiheter",n=n,p=p,m1=m1,m2=m2,tt=10,alpha=5,group=1)
ddata=dd$data$pre
subjects=fct_inorder(ddata$subject)
heterData=list(data=ddata)
saveRDS(heterData,file="heterOneStageData.rds")

Estimation

set.seed(1)
aa=lglasso(data=ddata,lambda = 0.01,random=TRUE,trace=T)
#> [1] "alpha estimate: 1"
#> [1] "iteration 1 precision difference: 1.084 /correlation alpha difference: 8.78"
#> [1] "alpha estimate: 9.78005103277255"
#> [1] "iteration 2 precision difference: 0.174 /correlation alpha difference: 4.705"
#> [1] "alpha estimate: 14.4853318046748"
#> [1] "iteration 3 precision difference: 0.149 /correlation alpha difference: 0.744"
#> [1] "alpha estimate: 13.7413125798929"
#> [1] "iteration 4 precision difference: 0.138 /correlation alpha difference: 1.738"
#> [1] "alpha estimate: 12.0030569326643"
#> [1] "iteration 5 precision difference: 0.109 /correlation alpha difference: 1.302"
#> [1] "alpha estimate: 10.7012686120633"
#> [1] "iteration 6 precision difference: 0.094 /correlation alpha difference: 1.025"
#> [1] "alpha estimate: 9.6763408841541"
#> [1] "iteration 7 precision difference: 0.104 /correlation alpha difference: 0.756"
#> [1] "alpha estimate: 8.92036845357043"
#> [1] "iteration 8 precision difference: 0.084 /correlation alpha difference: 0.504"
#> [1] "alpha estimate: 8.41674562359109"
#> [1] "iteration 9 precision difference: 0.053 /correlation alpha difference: 0.566"
#> [1] "alpha estimate: 7.8503105371897"
#> [1] "iteration 10 precision difference: 0.063 /correlation alpha difference: 0.429"
#> [1] "alpha estimate: 7.42140052806907"
#> [1] "iteration 11 precision difference: 0.099 /correlation alpha difference: 0.252"
#> [1] "alpha estimate: 7.16978192510701"
#> [1] "iteration 12 precision difference: 0.049 /correlation alpha difference: 0.257"
#> [1] "alpha estimate: 6.91265425374284"
#> [1] "iteration 13 precision difference: 0.044 /correlation alpha difference: 0.198"
#> [1] "alpha estimate: 6.71445489940757"
#> [1] "iteration 14 precision difference: 0.042 /correlation alpha difference: 0.076"
saveRDS(aa,file="heterOneStageResult.rds")
## estimated network  

estimates=lapply(aa$wi,function(ll){ifelse(abs(ll)>10^(-3),1,0)})
estimates
#> [[1]]
#>     V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
#> V1   1  1  0  1  1  0  1  1  0   1
#> V2   1  1  1  1  1  0  1  0  1   0
#> V3   0  1  1  0  1  1  0  1  0   0
#> V4   1  1  0  1  0  1  0  1  1   1
#> V5   1  1  1  0  1  1  0  1  1   1
#> V6   0  0  1  1  1  1  1  1  0   1
#> V7   1  1  0  0  0  1  1  0  1   1
#> V8   1  0  1  1  1  1  0  1  0   1
#> V9   0  1  0  1  1  0  1  0  1   0
#> V10  1  0  0  1  1  1  1  1  0   1
##  true network
dd$network
#> $pre
#>       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#>  [1,]    1    1    1    1    0    0    1    1    1     0
#>  [2,]    1    1    0    0    1    0    0    0    1     0
#>  [3,]    1    0    1    0    0    0    0    1    1     0
#>  [4,]    1    0    0    1    0    1    0    0    0     1
#>  [5,]    0    1    0    0    1    0    1    1    1     0
#>  [6,]    0    0    0    1    0    1    1    1    0     1
#>  [7,]    1    0    0    0    1    1    1    0    0     0
#>  [8,]    1    0    1    0    1    1    0    1    1     1
#>  [9,]    1    1    1    0    1    0    0    1    1     0
#> [10,]    0    0    0    1    0    1    0    1    0     1

## correlation parameter
 model <- lm(dd$tau ~ aa$tau)
plot(aa$tau,dd$tau,xlab = "estimated tau",ylab="true tau")
abline(model, col = "red", lwd = 2)

summary(model)
#> 
#> Call:
#> lm(formula = dd$tau ~ aa$tau)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -0.11611 -0.04781 -0.01320  0.01592  0.40745 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) -0.04517    0.04496  -1.005    0.328    
#> aa$tau       1.78607    0.25642   6.966 1.66e-06 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.1091 on 18 degrees of freedom
#> Multiple R-squared:  0.7294, Adjusted R-squared:  0.7144 
#> F-statistic: 48.52 on 1 and 18 DF,  p-value: 1.66e-06

Two-stage model

Estimate the networks based on homogeneous two-stage model

Generate data

set.seed(1)
dd=Simulate(type="longihomo",n=n,p=p,m1=m1,m2=m2,tau=c(2,1),tt=2)
ddata=do.call(rbind,dd$data)
subjects=fct_inorder(ddata$subject)
group=c(rep(0,nrow(ddata)/2),rep(1,nrow(ddata)/2))
homoTwoStageData=list(data=ddata,group=group)
saveRDS(homoTwoStageData,file="homoTwoStageData.rds")

Estimation

set.seed(1)
aa=lglasso(data=ddata,lambda = c(0.1,0.1),group = group,trace=T)
#> [1] "iteration 1 precision difference: 0.865 /correlation tau difference: 0.767"
#> [1] "iteration 2 precision difference: 0.713 /correlation tau difference: 1.322"
#> [1] "iteration 3 precision difference: 0.713 /correlation tau difference: 1.098"
#> [1] "iteration 4 precision difference: 0.442 /correlation tau difference: 1.031"
#> [1] "iteration 5 precision difference: 0.417 /correlation tau difference: 0.592"
#> [1] "iteration 6 precision difference: 0.199 /correlation tau difference: 0.581"
#> [1] "iteration 7 precision difference: 0.195 /correlation tau difference: 0.178"
#> [1] "iteration 8 precision difference: 0.048 /correlation tau difference: 0.176"
#> [1] "iteration 9 precision difference: 0.047 /correlation tau difference: 0.036"
saveRDS(aa,file="homoTwoStageResult.rds")
estimates=lapply(aa$wi,function(ll){ifelse(abs(ll)>10^(-3),1,0)})

## pre-treatment network estimate 
estimates[[1]]
#>     V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
#> V1   1  1  0  0  0  0  0  0  0   0
#> V2   1  1  0  0  0  0  0  0  0   0
#> V3   0  0  1  0  0  0  0  0  0   1
#> V4   0  0  0  1  0  0  0  0  0   0
#> V5   0  0  0  0  1  0  0  0  0   0
#> V6   0  0  0  0  0  1  0  0  0   0
#> V7   0  0  0  0  0  0  1  0  0   0
#> V8   0  0  0  0  0  0  0  1  0   0
#> V9   0  0  0  0  0  0  0  0  1   1
#> V10  0  0  1  0  0  0  0  0  1   1

##  true post-treatment network
dd$network$pre
#>     V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
#> V1   1  1  0  0  1  0  0  1  1   0
#> V2   1  1  1  0  0  0  0  0  0   0
#> V3   0  1  1  1  0  0  0  1  0   1
#> V4   0  0  1  1  0  1  1  1  0   0
#> V5   1  0  0  0  1  1  1  1  0   0
#> V6   0  0  0  1  1  1  1  1  1   0
#> V7   0  0  0  1  1  1  1  1  1   0
#> V8   1  0  1  1  1  1  1  1  0   0
#> V9   1  0  0  0  0  1  1  0  1   1
#> V10  0  0  1  0  0  0  0  0  1   1


## post-treatment network estimate 
estimates[[2]]
#>     V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
#> V1   1  1  0  0  0  0  0  0  0   0
#> V2   1  1  0  0  0  0  0  0  0   0
#> V3   0  0  1  0  0  0  0  0  0   1
#> V4   0  0  0  1  0  0  0  0  0   0
#> V5   0  0  0  0  1  0  0  0  0   0
#> V6   0  0  0  0  0  1  0  0  0   0
#> V7   0  0  0  0  0  0  1  0  0   0
#> V8   0  0  0  0  0  0  0  1  0   0
#> V9   0  0  0  0  0  0  0  0  1   1
#> V10  0  0  1  0  0  0  0  0  1   1

##  true post-treatment network
dd$network$post
#>     V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
#> V1   0  1  0  0  1  0  1  1  1   0
#> V2   1  0  1  0  0  0  0  0  0   0
#> V3   0  1  0  1  0  0  0  1  0   1
#> V4   0  0  1  0  0  1  0  1  0   0
#> V5   1  0  0  0  0  1  1  1  0   0
#> V6   0  0  0  1  1  0  1  1  1   0
#> V7   1  0  0  0  1  1  0  1  1   0
#> V8   1  0  1  1  1  1  1  0  0   0
#> V9   1  0  0  0  0  1  1  0  0   1
#> V10  0  0  1  0  0  0  0  0  1   0
## correlation parameters
aa$tau
#> [1] 1.373564 1.365533

Estimate the networks based on hetergeneous two-stage model

Generate data

set.seed(1)
dd=Simulate(type="longiheter",n=n,p=p,m1=m1,m2=m2,tt=20,alpha=0.5,group=2)
ddata=do.call(rbind,dd$data)
subjects=fct_inorder(ddata$subject)
group=c(rep(0,nrow(ddata)/2),rep(1,nrow(ddata)/2))
heterTwoStageData=list(data=ddata,group=group)
saveRDS(heterTwoStageData,file="heterTwoStageData.rds")

Estimation

set.seed(1)
aa=lglasso(data=ddata,lambda = c(0.1,0.1),random=T,group = group,trace=T)
#> [1] "alpha estimate: 1"
#> [1] "iteration 1 precision difference: 0.588 /correlation alpha difference: 0.826"
#> [1] "alpha estimate: 1.82577740354666"
#> [1] "iteration 2 precision difference: 0.116 /correlation alpha difference: 0.556"
#> [1] "alpha estimate: 1.27024216921394"
#> [1] "iteration 3 precision difference: 0.205 /correlation alpha difference: 0.436"
#> [1] "alpha estimate: 0.834246032526682"
#> [1] "iteration 4 precision difference: 0.14 /correlation alpha difference: 0.223"
#> [1] "alpha estimate: 0.611383893937042"
#> [1] "iteration 5 precision difference: 0.047 /correlation alpha difference: 0.116"
#> [1] "alpha estimate: 0.495392831565826"
#> [1] "iteration 6 precision difference: 0.026 /correlation alpha difference: 0.055"
saveRDS(aa,file="heterTwoStageResult.rds")
estimates=lapply(aa$wi,function(ll){ifelse(abs(ll)>10^(-3),1,0)})

## pre-treatment network estimate 
estimates[[1]]
#>     V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
#> V1   1  0  0  0  0  0  0  0  0   0
#> V2   0  1  0  0  0  0  0  0  0   0
#> V3   0  0  1  0  0  1  0  0  0   0
#> V4   0  0  0  1  0  0  0  0  0   0
#> V5   0  0  0  0  1  0  1  0  0   0
#> V6   0  0  1  0  0  1  0  0  0   0
#> V7   0  0  0  0  1  0  1  0  0   0
#> V8   0  0  0  0  0  0  0  1  0   0
#> V9   0  0  0  0  0  0  0  0  1   0
#> V10  0  0  0  0  0  0  0  0  0   1

##  true pre-treatment network
dd$network$pre
#>       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#>  [1,]    1    1    1    0    0    0    0    0    0     0
#>  [2,]    1    1    1    1    0    1    1    1    1     0
#>  [3,]    1    1    1    1    0    1    0    0    0     1
#>  [4,]    0    1    1    1    1    0    1    0    1     0
#>  [5,]    0    0    0    1    1    0    1    1    0     0
#>  [6,]    0    1    1    0    0    1    0    1    0     0
#>  [7,]    0    1    0    1    1    0    1    1    0     1
#>  [8,]    0    1    0    0    1    1    1    1    0     1
#>  [9,]    0    1    0    1    0    0    0    0    1     0
#> [10,]    0    0    1    0    0    0    1    1    0     1


## post-treatment network estimate 
estimates[[2]]
#>     V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
#> V1   1  0  0  0  0  0  0  0  0   0
#> V2   0  1  0  0  0  0  0  0  0   0
#> V3   0  0  1  0  0  1  0  0  0   0
#> V4   0  0  0  1  0  0  0  0  0   0
#> V5   0  0  0  0  1  0  1  0  0   0
#> V6   0  0  1  0  0  1  0  0  0   0
#> V7   0  0  0  0  1  0  1  0  0   0
#> V8   0  0  0  0  0  0  0  1  0   0
#> V9   0  0  0  0  0  0  0  0  1   0
#> V10  0  0  0  0  0  0  0  0  0   1

##  true post-treatment network
dd$network$post
#>       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#>  [1,]    0    1    1    1    0    0    0    0    0     0
#>  [2,]    1    0    1    1    0    1    1    1    1     0
#>  [3,]    1    1    0    1    0    1    0    0    0     1
#>  [4,]    1    1    1    0    1    0    1    0    1     0
#>  [5,]    0    0    0    1    0    0    1    0    0     0
#>  [6,]    0    1    1    0    0    0    0    1    0     0
#>  [7,]    0    1    0    1    1    0    0    1    0     1
#>  [8,]    0    1    0    0    0    1    1    0    0     1
#>  [9,]    0    1    0    1    0    0    0    0    0     0
#> [10,]    0    0    1    0    0    0    1    1    0     0

## correlation parameters
 model <- lm(dd$tau ~ aa$tau)
plot(aa$tau,dd$tau,xlab = "estimated tau",ylab="true tau")
abline(model, col = "red", lwd = 2)

summary(model)$coef
#>                Estimate Std. Error    t value     Pr(>|t|)
#> (Intercept) -0.09259641  0.4002578 -0.2313419 8.196586e-01
#> aa$tau       1.00005372  0.1408539  7.0999381 1.283921e-06

[1] Zhou J, Gui J, Viles WD, Chen H, Li S, Madan JC, Coker MO, Hoen AG. Identifying stationary microbial interaction networks based on irregularly spaced longitudinal 16S rRNA gene sequencing data. Front Microbiomes. 2024;3:1366948. doi: 10.3389/frmbi.2024.1366948. Epub 2024 Jun 2. PMID: 40687607; PMCID: PMC12276884.

[2] Friedman J., Hastie T., Tibshirani R. (2019) Graphical Lasso: Estimation of Gaussian Graphical Models, Version: 1.11. https://CRAN.R-project.org/package=glasso.