.Rhistory

}
sz.col <- unique(structzero.df[structzero.df$where ==
"col", -3])
szcolnum <- numeric(length(l$networks))
for (i in 1:length(l$networks)) {
szcolnum[i] <- nrow(sz.col[sz.col$time == i,
])
}
totrow <- sapply(l$networks, function(x) nrow(as.matrix(x)))
totcol <- sapply(l$networks, function(x) ncol(as.matrix(x)))
if (offset == TRUE) {
dimensions <- rbind(totrow, totcol, szrownum,
szcolnum, totrow - szrownum, totcol - szcolnum)
rownames(dimensions) <- c("total number of rows",
"total number of columns", "row-wise structural zeros",
"column-wise structural zeros", "remaining rows",
"remaining columns")
}
else {
dimensions <- rbind(szrownum, szcolnum, totrow -
szrownum, totcol - szcolnum)
rownames(dimensions) <- c("maximum deleted nodes (row)",
"maximum deleted nodes (col)", "remaining rows",
"remaining columns")
}
colnames(dimensions) <- paste0("t=", t.start:t.end)
if (nrow(structzero.df) > 0) {
if (offset == TRUE) {
message("\nNodes affected completely by structural zeros:")
}
else {
message("\nAbsent nodes:")
}
szcopy <- structzero.df
szcopy$time <- szcopy$time - 1 + t.start
print(unique(szcopy))
}
else {
message("\nAll nodes are retained.")
}
message("\nNumber of nodes per time step after adjustment:")
print(dimensions)
}
}
l$nvertices <- sapply(l$networks, function(x) c(nrow(as.matrix(x)),
ncol(as.matrix(x))))
rownames(l$nvertices) <- c("row", "col")
colnames(l$nvertices) <- paste0("t=", t.start:t.end)
l$offsmat <- list()
for (i in 1:l$time.steps) {
mat <- matrix(0, nrow = nrow(as.matrix(l$networks[[i]])),
ncol = ncol(as.matrix(l$networks[[i]])))
rownames(mat) <- rownames(as.matrix(l$networks[[i]]))
colnames(mat) <- colnames(as.matrix(l$networks[[i]]))
l$offsmat[[i]] <- mat
}
if (nrow(structzero.df) > 0) {
for (i in 1:nrow(structzero.df)) {
if (structzero.df$where[i] == "row") {
index <- which(rownames(l$offsmat[[structzero.df$time[i]]]) ==
structzero.df$label[i])
l$offsmat[[structzero.df$time[i]]][index, ] <- 1
}
else {
index <- which(colnames(l$offsmat[[structzero.df$time[i]]]) ==
structzero.df$label[i])
l$offsmat[[structzero.df$time[i]]][, index] <- 1
}
}
}
if (offset == TRUE) {
l$rhs.terms[length(l$rhs.terms) + 1] <- "offset(edgecov(offsmat[[i]]))"
rhs.operators[length(rhs.operators) + 1] <- "+"
}
else {
if (l$auto.adjust == TRUE) {
l$offsmat <- suppressMessages(btergm::handleMissings(l$offsmat,
na = 1, method = "remove"))
for (j in 1:length(l$covnames)) {
l[[l$covnames[j]]] <- btergm::adjust(l[[l$covnames[j]]],
l$offsmat)
}
}
}
if (verbose == TRUE && length(l$covnames) > 1) {
dimensions <- lapply(lapply(l$covnames, function(x) l[[x]]),
function(y) sapply(y, function(z) dim(as.matrix(z))))
rownames(dimensions[[1]]) <- paste(l$lhs.original, c("(row)",
"(col)"))
for (i in 2:length(dimensions)) {
rownames(dimensions[[i]]) <- c(paste(l$covnames[i],
"(row)"), paste(l$covnames[i], "(col)"))
}
dimensions <- do.call(rbind, dimensions)
colnames(dimensions) <- paste0("t=", t.start:t.end)
message("\nDimensions of the network and covariates after adjustment:")
print(dimensions)
}
rhs <- l$rhs.terms[1]
if (length(rhs.operators) > 0) {
for (i in 1:length(rhs.operators)) {
rhs <- paste(rhs, rhs.operators[i], l$rhs.terms[i +
1])
}
}
f <- paste(lhs, tilde, rhs)
l$form <- stats::as.formula(f, env = environment())
if (blockdiag == TRUE) {
if (l$bipartite == TRUE) {
stop(paste("MCMC estimation is currently only supported for one-mode",
"networks. Use the btergm function instead."))
}
l$form <- stats::update.formula(l$form, networks ~ .)
l$form <- paste(deparse(l$form), collapse = "")
l$form <- paste(l$form, "+ offset(edgecov(offsmat))")
l$form <- stats::as.formula(l$form, env = environment())
if (length(l$covnames) > 1) {
for (j in 2:length(l$covnames)) {
l[[l$covnames[j]]] <- as.matrix(Matrix::bdiag(lapply(l[[l$covnames[j]]],
as.matrix)))
}
}
l$offsmat <- as.matrix(Matrix::bdiag(l$offsmat))
bdoffset <- lapply(l$networks, as.matrix)
for (i in 1:length(bdoffset)) {
bdoffset[[i]][, ] <- 1
}
bdoffset <- as.matrix((Matrix::bdiag(bdoffset) - 1) *
-1)
l$offsmat <- l$offsmat + bdoffset
rm(bdoffset)
l$offsmat[l$offsmat > 0] <- 1
if (class(l$networks[[1]]) == "network") {
attrnames <- network::list.vertex.attributes(l$networks[[1]])
attributes <- list()
for (i in 1:length(l$networks)) {
attrib <- list()
for (j in 1:length(attrnames)) {
attrib[[j]] <- network::get.vertex.attribute(l$networks[[i]],
attrnames[j])
}
attributes[[i]] <- attrib
l$networks[[i]] <- as.matrix(l$networks[[i]])
}
l$networks <- network::network(as.matrix(Matrix::bdiag(l$networks)),
directed = l$directed, bipartite = l$bipartite)
for (i in 1:length(attrnames)) {
attrib <- unlist(lapply(attributes, function(x) x[[i]]))
network::set.vertex.attribute(l$networks, attrnames[i],
attrib)
}
}
else {
l$networks <- network::network(as.matrix(Matrix::bdiag(l$networks)),
directed = l$directed, bipartite = l$bipartite)
}
if (verbose == TRUE) {
cat("\n")
}
}
form3 <- paste(deparse(l$form[[3]]), collapse = "")
form3 <- gsub("\\s+", " ", form3)
l$form <- paste(deparse(l$form[[2]]), deparse(l$form[[1]]),
form3)
return(l)
}
#' Check if a matrix is a one-mode matrix
#'
#' Check if a matrix is a one-mode matrix.
#'
#' @param mat A matrix object containing zeros and ones.
#' @return \code{TRUE} if the input matrix \code{mat} represents a one-mode
#'   network and \code{FALSE} otherwise.
#'
#' @noRd
is.mat.onemode <- function(mat) {
if (nrow(mat) != ncol(mat)) {
return(FALSE)
} else if (!is.null(rownames(mat)) && !is.null(colnames(mat))
&& any(rownames(mat) != colnames(mat))) {
return(FALSE)
} else {
return(TRUE)
}
}
#' Check if a matrix represents a directed network
#'
#' Check if a matrix represents a directed network.
#'
#' @param mat A matrix object containing zeros and ones.
#' @return \code{TRUE} if the input matrix \code{mat} represents a directed
#'   network and \code{FALSE} otherwise.
#'
#' @noRd
is.mat.directed <- function(mat) {
if (nrow(mat) != ncol(mat)) {
return(FALSE)
} else if (!is.null(rownames(mat)) && !is.null(colnames(mat))
&& any(rownames(mat) != colnames(mat), na.rm = TRUE)) {
return(FALSE)
} else {
if (any(as.matrix(mat) != t(as.matrix(mat)), na.rm = TRUE)) {
return(TRUE)
} else {
return(FALSE)
}
}
}
ergm.AME(TERGM_1,var1="nodeifactor.gender.2")
ergm.AME(TERGM_1,var1="nodeofactor.gender.2")
ergm.AME(TERGM_1,var1="nodeifactor.gender.2",var2="nodeofactor.gender.2",inter="nodematch.gender")
ergm.AME(TERGM_1,var1="nodeifactor.smoking.2")
ergm.AME(TERGM_1,var1="nodeofactor.smoking.2")
ergm.AME(TERGM_1,var1="nodeifactor.smoking.2",var2="nodeofactor.smoking.2",inter="nodematch.smoking")
ergm.AME(TERGM_1,var1="nodeicov.program")
ergm.AME(TERGM_1,var1="nodeocov.program")
ergm.AME(TERGM_1,var1="nodeocov.program",var2="nodeicov.program",inter="absdiff.program")
ergm.AME(TERGM_1,var1="edgecov.delrecip")
ergm.AME(TERGM_1,var1="gwesp.fixed.0.7")
ergm.AME(TERGM_1,var1="gwidegree.fixed.0.6")
summary(TERGM_1)
###run the model
set.seed(21093) ###set seed so that results replicate
###estimate the model
TERGM_1<-mtergm(vandebunt~edges+         #equivalent to intercept
nodeifactor("gender")+  #receiver effect of gender
nodeofactor("gender")+  #sender effect of gender
nodematch("gender")+    #same gender
nodeifactor("smoking")+ #receiver effect of smoking
nodeofactor("smoking")+  #sender effect of smoking
nodematch("smoking")+    #same smoking behavior
nodeicov("program")+     ##receiver effect for program
nodeocov("program")+     ##sender effect for program
absdiff("program")+       ##absolute difference for program
delrecip+             #delayed reciprocity
gwesp(.7,fixed=T)+         ###geometrically weighted edgewise shared partnership
gwidegree(.6,fixed=T),##geometrically weighted indegree distribution
control=control.ergm(parallel=3, #run on 3 processors to speed up estimation
MCMC.samplesize = 3000, #increase MCMC sample size to improve convergence
MCMLE.maxit=7)
)
TERGM_1<-mtergm(vandebunt~edges+         #equivalent to intercept
nodeifactor("gender")+  #receiver effect of gender
nodeofactor("gender")+  #sender effect of gender
nodematch("gender")+    #same gender
nodeifactor("smoking")+ #receiver effect of smoking
nodeofactor("smoking")+  #sender effect of smoking
nodematch("smoking")+    #same smoking behavior
nodeicov("program")+     ##receiver effect for program
nodeocov("program")+     ##sender effect for program
absdiff("program")+       ##absolute difference for program
delrecip+             #delayed reciprocity
gwesp(.7,fixed=T)+         ###geometrically weighted edgewise shared partnership
gwidegree(.7,fixed=T),##geometrically weighted indegree distribution
control=control.ergm(parallel=3, #run on 3 processors to speed up estimation
MCMC.samplesize = 3000, #increase MCMC sample size to improve convergence
MCMLE.maxit=6)
)
###run the model
set.seed(21093) ###set seed so that results replicate
###estimate the model
TERGM_1<-mtergm(vandebunt~edges+         #equivalent to intercept
nodeifactor("gender")+  #receiver effect of gender
nodeofactor("gender")+  #sender effect of gender
nodematch("gender")+    #same gender
nodeifactor("smoking")+ #receiver effect of smoking
nodeofactor("smoking")+  #sender effect of smoking
nodematch("smoking")+    #same smoking behavior
nodeicov("program")+     ##receiver effect for program
nodeocov("program")+     ##sender effect for program
absdiff("program")+       ##absolute difference for program
delrecip+             #delayed reciprocity
gwesp(.7,fixed=T)+         ###geometrically weighted edgewise shared partnership
gwidegree(.7,fixed=T),##geometrically weighted indegree distribution
control=control.ergm(parallel=3, #run on 3 processors to speed up estimation
MCMC.samplesize = 3000, #increase MCMC sample size to improve convergence
MCMLE.maxit=5)
)
summary(TERGM_1)
library(ergMargins)
ergm.AME(TERGM_1,var1="nodeifactor.gender.2")
ergm.AME(TERGM_1,var1="nodeofactor.gender.2")
ergm.AME(TERGM_1,var1="nodeifactor.gender.2",var2="nodeofactor.gender.2",inter="nodematch.gender")
ergm.AME(TERGM_1,var1="nodeifactor.smoking.2")
ergm.AME(TERGM_1,var1="nodeofactor.smoking.2")
ergm.AME(TERGM_1,var1="nodeifactor.smoking.2",var2="nodeofactor.smoking.2",inter="nodematch.smoking")
ergm.AME(TERGM_1,var1="nodeicov.program")
ergm.AME(TERGM_1,var1="nodeocov.program")
ergm.AME(TERGM_1,var1="nodeocov.program",var2="nodeicov.program",inter="absdiff.program")
ergm.AME(TERGM_1,var1="edgecov.delrecip")
ergm.AME(TERGM_1,var1="gwesp.fixed.0.7")
ergm.AME(TERGM_1,var1="gwidegree.fixed.0.6")
devtools::build()
devtools::check()
vif.ergm<-function(my.ergm){
if(class(my.ergm)%in%"btergm"){
data_mat<-my.ergm@effects
corr5<-stats::cor(data_mat[!rownames(data_mat)%in%"edges",
!colnames(data_mat)%in%"edges"]) ##omit edges term
beta<-btergm:::coef(my.ergm)
}else{
#get correlation matrix
if(class(my.ergm)%in%"mlergm"){
cor.mat<-stats::cov2cor(solve(my.ergm$information_matrix))
beta<-my.ergm$theta
}else{
cor.mat<-stats::cov2cor(my.ergm$covar) #calculate correlation matrix
beta<-ergm:::coef.ergm(my.ergm)
}
#omit edges, assign names to matrix
rownames(cor.mat)<-colnames(cor.mat)<-names(beta)
corr5<-cor.mat[!rownames(cor.mat)%in%"edges",
!colnames(cor.mat)%in%"edges"]
}
corr5<-corr5[!is.na(corr5[1:nrow(corr5)]),]
corr5<-corr5[,which(!is.na(corr5[1,1:ncol(corr5)]))]
VIFS<-matrix(0,nrow=1,ncol=ncol(corr5))
for(i in 1:ncol(corr5)){
gvec<-as.vector(corr5[-c(i),i]) ##create vector of correlations between covariate of interest and other covariates in the model
tgvec<-t(gvec)
xcor<-solve(corr5[-c(i),-c(i)]) ##create square matrix of correlations between covariates in the model other than the one of interest
Rsq<-tgvec%*%xcor%*%gvec
VIFS[1,i]<-1/(1-Rsq)
}
colnames(VIFS)<-names(beta[!names(beta)%in%"edges"])
if(class(my.ergm)%in%"btergm"){
warning("VIFS for bootstrap TERGM based on model matrix, not the covariance matrix of the estimator. Benchmarks used for MCMC ML estimation may not apply.")
}else{
message("Higher values indicate greater correlation.\nVIF > 20 is concerning, VIF > 100 indicates severe collinearity.")
}
VIFS
}
library(btergm)
networks <- list()
for (i in 1:10) {              # create 10 random networks with 10 actors
mat <- matrix(rbinom(100, 1, .25), nrow = 10, ncol = 10)
diag(mat) <- 0               # loops are excluded
nw <- network::network(mat)  # create network object
networks[[i]] <- nw          # add network to the list
}
covariates <- list()
for (i in 1:10) {              # create 10 matrices as covariate
mat <- matrix(rnorm(100), nrow = 10, ncol = 10)
covariates[[i]] <- mat       # add matrix to the list
}
fit <- btergm(networks ~ edges + istar(2) + edgecov(covariates), R = 100)
summary(fit)                   # show estimation results
###VIF for btergm
VIF.ERGM_BS<-function(my.ergm){
require(ergm)
data_mat<-my.ergm@effects
corr5<-cor(data_mat[!rownames(data_mat)%in%"edges",
!colnames(data_mat)%in%"edges"]) ##omit edges term
corr5<-corr5[!is.na(corr5[1:nrow(corr5)]),]
corr5<-corr5[,which(!is.na(corr5[1,1:ncol(corr5)]))]
VIFS<-matrix(0,nr=1,nc=ncol(corr5))
for(i in 1:ncol(corr5)){
gvec<-as.vector(corr5[-c(i),i]) ##create vector of correlations between covariate of interest and other covariates in the model
tgvec<-t(gvec)
xcor<-solve(corr5[-c(i),-c(i)]) ##create square matrix of correlations between covariates in the model other than the one of interest
Rsq<-tgvec%*%xcor%*%gvec
VIFS[1,i]<-1/(1-Rsq)
}
colnames(VIFS)<-names(my.ergm@coef[!names(my.ergm@coef)%in%"edges"])
VIFS
}
##check the results
VIF.ERGM_BS(fit)
model2<-ergm(faux.mesa.high~edges+
nodematch("Race")+
nodecov("Grade"),#+
control = control.ergm(MCMLE.maxit = 2))
data("faux.mesa.high")
model2<-ergm(faux.mesa.high~edges+
nodematch("Race")+
nodecov("Grade"),#+
control = control.ergm(MCMLE.maxit = 2))
vif.ergm(model2)
stats::coef(model2)
model2<-ergm(faux.mesa.high~edges+
nodematch("Race")+
nodecov("Grade")+
gwesp(fixed=T,.7),
control = control.ergm(MCMLE.maxit = 2))
ergm:::coef.ergm(model2)
stats::coef(model2)
stats::coef(fit)
vif.ergm<-function(my.ergm){
if(class(my.ergm)%in%"btergm"){
data_mat<-my.ergm@effects
corr5<-stats::cor(data_mat[!rownames(data_mat)%in%"edges",
!colnames(data_mat)%in%"edges"]) ##omit edges term
beta<-btergm:::coef(my.ergm)
}else{
#get correlation matrix
if(class(my.ergm)%in%"mlergm"){
cor.mat<-stats::cov2cor(solve(my.ergm$information_matrix))
beta<-my.ergm$theta
}else{
cor.mat<-stats::cov2cor(my.ergm$covar) #calculate correlation matrix
beta<-stats::coef.ergm(my.ergm)
}
#omit edges, assign names to matrix
rownames(cor.mat)<-colnames(cor.mat)<-names(beta)
corr5<-cor.mat[!rownames(cor.mat)%in%"edges",
!colnames(cor.mat)%in%"edges"]
}
corr5<-corr5[!is.na(corr5[1:nrow(corr5)]),]
corr5<-corr5[,which(!is.na(corr5[1,1:ncol(corr5)]))]
VIFS<-matrix(0,nrow=1,ncol=ncol(corr5))
for(i in 1:ncol(corr5)){
gvec<-as.vector(corr5[-c(i),i]) ##create vector of correlations between covariate of interest and other covariates in the model
tgvec<-t(gvec)
xcor<-solve(corr5[-c(i),-c(i)]) ##create square matrix of correlations between covariates in the model other than the one of interest
Rsq<-tgvec%*%xcor%*%gvec
VIFS[1,i]<-1/(1-Rsq)
}
colnames(VIFS)<-names(beta[!names(beta)%in%"edges"])
if(class(my.ergm)%in%"btergm"){
warning("VIFS for bootstrap TERGM based on model matrix, not the covariance matrix of the estimator. Benchmarks used for MCMC ML estimation may not apply.")
}else{
message("Higher values indicate greater correlation.\nVIF > 20 is concerning, VIF > 100 indicates severe collinearity.")
}
VIFS
}
vif.ergm(model2)
vif.ergm<-function(my.ergm){
if(class(my.ergm)%in%"btergm"){
data_mat<-my.ergm@effects
corr5<-stats::cor(data_mat[!rownames(data_mat)%in%"edges",
!colnames(data_mat)%in%"edges"]) ##omit edges term
beta<-btergm:::coef(my.ergm)
}else{
#get correlation matrix
if(class(my.ergm)%in%"mlergm"){
cor.mat<-stats::cov2cor(solve(my.ergm$information_matrix))
beta<-my.ergm$theta
}else{
cor.mat<-stats::cov2cor(my.ergm$covar) #calculate correlation matrix
beta<-stats::coef(my.ergm)
}
#omit edges, assign names to matrix
rownames(cor.mat)<-colnames(cor.mat)<-names(beta)
corr5<-cor.mat[!rownames(cor.mat)%in%"edges",
!colnames(cor.mat)%in%"edges"]
}
corr5<-corr5[!is.na(corr5[1:nrow(corr5)]),]
corr5<-corr5[,which(!is.na(corr5[1,1:ncol(corr5)]))]
VIFS<-matrix(0,nrow=1,ncol=ncol(corr5))
for(i in 1:ncol(corr5)){
gvec<-as.vector(corr5[-c(i),i]) ##create vector of correlations between covariate of interest and other covariates in the model
tgvec<-t(gvec)
xcor<-solve(corr5[-c(i),-c(i)]) ##create square matrix of correlations between covariates in the model other than the one of interest
Rsq<-tgvec%*%xcor%*%gvec
VIFS[1,i]<-1/(1-Rsq)
}
colnames(VIFS)<-names(beta[!names(beta)%in%"edges"])
if(class(my.ergm)%in%"btergm"){
warning("VIFS for bootstrap TERGM based on model matrix, not the covariance matrix of the estimator. Benchmarks used for MCMC ML estimation may not apply.")
}else{
message("Higher values indicate greater correlation.\nVIF > 20 is concerning, VIF > 100 indicates severe collinearity.")
}
VIFS
}
vif.ergm(model2)
model2<-ergm(faux.mesa.high~edges+
nodematch("Race")+
nodecov("Grade"),+
# gwesp(fixed=T,.7),
control = control.ergm(MCMLE.maxit = 2))
model2<-ergm(faux.mesa.high~edges+
nodematch("Race")+
nodecov("Grade"),
# gwesp(fixed=T,.7),
control = control.ergm(MCMLE.maxit = 2))
vif.ergm(model2)
model_est <- mlergm(classes ~ edges +
mutual +
nodematch("sex"),# +
#gwesp(fixed = T,.7),
seed = 123,
options = set_options(number_cores = 2))
library(mlergm)
model_est <- mlergm(classes ~ edges +
mutual +
nodematch("sex"),# +
#gwesp(fixed = T,.7),
seed = 123,
options = set_options(number_cores = 2))
vif.ergm(model_est)
model2<-ergm(faux.mesa.high~
nodematch("Race")+
nodecov("Grade"),
# gwesp(fixed=T,.7),
control = control.ergm(MCMLE.maxit = 2))
vif.ergm(model2)
model2<-ergm(faux.mesa.high~edges+
nodematch("Race")+
nodecov("Grade")+
nodematch("Grade"),
# gwesp(fixed=T,.7),
control = control.ergm(MCMLE.maxit = 2))
vif.ergm(model2)
model2<-ergm(faux.mesa.high~
nodematch("Race")+
nodecov("Grade")+
nodematch("Grade"),
# gwesp(fixed=T,.7),
control = control.ergm(MCMLE.maxit = 2))
vif.ergm(model2)
devtools::build()
devtools::check()