.Rhistory

return(t(mat)); index = attempts}
}
}
updateMessage <- function(message) {
# Clear the previous message by printing '\r' (carriage return)
cat("\r")
# Print the updated message
cat(message)
# Flush the output immediately to ensure the message is displayed
flush.console()
}
# Function:
# nvar = Number of Variables
# auto = value of the AR coefficients
# cross = mean value of the cross-regressions
# sd.cross = SD of the cross-regressions
# bxy = The paths to be filled with non-zero cross-regressions
# attempts = The number of attempts to generate a VAR model before giving up
creat_pop_coef <- function(nvar, auto, cross, sd.cross, bxy, attempts = 50) {
index = 1
while(index <= attempts){
# Creating Dynamics Matrix
mat <- matrix(0, ncol = nvar, nrow = nvar)
diag(mat) <- auto
for (i in 1:nrow(bxy)){
mat[bxy[i,1], bxy[i,2]]  <- rnorm(1, cross,sd.cross)
}
# Stationarity Test
evB1 <- eigen(matrix(mat, nvar, nvar, byrow = T))$values
if(any(sqrt(abs(Re(evB1))^2 + abs(Im(evB1))^2) > 1)){
index = index+1;
if(index == attempts){updateMessage(paste0('Could not find a stationary VAR model with your desired settings'))
break}else{
updateMessage(paste0('Non-Stationary. Retrying. Attempt #',index,' of ', attempts))}
}else{
return(t(mat)); index = attempts}
}
}
# The actual code
set.seed(321993)
# Number of Variables
nvar = 10
# Number of Cross-Regressions desired
paths = 20
# Listing all possible matrix elements
cfs = matrix(c(rep(1:nvar, each=nvar), rep(1:nvar, nvar)), ncol=2, nrow=nvar^2)
# Eliminating AR's from the list of options for possible candidate paths
cfs = cfs[-which(cfs[,1] - cfs[,2] == 0),]
# Simulating a population VAR(1) model with cross-regressions set to rnorm(1, mean = 0.2, sd = 0.05)
mat = creat_pop_coef(nvar = nvar,
auto = 0.7,
cross = -0.2, sd.cross = 0.05,
bxy = cfs[sample(nrow(cfs),paths),])
# Printing the matrix
round(mat, 2)
# Showcasing a failed example:
mat = creat_pop_coef(nvar = nvar,
auto = 0.7,
cross = 0.5, sd.cross = 0.05,
bxy = cfs[sample(nrow(cfs),27),])
# Example usage in a for loop
for (i in 1:10) {
# Construct the message
message <- paste("Processing iteration", i)
# Update the message
updateMessage(paste0('Non-Stationary. Retrying.'))
# Perform some computation or task
# ...
# Add a slight delay to simulate work being done
Sys.sleep(1)
}
# Showcasing a failed example:
mat = creat_pop_coef(nvar = nvar,
auto = 0.7,
cross = 0.5, sd.cross = 0.05,
bxy = cfs[sample(nrow(cfs),27),])
# Showcasing a failed example:
mat = creat_pop_coef(nvar = nvar,
auto = 0.7,
cross = 0.5, sd.cross = 0.05,
bxy = cfs[sample(nrow(cfs),27),])
# Showcasing a failed example:
mat = creat_pop_coef(nvar = nvar,
auto = 0.7,
cross = 0.5, sd.cross = 0.05,
bxy = cfs[sample(nrow(cfs),27),])
# Showcasing a failed example with 27 positive Cross-Regressions with a mean-value of 0.3
mat = creat_pop_coef(nvar = nvar,
auto = 0.7,
cross = 0.3, sd.cross = 0.05,
bxy = cfs[sample(nrow(cfs),27),])
# The actual code
set.seed(321993)
# Number of Variables
nvar = 100
# Number of Cross-Regressions desired
paths = 200
# Listing all possible matrix elements
cfs = matrix(c(rep(1:nvar, each=nvar), rep(1:nvar, nvar)), ncol=2, nrow=nvar^2)
# Eliminating AR's from the list of options for possible candidate paths
cfs = cfs[-which(cfs[,1] - cfs[,2] == 0),]
# Simulating a population VAR(1) model with cross-regressions set to rnorm(1, mean = 0.2, sd = 0.05)
mat = creat_pop_coef(nvar = nvar,
auto = 0.7,
cross = -0.2, sd.cross = 0.05,
bxy = cfs[sample(nrow(cfs),paths),])
mat
sum(abs(mat)>0)
library(dynr); library(qgraph)
library(igraph); library(fclust)
# Sample Characteristics
nv = 5; times = 300; N = 8
set.seed(312489)
# Variables, times, sample size
# Params for Group 1
beta1 = matrix(c(0.62, 0.00, 0.00, 0.00, 0.00,
0.00, 0.53, 0.00, 0.00, 0.00,
0.00, 0.00, 0.61, 0.00, 0.00,
0.00, 0.00, 0.00, 0.52,-0.22,
0.32,-0.16, 0.00, 0.00, 0.52), nv, nv, byrow = TRUE)
psi1 = matrix(c(0.83,0.00,0.00,-0.40,0.40,
0.00,0.99,0.00,0.00,0.00,
0.00,0.00,0.88,0.00,0.00,
-0.40,0.00,0.00,1.01,0.00,
0.40,0.00,0.00,0.00,0.91), nv, nv, byrow = TRUE)
# Params for Group 2
beta2 = matrix(c(0.62,0.00, 0.00, 0.00,0.00,
0.00,0.53, 0.00, 0.00,0.00,
0.00,0.00, 0.61, 0.00,0.00,
0.00,0.00, 0.30, 0.52,0.00,
0.00,-0.16, 0.00, -0.35,0.52), nv, nv, byrow = TRUE)
psi2 = matrix(c(0.83,0.00,0.00,0.00,0.00,
0.00,0.99,0.40,0.00,0.00,
0.00,0.40,0.88,0.40,0.00,
0.00,0.00,0.40,1.01,0.00,
0.00,0.00,0.00,0.00,0.91), nv, nv, byrow = TRUE)
# Variables, times, sample size
# Params for Group 1
beta1 = matrix(c(0.62, 0.00, 0.00, 0.00, 0.00,
0.00, 0.53, 0.00, 0.00, 0.00,
0.00, 0.00, 0.61, 0.00, 0.00,
0.00, 0.00, 0.00, 0.52,-0.22,
0.32,-0.16, 0.00, 0.00, 0.52), nv, nv, byrow = TRUE)
psi1 = matrix(c(0.83,0.00,0.00,-0.40,0.40,
0.00,0.99,0.00,0.00,0.00,
0.00,0.00,0.88,0.00,0.00,
-0.40,0.00,0.00,1.01,0.00,
0.40,0.00,0.00,0.00,0.91), nv, nv, byrow = TRUE)
# Params for Group 2
beta2 = matrix(c(0.62,0.00, 0.00, 0.00,0.00,
0.00,0.53, 0.00, 0.00,0.00,
0.00,0.00, 0.61, 0.00,0.00,
0.00,0.00, 0.30, 0.52,0.00,
0.00,-0.16, 0.00, -0.35,0.52), nv, nv, byrow = TRUE)
psi2 = matrix(c(0.83,0.00,0.00,0.00,0.00,
0.00,0.99,0.40,0.00,0.00,
0.00,0.40,0.88,0.40,0.00,
0.00,0.00,0.40,1.01,0.00,
0.00,0.00,0.00,0.00,0.91), nv, nv, byrow = TRUE)
# Variables, times, sample size
# Params for Group 1
beta1 = matrix(c(0.62, 0.00, 0.00, 0.00, 0.00,
0.00, 0.53, 0.00, 0.00, 0.00,
0.00, 0.00, 0.61, 0.00, 0.00,
0.00, 0.00, 0.00, 0.52,-0.22,
0.32,-0.16, 0.00, 0.00, 0.52), nv, nv, byrow = TRUE)
psi1 = matrix(c(0.83,0.00,0.00,-0.40,0.40,
0.00,0.99,0.00,0.00,0.00,
0.00,0.00,0.88,0.00,0.00,
-0.40,0.00,0.00,1.01,0.00,
0.40,0.00,0.00,0.00,0.91), nv, nv, byrow = TRUE)
# Params for Group 2
beta2 = matrix(c(0.62,0.00, 0.00, 0.00,0.00,
0.00,0.53, 0.00, 0.00,0.00,
0.00,0.00, 0.61, 0.00,0.00,
0.00,0.00, 0.30, 0.52,0.00,
0.00,-0.16, 0.00, -0.35,0.52), nv, nv, byrow = TRUE)
psi2 = matrix(c(0.83,0.00,0.00,0.00,0.00,
0.00,0.99,0.40,0.00,0.00,
0.00,0.40,0.88,0.40,0.00,
0.00,0.00,0.40,1.01,0.00,
0.00,0.00,0.00,0.00,0.91), nv, nv, byrow = TRUE)
# Data Generation
# N/2 in Subgroup 1
# N/2 in Subgroup 2
true.dat = NULL
for(j in 1:N){
dat2fit = matrix(NA, times, nv)
dat2fit[1,] = rnorm(nv)
if(j <= N/2){
beta = beta1; kap = psi1
}else{beta = beta2; kap = psi2}
for(i in 2:times){
dat2fit[i,] = beta %*% dat2fit[i-1,] +
MASS::mvrnorm(1, mu = rep(0, 5), Sigma = kap)
}
true.dat = rbind(true.dat, dat2fit)
}
# Adding 2 Fuzzy Subjects
for(i in 1:2){
dat2fit = matrix(NA, times, nv)
dat2fit[1,] = rnorm(nv)
fuzz = rbinom(times, 1, prob = 0.5)
for(i in 2:times){
if(fuzz[i] == 1){
beta = beta1; kap = psi1
}else{beta = beta2; kap = psi2}
dat2fit[i,] = beta %*% dat2fit[i-1,] +
MASS::mvrnorm(1, mu = rep(0, 5), Sigma = kap)
}
true.dat = rbind(true.dat, dat2fit)
}
# Adding Labels
true.dat = data.frame(true.dat)
true.dat$id = rep(1:(N+2), each = times)
true.dat$time = rep(1:times, (N+2))
source('./dynrVAR.R')
varnames = paste0('Y', 1:4)
varnames
nv=4
ini.loadings <- diag(1, nv)
# Preparing Measurement
meas <- prep.measurement(
values.load = ini.loadings,
params.load = matrix(rep("fixed", (nv^2)), ncol = nv), # Flag for customization
state.names = paste0("X", 1:nv),
obs.names = varnames
)
meas
length(varnames)
source('./dynrVAR.R')
true.dat
head(true.dat)
true.dat[,true.dat$id == 1]
true.dat[true.dat$id == 1,]
poop = true.dat[true.dat$id == 1,]
# Fitting Auto-VAR
# COMMENTED OUT TO SAVE YOU TIME
func.test = dynr.VAR(data = poop, nv = nv, ID = 'id',
time = 'time', dir = '~/Desktop/test/')
# Fitting Auto-VAR
# COMMENTED OUT TO SAVE YOU TIME
func.test = dynr.VAR(data = poop, varnames = paste0('X', 1:5), ID = 'id',
time = 'time', dir = '~/Desktop/test/')
# Fitting Auto-VAR
# COMMENTED OUT TO SAVE YOU TIME
func.test = dynr.VAR(data = poop, varnames = paste0('X', 1:5), id = 'id',
time = 'time', dir = '~/Desktop/test/')
dataframe = true.dat
head(dataframe)
id
id='id'
unique(dataframe[,id])
# 1. Create Chained Time-Series
# Take each individual, i, and chain them with the i+1 subject
# Add a blank space between i and i+1
# Knowns
# 1. Dataset
# 2. ID Variable
chain.temp = NULL
chain.temp
i=1
dataframe[,id==i]
dataframe[id==i,]
id==i
dataframe[,id[i]]
id
dataframe[,id]
dataframe[,id]==id
dataframe[,id]==i
dataframe[dataframe[,id]==i,]
# 1. Create Chained Time-Series
# Take each individual, i, and chain them with the i+1 subject
# Add a blank space between i and i+1
# Knowns
# 1. Dataset
# 2. ID Variable
chain.temp = NULL
for(i in unique(dataframe[,id])){
chain.temp = rbind(chain.temp, dataframe[dataframe[,id]==i,])
}
chain.temp
dataframe[dataframe[,id]==i,]
temp1 = dataframe[dataframe[,id]==i,]
temp1 = dataframe[dataframe[,id]==i,]
tail(temp1)
nrow(temp1)
nrow(temp1)+1
temp1 = dataframe[dataframe[,id]==i,varnames]
dataframe[dataframe[,id]==i,varnames]
varnames
varnames = paste0('X', 1:5)
dataframe[dataframe[,id]==i,varnames]
temp1 = dataframe[dataframe[,id]==i,varnames]
rep(NA, length(varnames))
temp1[nrow(temp1)+1,] = rep(NA, length(varnames))
length(unique(dataframe[,id]))
unique(dataframe[,id])[length(unique(dataframe[,id]))]
# 1. Create Chained Time-Series
# Take each individual, i, and chain them with the i+1 subject
# Add a blank space between i and i+1
# Knowns
# 1. Dataset
# 2. ID Variable
chain.temp = NULL
for(i in unique(dataframe[,id])){
temp1 = dataframe[dataframe[,id]==i,varnames]
if(i != unique(dataframe[,id])[length(unique(dataframe[,id]))]){
temp1[nrow(temp1)+1,] = rep(NA, length(varnames))
}
chain.temp = rbind(chain.temp, temp1)
}
tail(chain.temp)
nrow(chain.temp)
chain.temp[299:305,]
1:nrow(chain.temp)
rep(1, nrow(chain.temp))
chain.temp$TIMES = 1:nrow(chain.temp)
chain.temp$ID = rep(1, nrow(chain.temp))
names(chain.temp)
id
times
time
time = 'time'
names(chain.temp) = c(varnames, time, id)
names(chain.temp)
# Calculating the PDC/PCCs
# See Wild, Eichler, et al (2010)
for (p in 1:nv) {
for (q in 1:nv) {
PCC[p, q] <- -((OG.kappa[p, q]) / (sqrt(OG.kappa[p, p] * OG.kappa[q, q])))
PDC[p, q] <- (OG.betas[p, q]) / (sqrt(OG.betas[p, p] * OG.kappa[q, q] + (OG.betas[p, q]^2)))
}
}
source('./dynrVAR.R')
# Sample Characteristics
nv = 5; times = 300; N = 8
set.seed(312489)
# Variables, times, sample size
# Params for Group 1
beta1 = matrix(c(0.62, 0.00, 0.00, 0.00, 0.00,
0.00, 0.53, 0.00, 0.00, 0.00,
0.00, 0.00, 0.61, 0.00, 0.00,
0.00, 0.00, 0.00, 0.52,-0.22,
0.32,-0.16, 0.00, 0.00, 0.52), nv, nv, byrow = TRUE)
psi1 = matrix(c(0.83,0.00,0.00,-0.40,0.40,
0.00,0.99,0.00,0.00,0.00,
0.00,0.00,0.88,0.00,0.00,
-0.40,0.00,0.00,1.01,0.00,
0.40,0.00,0.00,0.00,0.91), nv, nv, byrow = TRUE)
# Params for Group 2
beta2 = matrix(c(0.62,0.00, 0.00, 0.00,0.00,
0.00,0.53, 0.00, 0.00,0.00,
0.00,0.00, 0.61, 0.00,0.00,
0.00,0.00, 0.30, 0.52,0.00,
0.00,-0.16, 0.00, -0.35,0.52), nv, nv, byrow = TRUE)
psi2 = matrix(c(0.83,0.00,0.00,0.00,0.00,
0.00,0.99,0.40,0.00,0.00,
0.00,0.40,0.88,0.40,0.00,
0.00,0.00,0.40,1.01,0.00,
0.00,0.00,0.00,0.00,0.91), nv, nv, byrow = TRUE)
# Data Generation
# N/2 in Subgroup 1
# N/2 in Subgroup 2
true.dat = NULL
for(j in 1:N){
dat2fit = matrix(NA, times, nv)
dat2fit[1,] = rnorm(nv)
if(j <= N/2){
beta = beta1; kap = psi1
}else{beta = beta2; kap = psi2}
for(i in 2:times){
dat2fit[i,] = beta %*% dat2fit[i-1,] +
MASS::mvrnorm(1, mu = rep(0, 5), Sigma = kap)
}
true.dat = rbind(true.dat, dat2fit)
}
# Adding 2 Fuzzy Subjects
for(i in 1:2){
dat2fit = matrix(NA, times, nv)
dat2fit[1,] = rnorm(nv)
fuzz = rbinom(times, 1, prob = 0.5)
for(i in 2:times){
if(fuzz[i] == 1){
beta = beta1; kap = psi1
}else{beta = beta2; kap = psi2}
dat2fit[i,] = beta %*% dat2fit[i-1,] +
MASS::mvrnorm(1, mu = rep(0, 5), Sigma = kap)
}
true.dat = rbind(true.dat, dat2fit)
}
# Adding Labels
true.dat = data.frame(true.dat)
true.dat$id = rep(1:(N+2), each = times)
true.dat$time = rep(1:times, (N+2))
poop = true.dat[true.dat$id == c(1:3),]
poop
?subset
# Fitting Auto-VAR
# COMMENTED OUT TO SAVE YOU TIME
subset(true.dat, id == c(1:3))
# Fitting Auto-VAR
# COMMENTED OUT TO SAVE YOU TIME
subset(true.dat, id == c(1))
subset(true.dat, id == 1 | id == 2 | id == 3)
# Fitting Auto-VAR
# COMMENTED OUT TO SAVE YOU TIME
poop = subset(true.dat, id == 1 | id == 2 | id == 3)
func.test = dynr.VAR(data = poop, varnames = paste0('X', 1:5),
id = 'id', time = 'time', dir = '~/Desktop/test/')
source('./dynrVAR.R')
# Fitting Auto-VAR
# COMMENTED OUT TO SAVE YOU TIME
poop = subset(true.dat, id == 1 | id == 2)# | id == 3)
func.test = dynr.VAR(data = poop, varnames = paste0('X', 1:5),
id = 'id', time = 'time', dir = '~/Desktop/test/',
MLVAR = TRUE)
func.test
func.test[[2]]
func.test[[2]]$PCC
func.test[[1]][[1]]
func.test[[1]][[i]]$PCC
i=1
func.test[[1]][[i]]$PCC
c(func.test[[1]][[i]]$PCC  )
averagenet = NULL
c(func.test[[1]][[i]]$PCC)
temp.avg = NULL
averagenet = NULL
temp.avg = NULL
averagenet = NULL
for(i in 1:2){
temp.avg = rbind(temp.avg, c(func.test[[1]][[i]]$PCC))
}
temp.avg
colMeans(temp.avg)
matrix(colMeans(temp.avg), length(varnames), length(varnames), byrow=TRUE)
func.test[[1]][[i]]$PCC
func.test[[1]][[1]]$PCC
apply(temp.avg, 2, sd)
temp.avg = NULL
averagenet = NULL
temp.avg = NULL
averagenet = NULL
for(i in 1:2){
temp.avg = rbind(temp.avg, c(func.test[[1]][[i]]$PCC))
sdnet = matrix(apply(temp.avg, 2, sd), length(varnames), length(varnames), byrow = TRUE)
averagenet = matrix(colMeans(temp.avg), length(varnames), length(varnames), byrow = TRUE)
}
sdnet
func.test[[1]][[i]]$paste0('PCC')
dynrVAR.out=func.test
c(dynrVAR.out[[1]][[i]])
dynrVAR.out[[1]][[i]][1]
dynrVAR.out[[1]][[i]][[paste0('Betas')]]
length(dynrVAR.out[[1]])
avgnets = list()
avgnets = list()
list(averagenet, sdnet)
list('Average Network' = averagenet, 'SD Network' = sdnet)
avgnets$paste0(pararms)
j='Betas'
avgnets[paste0(params[j])] = list('Average Network' = averagenet, 'SD Network' = sdnet)
params = c('Betas', 'Resid','PDC', 'PCC')
avgnets = list()
avgnets[paste0(params[j])] = list('Average Network' = averagenet, 'SD Network' = sdnet)
paste0(params[j])
j
params
avgnets[paste0(j)] = list('Average Network' = averagenet, 'SD Network' = sdnet)
paste0(j)
avgnets[[paste0(j)]] = list('Average Network' = averagenet, 'SD Network' = sdnet)
avgnets
avgnets = list()
temp.avg = NULL
averagenet = NULL
for(j in params){
for(i in 1:length(dynrVAR.out[[1]])){
thenets = c(dynrVAR.out[[1]][[i]][[j]])
temp.avg = rbind(temp.avg, thenets)
sdnet = matrix(apply(temp.avg, 2, sd), length(varnames), length(varnames), byrow = TRUE)
averagenet = matrix(colMeans(temp.avg), length(varnames), length(varnames), byrow = TRUE)
}
avgnets[[paste0(j)]] = list('Average Network' = averagenet, 'SD Network' = sdnet)
}
avgnets
avgnets = list()
temp.avg = NULL
averagenet = NULL
for(j in params){
for(i in 1:length(dynrVAR.out[[1]])){
thenets = c(dynrVAR.out[[1]][[i]][[j]])
temp.avg = rbind(temp.avg, thenets)
sdnet = matrix(apply(temp.avg, 2, sd), length(varnames), length(varnames), byrow = TRUE)
averagenet = matrix(colMeans(temp.avg), length(varnames), length(varnames), byrow = TRUE)
}
avgnets[j] = list('Average Network' = averagenet, 'SD Network' = sdnet)
}
avgnets
avgnets = list()
temp.avg = NULL
averagenet = NULL
for(j in params){
for(i in 1:length(dynrVAR.out[[1]])){
thenets = c(dynrVAR.out[[1]][[i]][[j]])
temp.avg = rbind(temp.avg, thenets)
sdnet = matrix(apply(temp.avg, 2, sd), length(varnames), length(varnames), byrow = TRUE)
averagenet = matrix(colMeans(temp.avg), length(varnames), length(varnames), byrow = TRUE)
}
avgnets[[j]] = list('Average Network' = averagenet, 'SD Network' = sdnet)
}
avgnets