## ----setup, include=FALSE-----------------------------------------------------
rmarkdown::find_pandoc(version = '2.9.2.1')
knitr::opts_chunk$set(echo = TRUE)
## ----IV1, echo = TRUE, eval = TRUE--------------------------------------------
library(PartialNetwork)
set.seed(123)
# Number of groups
M <- 30
# size of each group
N <- rep(50,M)
# individual effects
beta <- c(2,1,1.5)
# endogenous effects
alpha <- 0.4
# std-dev errors
se <- 1
# network distribution
distr <- runif(sum(N*(N-1)))
distr <- vec.to.mat(distr, N, normalise = FALSE)
# covariates
X <- cbind(rnorm(sum(N),0,5),rpois(sum(N),7))
# true network
G0 <- sim.network(distr)
# normalise
G0norm <- norm.network(G0)
# simulate dependent variable use an external package
y <- CDatanet::simsar(~ X, Glist = G0norm, theta = c(alpha, beta, se))
y <- y$y
# generate instruments
instr <- sim.IV(distr, X, y, replication = 1, power = 1)
GY1c1 <- instr[[1]]$G1y # proxy for Gy (draw 1)
GXc1 <- instr[[1]]$G1X[,,1] # proxy for GX (draw 1)
GXc2 <- instr[[1]]$G2X[,,1] # proxy for GX (draw 2)
# build dataset
# keep only instrument constructed using a different draw than the one used to proxy Gy
dataset <- as.data.frame(cbind(y, X, GY1c1, GXc1, GXc2))
colnames(dataset) <- c("y","X1","X2","G1y", "G1X1", "G1X2", "G2X1", "G2X2")
## ----IV2, echo = TRUE, eval = TRUE, message=FALSE-----------------------------
library(AER)
# Same draws
out.iv1 <- ivreg(y ~ X1 + X2 + G1y | X1 + X2 + G1X1 + G1X2, data = dataset)
summary(out.iv1)
## ----IV3, echo = TRUE, eval = TRUE, message=FALSE-----------------------------
# Different draws
out.iv2 <- ivreg(y ~ X1 + X2 + G1y | X1 + X2 + G2X1 + G2X2, data = dataset)
summary(out.iv2)
## ----IV3bias------------------------------------------------------------------
ddS <- as.matrix(cbind(1, dataset[,c("X1", "X2", "G1y")])) #\ddot{S}
dZ <- as.matrix(cbind(1, dataset[,c("X1", "X2", "G2X1", "G2X2")]))#\dot{Z}
dZddS <- crossprod(dZ, ddS)/sum(N)
W <- solve(crossprod(dZ)/sum(N))
matM <- solve(crossprod(dZddS, W%*%dZddS), crossprod(dZddS, W))
maxbias <- apply(sapply(1:1000, function(...){
dddGy <- peer.avg(sim.network(distr, normalise = TRUE) , y)
abs(matM%*%crossprod(dZ, dddGy - dataset$G1y)/sum(N))
}), 1, max); names(maxbias) <- c("(Intercept)", "X1", "X2", "G1y")
{cat("Maximal absolute bias\n"); print(maxbias)}
## ----IV4, echo = TRUE, eval = TRUE--------------------------------------------
rm(list = ls())
library(PartialNetwork)
set.seed(123)
# Number of groups
M <- 30
# size of each group
N <- rep(50,M)
# individual effects
beta <- c(2,1,1.5)
# contextual effects
gamma <- c(5, -3)
# endogenous effects
alpha <- 0.4
# std-dev errors
se <- 1
# network distribution
distr <- runif(sum(N*(N-1)))
distr <- vec.to.mat(distr, N, normalise = FALSE)
# covariates
X <- cbind(rnorm(sum(N),0,5),rpois(sum(N),7))
# true network
G0 <- sim.network(distr)
# normalise
G0norm <- norm.network(G0)
# GX
GX <- peer.avg(G0norm, X)
# simulate dependent variable use an external package
y <- CDatanet::simsar(~ X + GX, Glist = G0norm,
theta = c(alpha, beta, gamma, se))
y <- y$y
# generate instruments
# we need power = 2 for models with contextual effetcs
instr <- sim.IV(distr, X, y, replication = 1, power = 2)
GY1c1 <- instr[[1]]$G1y # proxy for Gy (draw 1)
GXc1 <- instr[[1]]$G1X[,,1] # proxy for GX (draw 1)
GXc2 <- instr[[1]]$G2X[,,1] # proxy for GX (draw 2)
GXc2sq <- instr[[1]]$G2X[,,2] # proxy for G^2X (draw 2)
# build dataset
# keep only instrument constructed using a different draw than the one used to proxy Gy
dataset <- as.data.frame(cbind(y, X, GX, GY1c1, GXc1, GXc2, GXc2sq))
colnames(dataset) <- c("y","X1","X2", "GX1", "GX2", "G1y", "G1X1", "G1X2", "G2X1", "G2X2",
"G2X1sq", "G2X2sq")
## ----IV5, echo = TRUE, eval = TRUE, message=FALSE-----------------------------
# Different draws
out.iv2 <- ivreg(y ~ X1 + X2 + GX1 + GX2 + G1X1 + G1X2 + G1y | X1 + X2 + GX1 +
GX2 + G1X1 + G1X2 + G2X1 + G2X2 + G2X1sq + G2X2sq,
data = dataset)
summary(out.iv2)
## ----IV5bias------------------------------------------------------------------
ddS <- as.matrix(cbind(1, dataset[,c("X1", "X2", "GX1", "GX2", "G1X1", "G1X2",
"G1y")]))
dZ <- as.matrix(cbind(1, dataset[,c("X1", "X2", "GX1", "GX2", "G1X1",
"G1X2", "G2X1", "G2X2", "G2X1sq", "G2X2sq")]))
dZddS <- crossprod(dZ, ddS)/sum(N)
W <- solve(crossprod(dZ)/sum(N))
matM <- solve(crossprod(dZddS, W%*%dZddS), crossprod(dZddS, W))
maxbias <- apply(sapply(1:1000, function(...){
dddGy <- peer.avg(sim.network(distr, normalise = TRUE) , y)
abs(matM%*%crossprod(dZ, dddGy - dataset$G1y)/sum(N))
}), 1, max); names(maxbias) <- c("(Intercept)", "X1", "X2", "GX1", "GX2", "G1X1",
"G1X2", "G1y")
{cat("Maximal absolute bias\n"); print(maxbias)}
## ----smm1, echo = TRUE, eval = TRUE-------------------------------------------
rm(list = ls())
library(PartialNetwork)
set.seed(123)
# Number of groups
M <- 100
# size of each group
N <- rep(30,M)
# individual effects
beta <- c(2, 1, 1.5, 5, -3)
# endogenous effects
alpha <- 0.4
# std-dev errors
se <- 1
# network distribution
distr <- runif(sum(N*(N-1)))
distr <- vec.to.mat(distr, N, normalise = FALSE)
# covariates
X <- cbind(rnorm(sum(N),0,5),rpois(sum(N),7))
# true network
G0 <- sim.network(distr)
# normalise
G0norm <- norm.network(G0)
# Matrix GX
GX <- peer.avg(G0norm, X)
# simulate dependent variable use an external package
y <- CDatanet::simsar(~ X + GX, Glist = G0norm, theta = c(alpha, beta, se))
Gy <- y$Gy
y <- y$y
# build dataset
dataset <- as.data.frame(cbind(y, X, Gy, GX))
colnames(dataset) <- c("y","X1","X2", "Gy", "GX1", "GX2")
## ----Smmload, echo = FALSE, eval = TRUE, message=FALSE------------------------
load(url('https://github.com/ahoundetoungan/PartialNetwork/blob/master/datavignettes/smm.rda?raw=true'))
## ----Smm2, echo = TRUE, eval = FALSE, message=FALSE---------------------------
# out.smm1 <- smmSAR(y ~ X1 + X2 | Gy | GX1 + GX2, dnetwork = distr, contextual = T,
# smm.ctr = list(R = 1, print = F), data = dataset)
# summary(out.smm1)
## ----Smm2a, echo = FALSE, eval = TRUE, message=FALSE--------------------------
print(smm$smm1)
## ----Smm3, echo = TRUE, eval = FALSE, message=FALSE---------------------------
# out.smm2 <- smmSAR(y ~ X1 + X2 || GX1 + GX2, dnetwork = distr, contextual = T,
# smm.ctr = list(R = 1, print = F), data = dataset)
# summary(out.smm2)
## ----Smm3a, echo = FALSE, eval = TRUE, message=FALSE--------------------------
print(smm$smm2)
## ----Smm4, echo = TRUE, eval = FALSE, message=FALSE---------------------------
# out.smm3 <- smmSAR(y ~ X1 + X2 | Gy, dnetwork = distr, contextual = T,
# smm.ctr = list(R = 100, print = F), data = dataset)
# summary(out.smm3)
## ----Smm4a, echo = FALSE, eval = TRUE, message=FALSE--------------------------
print(smm$smm3)
## ----Smm5, echo = TRUE, eval = FALSE, message=FALSE---------------------------
# out.smm4 <- smmSAR(y ~ X1 + X2, dnetwork = distr, contextual = T,
# smm.ctr = list(R = 100, print = F), data = dataset)
# summary(out.smm4)
## ----Smm5a, echo = FALSE, eval = TRUE, message=FALSE--------------------------
print(smm$smm4)
## ----smm6a, echo = TRUE, eval = FALSE-----------------------------------------
# rm(list = ls())
## ----smm6b, echo = TRUE, eval = TRUE------------------------------------------
library(PartialNetwork)
set.seed(123)
# Number of groups
M <- 200
# size of each group
N <- rep(30,M)
# individual effects
beta <- c(1, 1, 1.5, 5, -3)
# endogenous effects
alpha <- 0.4
# std-dev errors
se <- 1
# network distribution
distr <- runif(sum(N*(N-1)))
distr <- vec.to.mat(distr, N, normalise = FALSE)
# covariates
X <- cbind(rnorm(sum(N),0,5), rpois(sum(N),7))
# Groups' fixed effects
# In order to have groups' heterogeneity correlated to X (fixed effects),
# We consider the quantile of X2 at 25% in the group
eff <- unlist(lapply(1:M, function(x)
rep(quantile(X[(c(0, cumsum(N))[x]+1):(cumsum(N)[x]),2], probs = 0.25), each = N[x])))
print(c("cor(eff, X1)" = cor(eff, X[,1]), "cor(eff, X2)" = cor(eff, X[,2])))
# We can see that eff is correlated to X2. We can confirm that the correlation is
# strongly significant.
print(c("p.value.cor(eff, X1)" = cor.test(eff, X[,1])$p.value,
"p.value.cor(eff, X2)" = cor.test(eff, X[,2])$p.value))
# true network
G0 <- sim.network(distr)
# normalise
G0norm <- norm.network(G0)
# Matrix GX
GX <- peer.avg(G0norm, X)
# simulate dependent variable use an external package
y <- CDatanet::simsar(~ -1 + eff + X + GX, Glist = G0norm,
theta = c(alpha, beta, se))
Gy <- y$Gy
y <- y$y
# build dataset
dataset <- as.data.frame(cbind(y, X, Gy, GX))
colnames(dataset) <- c("y","X1","X2", "Gy", "GX1", "GX2")
## ----Smm7, echo = TRUE, eval = FALSE, message=FALSE---------------------------
# out.smmeff1 <- smmSAR(y ~ X1 + X2 || GX1 + GX2, dnetwork = distr, contextual = T,
# fixed.effects = T, smm.ctr = list(R = 1, print = F),
# data = dataset)
# summary(out.smmeff1)
## ----Smm7a, echo = FALSE, eval = TRUE, message=FALSE--------------------------
print(smm$smme1)
## ----Smm8, echo = TRUE, eval = FALSE, message=FALSE---------------------------
# out.smmeff2 <- smmSAR(y ~ X1 + X2 | Gy, dnetwork = distr, contextual = T,
# fixed.effects = T, smm.ctr = list(R = 100, print = F),
# data = dataset)
# summary(out.smmeff2)
## ----Smm8a, echo = FALSE, eval = TRUE, message=FALSE--------------------------
print(smm$smme2)
## ----Smm9, echo = TRUE, eval = FALSE, message=FALSE---------------------------
# out.smmeff3 <- smmSAR(y ~ X1 + X2, dnetwork = distr, contextual = T, fixed.effects = T,
# smm.ctr = list(R = 100, print = F), data = dataset)
# summary(out.smmeff3)
## ----Smm9a, echo = FALSE, eval = TRUE, message=FALSE--------------------------
print(smm$smme3)
## ----Smm10, echo = TRUE, eval = FALSE, message=FALSE--------------------------
# fMC <- function(...){
# # Number of groups
# M <- 200
# # size of each group
# N <- rep(30,M)
# # individual effects
# beta <- c(1, 1, 1.5, 5, -3)
# # endogenous effects
# alpha <- 0.4
# # std-dev errors
# se <- 1
# # network distribution
# distr <- runif(sum(N*(N-1)))
# distr <- vec.to.mat(distr, N, normalise = FALSE)
# # covariates
# X <- cbind(rnorm(sum(N),0,5), rpois(sum(N),7))
# # Groups' fixed effects
# # We defined the groups' fixed effect as the quantile at 25% of X2 in the group
# # This implies that the effects are correlated with X
# eff <- unlist(lapply(1:M, function(x)
# rep(quantile(X[(c(0, cumsum(N))[x]+1):(cumsum(N)[x]),2], probs = 0.25), each = N[x])))
# # true network
# G0 <- sim.network(distr)
# # normalise
# G0norm <- norm.network(G0)
# # Matrix GX
# GX <- peer.avg(G0norm, X)
# # simulate dependent variable use an external package
# y <- CDatanet::simsar(~ -1 + eff + X + GX, Glist = G0norm,
# theta = c(alpha, beta, se))
# Gy <- y$Gy
# y <- y$y
# # build dataset
# dataset <- as.data.frame(cbind(y, X, Gy, GX))
# colnames(dataset) <- c("y","X1","X2", "Gy", "GX1", "GX2")
# out.smmeff1 <- smmSAR(y ~ X1 + X2 || GX1 + GX2, dnetwork = distr, contextual = T,
# fixed.effects = T, smm.ctr = list(R = 1, print = F),
# data = dataset)
# out.smmeff2 <- smmSAR(y ~ X1 + X2 | Gy, dnetwork = distr, contextual = T,
# fixed.effects = T, smm.ctr = list(R = 100, print = F),
# data = dataset)
# out.smmeff3 <- smmSAR(y ~ X1 + X2, dnetwork = distr, contextual = T, fixed.effects = T,
# smm.ctr = list(R = 100, print = F), data = dataset)
# out <- data.frame("GX.observed" = out.smmeff1$estimates,
# "Gy.observed" = out.smmeff2$estimates,
# "None.observed" = out.smmeff3$estimates)
# out
# }
## ----Smm10a, echo = TRUE, eval = FALSE, message=FALSE-------------------------
# smm.Monte.C <- lapply(1:250, fMC)
## ----Smm10b, echo = TRUE, eval = FALSE, message=FALSE-------------------------
# Reduce('+', smm.Monte.C)/250
## ----Smm10c, echo = FALSE, eval = TRUE, message=FALSE-------------------------
print(smm$smmmc)
## ----smmsave, echo = FALSE, eval = FALSE--------------------------------------
# smm <- list(smm1 = out.smm1, smm2 = out.smm2,
# smm3 = out.smm3, smm4 = out.smm4,
# smme1 = out.smmeff1, smme2 = out.smmeff2,
# smme3 = out.smmeff3, smmmc = smm$smmmc)
# save(smm, file = "~/Dropbox/Papers - In progress/Partial Network/Package/AH/PartialNetwork/datavignettes/smm.rda")
## ----Smmlogitload, echo = FALSE, eval = TRUE, message=FALSE-------------------
load(url('https://github.com/ahoundetoungan/PartialNetwork/blob/master/datavignettes/smmlogit.rda?raw=true'))
## ----smmp1, echo = TRUE, eval = FALSE-----------------------------------------
# library(PartialNetwork)
# rm(list = ls())
# set.seed(123)
# # Number of groups
# M <- 100
# # size of each group
# N <- rep(30,M)
# # covariates
# X <- cbind(rnorm(sum(N),0,5),rpois(sum(N),7))
# # network formation model parameter
# rho <- c(-0.8, 0.2, -0.1)
# # individual effects
# beta <- c(2, 1, 1.5, 5, -3)
# # endogenous effects
# alpha <- 0.4
# # std-dev errors
# se <- 1
# # network
# tmp <- c(0, cumsum(N))
# X1l <- lapply(1:M, function(x) X[c(tmp[x] + 1):tmp[x+1],1])
# X2l <- lapply(1:M, function(x) X[c(tmp[x] + 1):tmp[x+1],2])
# dist.net <- function(x, y) abs(x - y)
# X1.mat <- lapply(1:M, function(m) {
# matrix(kronecker(X1l[[m]], X1l[[m]], FUN = dist.net), N[m])})
# X2.mat <- lapply(1:M, function(m) {
# matrix(kronecker(X2l[[m]], X2l[[m]], FUN = dist.net), N[m])})
# Xnet <- as.matrix(cbind("Const" = 1,
# "dX1" = mat.to.vec(X1.mat),
# "dX2" = mat.to.vec(X2.mat)))
# ynet <- Xnet %*% rho
# ynet <- c(1*((ynet + rlogis(length(ynet))) > 0))
# G0 <- vec.to.mat(ynet, N, normalise = FALSE)
# # normalise
# G0norm <- norm.network(G0)
# # Matrix GX
# GX <- peer.avg(G0norm, X)
# # simulate dependent variable use an external package
# y <- CDatanet::simsar(~ X + GX, Glist = G0norm, theta = c(alpha, beta, se))
# Gy <- y$Gy
# y <- y$y
# # build dataset
# dataset <- as.data.frame(cbind(y, X, Gy, GX))
# colnames(dataset) <- c("y","X1","X2", "Gy", "GX1", "GX2")
## ----smmp2, echo = TRUE, eval = FALSE-----------------------------------------
# nNet <- nrow(Xnet) # network formation model sample size
# Aobs <- sample(1:nNet, round(0.3*nNet)) # We observed 30%
# # We can estimate rho using the gml function from the stats package
# logestim <- glm(ynet[Aobs] ~ -1 + Xnet[Aobs,], family = binomial(link = "logit"))
# slogestim <- summary(logestim)
# rho.est <- logestim$coefficients
# rho.var <- slogestim$cov.unscaled # we also need the covariance of the estimator
## ----smmp3, echo = TRUE, eval = FALSE-----------------------------------------
# d.logit <- lapply(1:M, function(x) {
# out <- 1/(1 + exp(-rho.est[1] - rho.est[2]*X1.mat[[x]] -
# rho.est[3]*X2.mat[[x]]))
# diag(out) <- 0
# out})
## ----Smmp4, echo = TRUE, eval = FALSE, message=FALSE--------------------------
# smm.logit <- smmSAR(y ~ X1 + X2, dnetwork = d.logit, contextual = T,
# smm.ctr = list(R = 100, print = F), data = dataset)
# summary(smm.logit)
## ----Smmp4b, echo = FALSE, eval = TRUE, message=FALSE-------------------------
print(smmlo$smm1)
## ----Smmp5, echo = TRUE, eval = FALSE, message=FALSE--------------------------
# fdist <- function(rho.est, rho.var, M, X1.mat, X2.mat){
# rho.est1 <- MASS::mvrnorm(mu = rho.est, Sigma = rho.var)
# lapply(1:M, function(x) {
# out <- 1/(1 + exp(-rho.est1[1] - rho.est1[2]*X1.mat[[x]] -
# rho.est1[3]*X2.mat[[x]]))
# diag(out) <- 0
# out})
# }
## ----Smmp6a, echo = TRUE, eval = FALSE, message=FALSE-------------------------
# fdist_args <- list(rho.est = rho.est, rho.var = rho.var, M = M, X1.mat = X1.mat,
# X2.mat = X2.mat)
# summary(smm.logit, dnetwork = d.logit, data = dataset, .fun = fdist, .args = fdist_args,
# sim = 500, ncores = 8) # ncores performs simulations in parallel
## ----Smmp6c, echo = FALSE, eval = TRUE, message=FALSE-------------------------
print(smmlo$smm2)
## ----smmlogitsave, echo = FALSE, eval = FALSE---------------------------------
# smm2 <- summary(smm.logit, dnetwork = d.logit, data = dataset, .fun = fdist, .args = fdist_args, sim = 500, ncores = 8)
# smmlo <- list(smm1 = smm.logit, smm2 = smm2)
# save(smmlo, file = "~/Dropbox/Papers - In progress/Partial Network/Package/AH/PartialNetwork/datavignettes/smmlogit.rda")
## ----BayesNone0, echo=FALSE---------------------------------------------------
load(url('https://github.com/ahoundetoungan/PartialNetwork/blob/master/datavignettes/out.none.rda?raw=true'))
#save(out.none1, out.none2.1, out.none2.2, out.none3.1, out.none3.2, file = "~/Dropbox/Papers - In progress/Partial Network/Package/AH/PartialNetwork/datavignettes/out.none.rda")
## ----BayesNone1as, eval=FALSE-------------------------------------------------
# rm(list = ls())
# library(PartialNetwork)
# set.seed(123)
# # EXAMPLE I: WITHOUT NETWORK FORMATION MODEL
# # Number of groups
# M <- 50
# # size of each group
# N <- rep(30,M)
# # individual effects
# beta <- c(2,1,1.5)
# # contextual effects
# gamma <- c(5,-3)
# # endogenous effects
# alpha <- 0.4
# # std-dev errors
# se <- 1
# # network distribution
# distr <- runif(sum(N*(N-1)))
# distr <- vec.to.mat(distr, N, normalise = FALSE)
# # covariates
# X <- cbind(rnorm(sum(N),0,5),rpois(sum(N),7))
# # true network
# G0 <- sim.network(distr)
# # normalize
# G0norm <- norm.network(G0)
# # GX
# GX <- peer.avg(G0norm, X)
# # simulate dependent variable use an external package
# y <- CDatanet::simsar(~ X + GX, Glist = G0norm,
# theta = c(alpha, beta, gamma, se))
# y <- y$y
# # dataset
# dataset <- as.data.frame(cbind(y, X1 = X[,1], X2 = X[,2]))
## ----BayesNone1ae, eval=FALSE-------------------------------------------------
# # Example I-1: When the network is fully observed
# out.none1 <- mcmcSAR(formula = y ~ X1 + X2, contextual = TRUE, G0.obs = "all",
# G0 = G0, data = dataset, iteration = 2e4)
# summary(out.none1)
## ----BayesNone1b, echo=FALSE--------------------------------------------------
summary(out.none1)
## ----BayesNone21s, eval=FALSE-------------------------------------------------
# # Example I-2: When a part of the network is observed
# # 60% of the network data is observed
# G0.obs <- lapply(N, function(x) matrix(rbinom(x^2, 1, 0.6), x))
## ----BayesNone21e, eval=FALSE-------------------------------------------------
# # replace the non-observed part of the network by 0 (missing links)
# G0.start <- lapply(1:M, function(x) G0[[x]]*G0.obs[[x]])
# # Use network with missing data as the true network
# out.none2.1 <- mcmcSAR(formula = y ~ X1 + X2, contextual = TRUE, G0.obs = "all",
# G0 = G0.start, data = dataset, iteration = 2e4)
# summary(out.none2.1) # the peer effets seem overestimated
## ----BayesNone21b, echo=FALSE-------------------------------------------------
summary(out.none2.1) # the peer effets seem overestimated
## ----BayesNone22a, eval=FALSE-------------------------------------------------
# out.none2.2 <- mcmcSAR(formula = y ~ X1 + X2, contextual = TRUE, G0.obs = G0.obs,
# G0 = G0.start, data = dataset,
# mlinks = list(dnetwork = distr), iteration = 2e4)
# summary(out.none2.2)
## ----BayesNone22b, echo=FALSE-------------------------------------------------
summary(out.none2.2)
## ----BayesNone31as, eval=FALSE------------------------------------------------
# # Example I-3: When only the network distribution is available
# # Simulate a fictitious network and use as true network
# G0.tmp <- sim.network(distr)
# out.none3.1 <- mcmcSAR(formula = y ~ X1 + X2, contextual = TRUE, G0.obs = "all",
# G0 = G0.tmp, data = dataset, iteration = 2e4)
# summary(out.none3.1) # the peer effets seem overestimated
## ----BayesNone31b, echo=FALSE-------------------------------------------------
summary(out.none3.1) # the peer effets seem overestimated
## ----BayesNone32a, eval=FALSE-------------------------------------------------
# out.none3.2 <- mcmcSAR(formula = y ~ X1 + X2, contextual = TRUE, G0.obs = "none",
# data = dataset, mlinks = list(dnetwork = distr), iteration = 2e4)
# summary(out.none3.2)
## ----BayesNone32b, echo=FALSE-------------------------------------------------
summary(out.none3.2)
## ----BayesLog0, echo=FALSE----------------------------------------------------
load(url('https://github.com/ahoundetoungan/PartialNetwork/blob/master/datavignettes/out.logi.rda?raw=true'))
# save(out.logi2.2, out.logi3.2, slogestim, sout.logi4.1, sout.logi4.2, sout.logi4.3, sout.selb1, sout.selb2, file = "~/Dropbox/Papers - In progress/Partial Network/Package/AH/PartialNetwork/datavignettes/out.logi.rda")
## ----BayesLog1as, eval=FALSE--------------------------------------------------
# # EXAMPLE II: NETWORK FORMATION MODEL: LOGIT
# rm(list = ls())
# library(PartialNetwork)
# set.seed(123)
# # Number of groups
# M <- 50
# # size of each group
# N <- rep(30,M)
# # individual effects
# beta <- c(2,1,1.5)
# # contextual effects
# gamma <- c(5,-3)
# # endogenous effects
# alpha <- 0.4
# # std-dev errors
# se <- 2
# # parameters of the network formation model
# rho <- c(-2, -.5, .2)
# # covariates
# X <- cbind(rnorm(sum(N),0,5),rpois(sum(N),7))
# # compute distance between individuals
# tmp <- c(0, cumsum(N))
# X1l <- lapply(1:M, function(x) X[c(tmp[x] + 1):tmp[x+1],1])
# X2l <- lapply(1:M, function(x) X[c(tmp[x] + 1):tmp[x+1],2])
# dist.net <- function(x, y) abs(x - y)
# X1.mat <- lapply(1:M, function(m) {
# matrix(kronecker(X1l[[m]], X1l[[m]], FUN = dist.net), N[m])})
# X2.mat <- lapply(1:M, function(m) {
# matrix(kronecker(X2l[[m]], X2l[[m]], FUN = dist.net), N[m])})
# # true network
# Xnet <- as.matrix(cbind("Const" = 1,
# "dX1" = mat.to.vec(X1.mat),
# "dX2" = mat.to.vec(X2.mat)))
# ynet <- Xnet %*% rho
# ynet <- 1*((ynet + rlogis(length(ynet))) > 0)
# G0 <- vec.to.mat(ynet, N, normalise = FALSE)
# G0norm <- norm.network(G0)
# # GX
# GX <- peer.avg(G0norm, X)
# # simulate dependent variable use an external package
# y <- CDatanet::simsar(~ X + GX, Glist = G0norm,
# theta = c(alpha, beta, gamma, se))
# y <- y$y
# # dataset
# dataset <- as.data.frame(cbind(y, X1 = X[,1], X2 = X[,2]))
## ----BayesLog21as, eval=FALSE-------------------------------------------------
# # Example II-1: When a part of the network is observed
# # 60% of the network data is observed
# G0.obs <- lapply(N, function(x) matrix(rbinom(x^2, 1, 0.6), x))
# # replace the non-observed part of the network by 0
# G0.start <- lapply(1:M, function(x) G0[[x]]*G0.obs[[x]])
# # Infer the missing links in the network data
# mlinks <- list(model = "logit", mlinks.formula = ~ dX1 + dX2,
# mlinks.data = as.data.frame(Xnet))
## ----BayesLog21ae, eval=FALSE-------------------------------------------------
# out.logi2.2 <- mcmcSAR(formula = y ~ X1 + X2, contextual = TRUE, G0.obs = G0.obs,
# G0 = G0.start, data = dataset, mlinks = mlinks,
# iteration = 2e4)
# summary(out.logi2.2)
## ----BayesLog22b, echo=FALSE--------------------------------------------------
summary(out.logi2.2)
## ----BayesLog22c, fig.height = 3, fig.align = "center"------------------------
plot(out.logi2.2, plot.type = "sim", mar = c(3, 2.1, 1, 1))
## ----BayesLog22cc, fig.height = 3, fig.align = "center"-----------------------
plot(out.logi2.2, plot.type = "sim", which.parms = "rho", mar = c(3, 2.1, 1, 1))
## ----BayesLog31as, eval=FALSE-------------------------------------------------
# # Example II-2: When only the network distribution is available
# # Infer the network data
# # We only provide estimate of rho and its variance
# Gvec <- mat.to.vec(G0, ceiled = TRUE)
# logestim <- glm(Gvec ~ -1 + Xnet, family = binomial(link = "logit"))
# slogestim <- summary(logestim)
# estimates <- list("rho" = logestim$coefficients,
# "var.rho" = slogestim$cov.unscaled,
# "N" = N)
# mlinks <- list(model = "logit", mlinks.formula = ~ dX1 + dX2,
# mlinks.data = as.data.frame(Xnet), estimates = estimates)
# out.logi3.2 <- mcmcSAR(formula = y ~ X1 + X2, contextual = TRUE, G0.obs = "none",
# data = dataset, mlinks = mlinks, iteration = 2e4)
# summary(out.logi3.2)
## ----BayesLog32b, echo=FALSE--------------------------------------------------
summary(out.logi3.2)
## ----BayesLog32c, fig.height = 3, fig.align = "center"------------------------
plot(out.logi3.2, plot.type = "sim", mar = c(3, 2.1, 1, 1))
## ----BayesLog32cc, fig.height = 3, fig.align = "center"-----------------------
plot(out.logi3.2, plot.type = "sim", which.parms = "rho", mar = c(3, 2.1, 1, 1))
## ----BayesLog21ae3, eval=FALSE------------------------------------------------
# estimates <- list("rho" = logestim$coefficients,
# "var.rho" = slogestim$cov.unscaled)
# mlinks <- list(model = "logit", mlinks.formula = ~ dX1 + dX2,
# mlinks.data = as.data.frame(Xnet), estimates = estimates)
# out.logi4.1 <- mcmcSAR(formula = y ~ X1 + X2, contextual = TRUE, G0.obs = G0.obs,
# G0 = G0.start, data = dataset, mlinks = mlinks,
# iteration = 2e4)
# summary(out.logi4.1)
## ----BayesLog21ae4, eval=TRUE, echo=FALSE-------------------------------------
print(sout.logi4.1)
## ----BayesLog21ae1, eval=TRUE-------------------------------------------------
print(slogestim)
## ----BayesLog21ae5, eval=FALSE------------------------------------------------
# prior <- list("rho" = c(0, 0, 0),
# "var.rho" = diag(3))
# mlinks <- list(model = "logit", mlinks.formula = ~ dX1 + dX2,
# mlinks.data = as.data.frame(Xnet), prior = prior)
# out.logi4.2 <- mcmcSAR(formula = y ~ X1 + X2, contextual = TRUE, G0.obs = G0.obs,
# G0 = G0.start, data = dataset, mlinks = mlinks,
# iteration = 2e4)
# summary(out.logi4.2)
## ----BayesLog21ae6, echo=FALSE------------------------------------------------
print(sout.logi4.2)
## ----BayesLog21ae8, echo=FALSE------------------------------------------------
print(sout.logi4.3)
## ----Bayesard0, echo=FALSE----------------------------------------------------
load(url('https://github.com/ahoundetoungan/PartialNetwork/blob/master/datavignettes/out.ard.rda?raw=true'))
# save(data.plot1, data.plot2, genz, genv, zi, vk, file = "~/Dropbox/Papers - In progress/Partial Network/Package/AH/PartialNetwork/datavignettes/out.ard.rda")
## ----ard1, eval=FALSE---------------------------------------------------------
# rm(list = ls())
# library(PartialNetwork)
# set.seed(123)
# # LATENT SPACE MODEL
# N <- 500
# genzeta <- 1
# mu <- -1.35
# sigma <- 0.37
# K <- 12 # number of traits
# P <- 3 # Sphere dimension
# # ARD parameters
# # Generate z (spherical coordinates)
# genz <- rvMF(N, rep(0,P))
# # Generate nu from a Normal(mu, sigma^2) (The gregariousness)
# gennu <- rnorm(N, mu, sigma)
# # compute degrees
# gend <- N*exp(gennu)*exp(mu+0.5*sigma^2)*exp(logCpvMF(P,0) - logCpvMF(P,genzeta))
# # Link probabilities
# distr <- sim.dnetwork(gennu, gend, genzeta, genz)
# # Adjacency matrix
# G <- sim.network(distr)
# # Generate vk, the trait location
# genv <- rvMF(K, rep(0, P))
# # set fixed some vk distant
# genv[1,] <- c(1, 0, 0)
# genv[2,] <- c(0, 1, 0)
# genv[3,] <- c(0, 0, 1)
# # eta, the intensity parameter
# geneta <- abs(rnorm(K, 2, 1))
# # Build traits matrix
# densityatz <- matrix(0, N, K)
# for(k in 1:K){
# densityatz[,k] <- dvMF(genz, genv[k,]*geneta[k])
# }
# trait <- matrix(0, N, K)
# NK <- floor(runif(K, 0.8, 0.95)*colSums(densityatz)/apply(densityatz, 2, max))
# for (k in 1:K) {
# trait[,k] <- rbinom(N, 1, NK[k]*densityatz[,k]/sum(densityatz[,k]))
# }
# # Build ADR
# ARD <- G %*% trait
# # generate b
# genb <- numeric(K)
# for(k in 1:K){
# genb[k] <- sum(G[,trait[,k]==1])/sum(G)
# }
## ----ard1as, eval=FALSE-------------------------------------------------------
# # Example1: ARD is observed for the whole population
# # initialization
# d0 <- exp(rnorm(N)); b0 <- exp(rnorm(K)); eta0 <- rep(1,K)
# zeta0 <- 2; z0 <- matrix(rvMF(N, rep(0,P)), N); v0 <- matrix(rvMF(K,rep(0, P)), K)
# # We should fix some vk and bk
# vfixcolumn <- 1:5
# bfixcolumn <- c(3, 7, 9)
# b0[bfixcolumn] <- genb[bfixcolumn]
# v0[vfixcolumn,] <- genv[vfixcolumn,]
# start <- list("z" = z0, "v" = v0, "d" = d0, "b" = b0, "eta" = eta0,
# "zeta" = zeta0)
## ----ard1e, eval=FALSE--------------------------------------------------------
# # MCMC
# estim.ard1 <- mcmcARD(Y = ARD, traitARD = trait, start = start, fixv = vfixcolumn,
# consb = bfixcolumn, iteration = 5000)
## ----ardaa1, eval=FALSE-------------------------------------------------------
# # plot coordinates of individual 123
# i <- 123
# zi <- estim.ard1$simulations$z[i,,]
# par(mfrow = c(3, 1), mar = c(2.1, 2.1, 1, 1))
# invisible(lapply(1:3, function(x) {
# plot(zi[x,], type = "l", ylab = "", col = "blue", ylim = c(-1, 1))
# abline(h = genz[i, x], col = "red")
# }))
## ----ardaa1a, echo=FALSE, fig.height = 4, fig.align = "center"----------------
i <- 123
par(mfrow = c(3, 1), mar = c(2.1, 2.1, 1, 1))
invisible(lapply(1:3, function(x) {
plot(zi[x,], type = "l", ylab = "", col = "blue", ylim = c(-1, 1))
abline(h = genz[i, x], col = "red")
}))
## ----ardaa2, eval=FALSE-------------------------------------------------------
# # plot coordinates of the trait 8
# k <- 8
# vk <- estim.ard1$simulations$v[k,,]
# par(mfrow = c(3, 1), mar = c(2.1, 2.1, 1, 1))
# invisible(lapply(1:3, function(x) {
# plot(vk[x,], type = "l", ylab = "", col = "blue", ylim = c(-1, 1))
# abline(h = genv[k, x], col = "red")
# }))
## ----ardaa2a, echo=FALSE, fig.height = 4, fig.align = "center"----------------
k <- 8
par(mfrow = c(3, 1), mar = c(2.1, 2.1, 1, 1))
invisible(lapply(1:3, function(x) {
plot(vk[x,], type = "l", ylab = "", col = "blue", ylim = c(-1, 1))
abline(h = genv[k, x], col = "red")
}))
## ----ardb, message=FALSE, eval=FALSE------------------------------------------
# # plot degree
# library(ggplot2)
# data.plot1 <- data.frame(True_degree = gend,
# Estim_degree = colMeans(tail(estim.ard1$simulations$d, 2500)))
# ggplot(data = data.plot1, aes(x = True_degree, y = Estim_degree)) +
# geom_abline(col = "red") + geom_point(col = "blue")
## ----ardbb, message=FALSE, echo=FALSE, fig.height = 3, fig.align = "center"----
library(ggplot2)
ggplot(data = data.plot1, aes(x = True_degree, y = Estim_degree)) +
geom_abline(col = "red") + geom_point(col = "blue")
## ----ard2as, eval=FALSE-------------------------------------------------------
# # Example2: ARD is observed for 70% population
# # sample with ARD
# n <- round(0.7*N)
# # individual with ARD
# iselect <- sort(sample(1:N, n, replace = FALSE))
# ARDs <- ARD[iselect,]
# traits <- trait[iselect,]
# # initialization
# d0 <- d0[iselect]; z0 <- z0[iselect,]
# start <- list("z" = z0, "v" = v0, "d" = d0, "b" = b0, "eta" = eta0, "zeta" = zeta0)
## ----ard2ae, eval=FALSE-------------------------------------------------------
# # MCMC
# estim.ard2 <- mcmcARD(Y = ARDs, traitARD = traits, start = start, fixv = vfixcolumn,
# consb = bfixcolumn, iteration = 5000)
## ----ard2b, eval=FALSE--------------------------------------------------------
# # estimation for non ARD
# # we need a logical vector indicating if the i-th element has ARD
# hasARD <- (1:N) %in% iselect
# # we use the matrix of traits to estimate distance between individuals
# estim.nard2 <- fit.dnetwork(estim.ard2, X = trait, obsARD = hasARD, m = 1)
## ----arde, eval=FALSE---------------------------------------------------------
# # estimated degree
# estd <- estim.nard2$degree
# data.plot2 <- data.frame(True_degree = gend,
# Estim_degree = estd,
# Has_ARD = ifelse(hasARD, "YES", "NO"))
# ggplot(data = data.plot2, aes(x = True_degree, y = Estim_degree, colour = Has_ARD)) +
# geom_abline(col = "red") + geom_point()
## ----ardee, message=FALSE, echo=FALSE, fig.height = 3, fig.align = "center"----
ggplot(data = data.plot2, aes(x = True_degree, y = Estim_degree, colour = Has_ARD)) +
geom_abline(col = "red") + geom_point()
## ----Bayesard00, echo=FALSE---------------------------------------------------
load(url('https://github.com/ahoundetoungan/PartialNetwork/blob/master/datavignettes/out.lsp.rda?raw=true'))
#save(out.lspa1, out.lspa2, file = "~/Dropbox/Papers - In progress/Partial Network/Package/AH/PartialNetwork/datavignettes/out.lsp.rda")
## ----bayesard1, eval=FALSE----------------------------------------------------
# rm(list = ls())
# library(PartialNetwork)
# set.seed(123)
# M <- 30
# N <- rep(60, M)
# genzeta <- 3
# mu <- -1.35
# sigma <- 0.37
# K <- 12 # number of traits
# P <- 3 # Sphere dimension
#
# # IN THIS LOOP, WE GENERATE DATA FOLLOWING BREZA ET AL. (2020) AND
# # ESTIMATE THEIR LATENT SPACE MODEL FOR EACH SUB-NETWORK.
# estimates <- list()
# list.trait <- list()
# G0 <- list()
# for (m in 1:M) {
# #######################################################################################
# ####### SIMULATION STAGE ########
# #######################################################################################
# # ARD parameters
# # Generate z (spherical coordinates)
# genz <- rvMF(N[m], rep(0,P))
# # Generate nu from a Normal(mu, sigma^2) (The gregariousness)
# gennu <- rnorm(N[m],mu,sigma)
# # compute degrees
# gend <- N[m]*exp(gennu)*exp(mu+0.5*sigma^2)*exp(logCpvMF(P,0) - logCpvMF(P,genzeta))
# # Link probabilities
# distr <- sim.dnetwork(gennu, gend, genzeta, genz)
# # Adjacency matrix
# G <- sim.network(distr)
# G0[[m]] <- G
# # Generate vk, the trait location
# genv <- rvMF(K, rep(0, P))
# # set fixed some vk distant
# genv[1,] <- c(1, 0, 0)
# genv[2,] <- c(0, 1, 0)
# genv[3,] <- c(0, 0, 1)
# # eta, the intensity parameter
# geneta <-abs(rnorm(K, 2, 1))
# # Build traits matrix
# densityatz <- matrix(0, N[m], K)
# for(k in 1:K){
# densityatz[,k] <- dvMF(genz, genv[k,]*geneta[k])
# }
# trait <- matrix(0, N[m], K)
# NK <- floor(runif(K, .8, .95)*colSums(densityatz)/apply(densityatz, 2, max))
# for (k in 1:K) {
# trait[,k] <- rbinom(N[m], 1, NK[k]*densityatz[,k]/sum(densityatz[,k]))
# }
# list.trait[[m]] <- trait
# # Build ADR
# ARD <- G %*% trait
# # generate b
# genb <- numeric(K)
# for(k in 1:K){
# genb[k] <- sum(G[,trait[,k]==1])/sum(G) + 1e-8
# }
#
# #######################################################################################
# ####### ESTIMATION STAGE ########
# #######################################################################################
# # initialization
# d0 <- gend; b0 <- exp(rnorm(K)); eta0 <- rep(1,K); zeta0 <- genzeta
# z0 <- matrix(rvMF(N[m], rep(0,P)), N[m]); v0 <- matrix(rvMF(K,rep(0, P)), K)
# # We should fix some vk and bk
# vfixcolumn <- 1:5
# bfixcolumn <- c(1, 3, 5, 7, 9, 11)
# b0[bfixcolumn] <- genb[bfixcolumn]
# v0[vfixcolumn,] <- genv[vfixcolumn,]
# start <- list("z" = z0, "v" = v0, "d" = d0, "b" = b0, "eta" = eta0,
# "zeta" = zeta0)
# estimates[[m]] <- mcmcARD(Y = ARD, traitARD = trait, start = start, fixv = vfixcolumn,
# consb = bfixcolumn, sim.d = FALSE, sim.zeta = FALSE,
# iteration = 5000, ctrl.mcmc = list(print = FALSE))
# }
#
# # SIMULATE DATA FOR THE OUTCOME MODEL
# # individual effects
# beta <- c(2,1,1.5)
# # contextual effects
# gamma <- c(5,-3)
# # endogenous effects
# alpha <- 0.4
# # std-dev errors
# se <- 1
# # covariates
# X <- cbind(rnorm(sum(N),0,5),rpois(sum(N),7))
# # Normalise G0
# G0norm <- norm.network(G0)
# # GX
# GX <- peer.avg(G0norm, X)
# # simulate dependent variable use an external package
# y <- CDatanet::simsar(~ X + GX, Glist = G0norm,
# theta = c(alpha, beta, gamma, se))
# y <- y$y
# # dataset
# dataset <- as.data.frame(cbind(y, X1 = X[,1], X2 = X[,2]))
## ----bayesard2s, eval=FALSE---------------------------------------------------
# mlinks <- list(model = "latent space", estimates = estimates)
# out.lspa1 <- mcmcSAR(formula = y ~ X1 + X2, contextual = TRUE, G0.obs = "none",
# data = dataset, mlinks = mlinks, iteration = 2e4)
# summary(out.lspa1)
## ----Bayesard2b, echo=FALSE---------------------------------------------------
summary(out.lspa1)
## ----Bayesard2c, fig.height = 3, fig.align = "center"-------------------------
plot(out.lspa1, plot.type = "sim", mar = c(3, 2.1, 1, 1))
## ----bayesard3, eval=FALSE----------------------------------------------------
# rm(list = ls())
# library(PartialNetwork)
# set.seed(123)
# M <- 30
# N <- rep(60, M)
# genzeta <- 3
# mu <- -1.35
# sigma <- 0.37
# K <- 12
# P <- 3
#
# # IN THIS LOOP, WE GENERATE DATA FOLLOWING BREZA ET AL. (2020) AND
# # ESTIMATE THEIR LATENT SPACE MODEL FOR EACH SUB-NETWORK.
# estimates <- list()
# list.trait <- list()
# obARD <- list()
# G0 <- list()
# for (m in 1:M) {
# #######################################################################################
# ####### SIMULATION STAGE ########
# #######################################################################################
# # ARD parameters
# # Generate z (spherical coordinates)
# genz <- rvMF(N[m], rep(0,P))
# # Generate nu from a Normal(mu, sigma^2) (The gregariousness)
# gennu <- rnorm(N[m],mu,sigma)
# # compute degrees
# gend <- N[m]*exp(gennu)*exp(mu+0.5*sigma^2)*exp(logCpvMF(P,0) - logCpvMF(P,genzeta))
# # Link probabilities
# distr <- sim.dnetwork(gennu, gend, genzeta, genz)
# # Adjacency matrix
# G <- sim.network(distr)
# G0[[m]] <- G
# # Generate vk, the trait location
# genv <- rvMF(K, rep(0, P))
# # set fixed some vk distant
# genv[1,] <- c(1, 0, 0)
# genv[2,] <- c(0, 1, 0)
# genv[3,] <- c(0, 0, 1)
# # eta, the intensity parameter
# geneta <-abs(rnorm(K, 2, 1))
# # Build traits matrix
# densityatz <- matrix(0, N[m], K)
# for(k in 1:K){
# densityatz[,k] <- dvMF(genz, genv[k,]*geneta[k])
# }
# trait <- matrix(0, N[m], K)
# NK <- floor(runif(K, .8, .95)*colSums(densityatz)/apply(densityatz, 2, max))
# for (k in 1:K) {
# trait[,k] <- rbinom(N[m], 1, NK[k]*densityatz[,k]/sum(densityatz[,k]))
# }
# list.trait[[m]] <- trait
# # Build ADR
# ARD <- G %*% trait
# # generate b
# genb <- numeric(K)
# for(k in 1:K){
# genb[k] <- sum(G[,trait[,k]==1])/sum(G) + 1e-8
# }
# # sample with ARD
# n <- round(runif(1, .7, 1)*N[m])
# # individual with ARD
# iselect <- sort(sample(1:N[m], n, replace = FALSE))
# hasARD <- (1:N[m]) %in% iselect
# obARD[[m]] <- hasARD
# ARDs <- ARD[iselect,]
# traits <- trait[iselect,]
# #######################################################################################
# ####### ESTIMATION STAGE ########
# #######################################################################################
# # initialization
# d0 <- gend[iselect]; b0 <- exp(rnorm(K)); eta0 <- rep(1,K); zeta0 <- genzeta
# z0 <- matrix(rvMF(n, rep(0,P)), n); v0 <- matrix(rvMF(K, rep(0, P)), K)
# # We should fix some vk and bk
# vfixcolumn <- 1:5
# bfixcolumn <- c(1, 3, 5, 7, 9, 11)
# b0[bfixcolumn] <- genb[bfixcolumn]; v0[vfixcolumn,] <- genv[vfixcolumn,]
# start <- list("z" = z0, "v" = v0, "d" = d0, "b" = b0, "eta" = eta0,
# "zeta" = zeta0)
# estimates[[m]] <- mcmcARD(Y = ARDs, traitARD = traits, start = start, fixv = vfixcolumn,
# consb = bfixcolumn, sim.d = FALSE, sim.zeta = FALSE,
# iteration = 5000, ctrl.mcmc = list(print = FALSE))
# }
#
# # SIMULATE DATA FOR THE OUTCOME MODEL
# # individual effects
# beta <- c(2,1,1.5)
# # contextual effects
# gamma <- c(5,-3)
# # endogenous effects
# alpha <- 0.4
# # std-dev errors
# se <- 1
# # covariates
# X <- cbind(rnorm(sum(N),0,5),rpois(sum(N),7))
# # Normalise G0
# G0norm <- norm.network(G0)
# # GX
# GX <- peer.avg(G0norm, X)
# # simulate dependent variable use an external package
# y <- CDatanet::simsar(~ X + GX, Glist = G0norm,
# theta = c(alpha, beta, gamma, se))
# y <- y$y
# # dataset
# dataset <- as.data.frame(cbind(y, X1 = X[,1], X2 = X[,2]))
## ----bayesard4s, eval=FALSE---------------------------------------------------
# mlinks <- list(model = "latent space", estimates = estimates,
# mlinks.data = list.trait, obsARD = obARD)
# out.lspa2 <- mcmcSAR(formula = y ~ X1 + X2, contextual = TRUE, G0.obs = "none",
# data = dataset, mlinks = mlinks, iteration = 2e4)
# summary(out.lspa2)
## ----Bayesard4b, echo=FALSE---------------------------------------------------
summary(out.lspa2)
## ----Bayesard4c, fig.height = 3, fig.align = "center"-------------------------
plot(out.lspa2, plot.type = "sim", mar = c(3, 2.1, 1, 1))
## ----selb1, echo = TRUE, eval=FALSE-------------------------------------------
# rm(list = ls())
# library(PartialNetwork)
# library(dplyr)
# set.seed(123)
# # Number of groups
# M <- 50
# # size of each group
# N <- rep(30,M)
# # individual effects
# beta <- c(2, 1, 1.5)
# # contextual effects
# gamma <- c(5,-3)
# # endogenous effects
# alpha <- 0.4
# # std-dev errors
# se <- 2
# # parameters of the network formation model
# rho <- c(-0.5, -.5, .4)
# # covariates
# X <- cbind(rnorm(sum(N),0,5),rpois(sum(N),7))
# # compute distance between individuals
# tmp <- c(0, cumsum(N))
# X1l <- lapply(1:M, function(x) X[c(tmp[x] + 1):tmp[x+1],1])
# X2l <- lapply(1:M, function(x) X[c(tmp[x] + 1):tmp[x+1],2])
# dist.net <- function(x, y) abs(x - y)
# X1.mat <- lapply(1:M, function(m) {
# matrix(kronecker(X1l[[m]], X1l[[m]], FUN = dist.net), N[m])})
# X2.mat <- lapply(1:M, function(m) {
# matrix(kronecker(X2l[[m]], X2l[[m]], FUN = dist.net), N[m])})
# # true network
# Xnet <- as.matrix(cbind("Const" = 1,
# "dX1" = mat.to.vec(X1.mat),
# "dX2" = mat.to.vec(X2.mat)))
# ynet <- Xnet %*% rho
# ynet <- c(1*((ynet + rlogis(length(ynet))) > 0))
# G0 <- vec.to.mat(ynet, N, normalise = FALSE)
# # number of friends
# nfriends <- unlist(lapply(G0, function(x) rowSums(x)))
# # number of missing links
# nmislink <- sapply(nfriends, function(x) sample(0:x, 1))
## ----selb2, echo = TRUE, eval=FALSE-------------------------------------------
# Gobs <- list(M) # The observed network
# G0.obs <- list(M) # Which information is true and doubtful
# for(x in 1:M){
# Gx <- G0[[x]]
# G0.obsx <- matrix(1, N[x], N[x]); diag(G0.obsx) <- 0
# csum <- cumsum(c(0, N))
# nmis <- nmislink[(csum[x] + 1):csum[x + 1]]
# for (i in 1:N[x]) {
# if(nmis[i] > 0){
# tmp <- which(c(Gx[i,]) == 1)
# if(length(which(c(Gx[i,]) == 1)) > 1) {
# tmp <- sample(which(c(Gx[i,]) == 1), nmis[i])
# }
# Gx[i,tmp] <- 0
# G0.obsx[i,] <- 0
# }
# }
# Gobs[[x]] <- Gx
# G0.obs[[x]] <- G0.obsx
# }
# G0norm <- norm.network(G0)
# # GX
# GX <- peer.avg(G0norm, X)
# # simulate dependent variable use an external package
# y <- CDatanet::simsar(~ X + GX, Glist = G0norm,
# theta = c(alpha, beta, gamma, se))
# y <- y$y
# # data set
# dataset <- as.data.frame(cbind(y, X1 = X[,1], X2 = X[,2]))
## ---- selb3a, eval=FALSE, echo=TRUE-------------------------------------------
# mlinks <- list(model = "logit", mlinks.formula = ~ dX1 + dX2,
# mlinks.data = as.data.frame(Xnet))
# out.selb1 <- mcmcSAR(formula = y ~ X1 + X2, contextual = TRUE, G0 = Gobs,
# G0.obs = G0.obs, data = dataset, mlinks = mlinks,
# iteration = 2e4)
# summary(out.selb1)
## ---- selb3b, echo=FALSE------------------------------------------------------
print(sout.selb1)
## ---- selb3c, eval=FALSE, echo=TRUE-------------------------------------------
# G0.obsvec <- as.logical(mat.to.vec(G0.obs))
# Gvec <- mat.to.vec(Gobs, ceiled = TRUE)[G0.obsvec]
# W <- unlist(data.frame(nfriends = nfriends, nmislink = nmislink) %>%
# group_by(nfriends) %>%
# summarise(w = length(nmislink)/sum(nmislink == 0)) %>%
# select(w))
# W <- lapply(1:M, function(x){
# matrix(rep(W[rowSums(G0[[x]]) + 1], each = N[x]), N[x], byrow = TRUE)})
# weights <- mat.to.vec(W)[G0.obsvec]
# mlinks <- list(model = "logit", mlinks.formula = ~ dX1 + dX2,
# mlinks.data = as.data.frame(Xnet), weights = weights)
# out.selb2 <- mcmcSAR(formula = y ~ X1 + X2, contextual = TRUE, G0 = Gobs,
# G0.obs = G0.obs, data = dataset, mlinks = mlinks,
# iteration = 2e4)
# summary(out.selb2)
## ---- selb3d, echo=FALSE------------------------------------------------------
print(sout.selb2)