Identify networks from longitudinal data (1)

# number of nodes
p=200
# number of edge in general network
m1=200
# the difference between the number of edges in individual networks and general network
m2=10
# number of subjects
n=100
## number of time points
tt=10
set.seed(1)

One-stage model

Estimate the network based on homogeneous one-stage model

data

dd=Simulate(type="homo",n=n,p=p,m1=m1,m2=m2,tau=2,tt=tt,alpha=1,group=2)
ddata=dd$data
dim(ddata)
#> [1] 570 202
ddata[1:2,1:5]
#>    subject     time        V1        V2        V3
#> 1 subject1 1.181643 0.6543404 -0.329604 0.5008536
#> 2 subject1 1.327350 0.9737273 -0.425873 0.7668399

Estimation

aa=lglasso(data=ddata,lambda = 0.01)
system.time(lglasso(data=ddata,lambda = 0.01))
#>    user  system elapsed 
#>   4.275   0.382   4.308
## estimated network 
estimates=lapply(aa$wi,function(ll){ifelse(abs(ll)>10^(-5),1,0)})
## estimated network
estimates[[1]][1:10,1:10]
#>     V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
#> V1   1  0  0  0  0  0  1  1  1   0
#> V2   0  1  0  0  1  0  1  0  1   1
#> V3   0  0  1  1  0  0  1  0  1   0
#> V4   0  0  1  1  1  0  0  1  1   0
#> V5   0  1  0  1  1  0  0  1  0   0
#> V6   0  0  0  0  0  1  0  0  0   1
#> V7   1  1  1  0  0  0  1  0  0   1
#> V8   1  0  0  1  1  0  0  1  0   0
#> V9   1  1  1  1  0  0  0  0  1   1
#> V10  0  1  0  0  0  1  1  0  1   1
## true network 
dd$network[1:10,1:10]
#>     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  0  0  0  0   0
#> 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   0
#> V10  0  0  0  0  0  0  0  0  0   1

## correlation parameter
aa$tau
#> [1] 1.702791
## likelihood
aa$ll
#> [1] 48.80352

Estimate the network based on heterogeneous one-stage model

data

dd=Simulate(type="heter",n=n,p=p,m1=m1,m2=m2,tt=tt,alpha=5,group = 1)
ddata=dd$data$pre
dim(ddata)
#> [1] 562 202
ddata[1:2,1:5]
#>    subject      time        V1        V2         V3
#> 1 subject1 0.9166829 0.9866552 -1.480360 -0.5925246
#> 2 subject1 1.4237719 0.9700890 -1.331953 -0.4729975

Estimation

aa=lglasso(data=ddata,lambda = 0.01,random=TRUE,N=100)
system.time(lglasso(data=ddata,lambda = 0.01,random=TRUE,N=100))
#>    user  system elapsed 
#>   1.658   0.179   1.586
estimates=lapply(aa$wi,function(ll){ifelse(abs(ll)>10^(-5),1,0)})
## estimated network
estimates[[1]][1:10,1:10]
#>     V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
#> V1   1  1  0  1  1  1  0  1  1   1
#> V2   1  1  1  1  1  1  1  1  1   1
#> V3   0  1  1  0  1  1  1  1  1   1
#> V4   1  1  0  1  0  1  1  0  1   1
#> V5   1  1  1  0  1  1  1  0  0   1
#> V6   1  1  1  1  1  1  1  1  0   0
#> V7   0  1  1  1  1  1  1  1  1   0
#> V8   1  1  1  0  0  1  1  1  0   1
#> V9   1  1  1  1  0  0  1  0  1   1
#> V10  1  1  1  1  1  0  0  1  1   1
##  true network
dd$network$pre[1:10,1:10]
#>       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#>  [1,]    1    0    0    0    0    0    0    0    0     0
#>  [2,]    0    1    0    0    0    0    0    0    0     0
#>  [3,]    0    0    1    0    0    0    0    0    0     0
#>  [4,]    0    0    0    1    0    0    0    0    0     0
#>  [5,]    0    0    0    0    1    0    0    0    0     0
#>  [6,]    0    0    0    0    0    1    0    0    0     0
#>  [7,]    0    0    0    0    0    0    1    0    0     0
#>  [8,]    0    0    0    0    0    0    0    1    0     0
#>  [9,]    0    0    0    0    0    0    0    0    1     0
#> [10,]    0    0    0    0    0    0    0    0    0     1
## likelihood
aa$ll
#> [1] 194.4138
## 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)

Two-stage model

Estimate the networks based on homogeneous two-stage model

Generate data

dd=Simulate(type="homo",n=n,p=p,m1=m1,m2=m2,tau=c(2,1),tt=tt)
ddata=do.call(rbind,dd$data)
group=c(rep(0,nrow(ddata)/2),rep(1,nrow(ddata)/2))
dim(ddata)
#> [1] 1198  202
ddata[1:2,1:5]
#>        subject     time         V1        V2        V3
#> pre.1 subject1 1.456898 -0.1694917 -0.194282 0.4388738
#> pre.2 subject1 3.386901 -0.1696782  0.498542 0.4995976

Estimation

aa=lglasso(data=ddata,lambda = c(0.01,0.01),group = group)
system.time(lglasso(data=ddata,lambda = c(0.01,0.01),group = group))
#>    user  system elapsed 
#>  17.388   1.395  17.778
estimates=lapply(aa$wi,function(ll){ifelse(abs(ll)>10^(-5),1,0)})

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

##  true pre-treatment network
dd$network$pre[1:10,1:10]
#>     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   1
#> V3   0  0  1  0  0  0  0  0  0   0
#> 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   0
#> V10  0  1  0  0  0  0  0  0  0   1


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

##  true post-treatment network
dd$network$post[1:10,1:10]
#>     V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
#> V1   0  0  0  0  0  0  0  0  0   0
#> V2   0  0  0  0  0  0  0  0  0   1
#> V3   0  0  0  0  0  0  0  0  0   0
#> V4   0  0  0  0  0  0  0  0  0   0
#> V5   0  0  0  0  0  0  0  0  0   0
#> V6   0  0  0  0  0  0  0  0  0   0
#> V7   0  0  0  0  0  0  0  0  0   0
#> V8   0  0  0  0  0  0  0  0  0   0
#> V9   0  0  0  0  0  0  0  0  0   0
#> V10  0  1  0  0  0  0  0  0  0   0
## correlation parameters
aa$tau
#> [1] 1.8395236 0.8629234

## likelihood
aa$ll
#> [1] 15.28526

Estimate the networks based on hetergeneous two-stage model

Generate data

dd=Simulate(type="heter",n=n,p=p,m1=m1,m2=m2,tt=tt,alpha=0.5,group=2)
ddata=do.call(rbind,dd$data)
group=c(rep(0,nrow(ddata)/2),rep(1,nrow(ddata)/2))
dim(ddata)
#> [1] 1094  202
ddata[1:2,1:5]
#>        subject      time          V1          V2         V3
#> pre.1 subject1 0.7008538 -0.13716301 -0.03287373  0.5416496
#> pre.2 subject1 1.0540184 -0.04239683  0.87577873 -0.6492479

Estimation

aa=lglasso(data=ddata,lambda = c(0.01,0.01),random=TRUE,group = group,N=100)
system.time(lglasso(data=ddata,lambda = c(0.01,0.01),random=TRUE,group = group,N=100))
#>    user  system elapsed 
#>   9.402   0.881   9.407
estimates=lapply(aa$wi,function(ll){ifelse(abs(ll)>10^(-5),1,0)})

## estimated pre-treatment network  
estimates[[1]][1:10,1:10]
#>     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  1  0  0  0  1  1   0
#> V3   0  0  1  1  0  0  0  0  0   0
#> V4   0  1  1  1  1  0  0  1  0   0
#> V5   0  0  0  1  1  0  0  0  1   1
#> V6   0  0  0  0  0  1  0  0  1   0
#> V7   0  0  0  0  0  0  1  1  0   0
#> V8   0  1  0  1  0  0  1  1  0   0
#> V9   0  1  0  0  1  1  0  0  1   1
#> V10  0  0  0  0  1  0  0  0  1   1

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


## estimated post-treatment network  
estimates[[2]][1:10,1:10]
#>     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  1  0  1  1  1   0
#> V3   0  0  1  1  0  0  0  0  0   0
#> V4   0  0  1  1  1  0  0  1  0   1
#> V5   0  1  0  1  1  0  0  0  1   1
#> V6   0  0  0  0  0  1  1  0  1   0
#> V7   0  1  0  0  0  1  1  1  0   1
#> V8   0  1  0  1  0  0  1  1  0   1
#> V9   0  1  0  0  1  1  0  0  1   0
#> V10  0  0  0  1  1  0  1  1  0   1

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

## likelihood
aa$ll
#> [1] 57.69259

## 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)