Nothing
#' @title Estimate Count Data Model with Social Interactions using NPL Method
#' @param formula an object of class \link[stats]{formula}: a symbolic description of the model. The `formula` should be as for example \code{y ~ x1 + x2 | x1 + x2}
#' where `y` is the endogenous vector, the listed variables before the pipe, `x1`, `x2` are the individual exogenous variables and
#' the listed variables after the pipe, `x1`, `x2` are the contextual observable variables. Other formulas may be
#' \code{y ~ x1 + x2} for the model without contextual effects, \code{y ~ -1 + x1 + x2 | x1 + x2} for the model
#' without intercept or \code{ y ~ x1 + x2 | x2 + x3} to allow the contextual variable to be different from the individual variables.
#' @param contextual (optional) logical; if true, this means that all individual variables will be set as contextual variables. Set the
#' the `formula` as `y ~ x1 + x2` and `contextual` as `TRUE` is equivalent to set the formula as `y ~ x1 + x2 | x1 + x2`.
#' @param Glist the adjacency matrix or list sub-adjacency matrix.
#' @param Rbar the value of Rbar. If not provided, it is automatically set at \code{quantile(y, 0.9)}.
#' @param estim.rho indicates if the parameter \eqn{\rho} should be estimated or set to zero.
#' @param starting (optional) starting value of \eqn{\theta = (\lambda, \beta', \gamma')'}, \eqn{\bar{\delta}}{deltabar}, \eqn{\delta = (\delta_2, ..., \delta_{\bar{R}})}{\delta = (\delta_2, ..., \delta_{Rbar})}, and \eqn{\rho}. The parameter \eqn{\gamma} should be removed if the model
#' does not contain contextual effects (see details).
#' @param yb0 (optional) expectation of y.
#' @param optimizer is either `fastlbfgs` (L-BFGS optimization method of the package \pkg{RcppNumerical}), `nlm` (referring to the function \link[stats]{nlm}), or `optim` (referring to the function \link[stats]{optim}).
#' Other arguments
#' of these functions such as, `control` and `method` can be defined through the argument `opt.ctr`.
#' @param npl.ctr list of controls for the NPL method (see details).
#' @param opt.ctr list of arguments to be passed in `optim_lbfgs` of the package \pkg{RcppNumerical}, \link[stats]{nlm} or \link[stats]{optim} (the solver set in `optimizer`), such as `maxit`, `eps_f`, `eps_g`, `control`, `method`, ...
#' @param cov a Boolean indicating if the covariance should be computed.
#' @param data an optional data frame, list or environment (or object coercible by \link[base]{as.data.frame} to a data frame) containing the variables
#' in the model. If not found in data, the variables are taken from \code{environment(formula)}, typically the environment from which `cdnet` is called.
#' @return A list consisting of:
#' \item{info}{list of general information about the model.}
#' \item{estimate}{NPL estimator.}
#' \item{yb}{ybar (see details), expectation of y.}
#' \item{Gyb}{average of the expectation of y among friends.}
#' \item{cov}{list of covariance matrices.}
#' \item{details}{step-by-step output as returned by the optimizer.}
#' @description
#' `cdnet` is used to estimate peer effects on counting data with rational expectations (see details). The model is presented in Houndetoungan (2022).
#' @details
#' ## Model
#' Following Houndetoungan (2022), the count data \eqn{\mathbf{y}}{y} is generated from a latent variable \eqn{\mathbf{y}^*}{ys}.
#' The latent variable is given for all i as
#' \deqn{y_i^* = \lambda \mathbf{g}_i \mathbf{E}(\bar{\mathbf{y}}|\mathbf{X},\mathbf{G}) + \mathbf{x}_i'\beta + \mathbf{g}_i\mathbf{X}\gamma + \epsilon_i,}{ys_i = \lambda g_i*E(ybar|X, G) + x_i'\beta + g_i*X\gamma + \epsilon_i,}
#' where \eqn{\epsilon_i \sim N(0, 1)}{\epsilon_i --> N(0, 1)}.\cr
#' Then, \eqn{y_i = r} iff \eqn{a_r \leq y_i^* \leq a_{r+1}}{a_r \le ys_i \le a_{r + 1}}, where
#' \eqn{a_0 = -\inf}{a_0 = -Inf}, \eqn{a_1 = 0}, \eqn{a_r = \sum_{k = 1}^r\delta_k}{a_r = \delta_1 + ... + \delta_r}.
#' The parameter are subject to the constraints \eqn{\delta_r \geq \lambda}{\delta_r \ge \lambda} if \eqn{1 \leq r \leq \bar{R}}{1 \le r \le Rbar}, and
#' \eqn{\delta_r = (r - \bar{R})^{\rho}\bar{\delta} + \lambda}{a_r = deltabar*(r - Rbar)^{\rho} + \lambda} if \eqn{r \geq \bar{R} + 1}{r \ge Rbar + 1}. The unknown parameters to be estimated are
#' \eqn{\lambda}, \eqn{\beta}, \eqn{\gamma}, \eqn{\delta_2}, ..., \eqn{\delta_{\bar{R}}}{\delta_{Rbar}}, \eqn{\bar{\delta}}{deltabar}, and \eqn{\rho}.
#' ## \code{npl.ctr}
#' The model parameters is estimated using the Nested Partial Likelihood (NPL) method. This approach
#' starts with a guess of \eqn{\theta} and \eqn{\bar{y}}{yb} and constructs iteratively a sequence
#' of \eqn{\theta} and \eqn{\bar{y}}{yb}. The solution converges when the \eqn{L_1}{L} distance
#' between two consecutive \eqn{\theta} and \eqn{\bar{y}}{yb} is less than a tolerance. \cr
#' The argument \code{npl.ctr} is an optional list which contain
#' \itemize{
#' \item{tol}{ the tolerance of the NPL algorithm (default 1e-4),}
#' \item{maxit}{ the maximal number of iterations allowed (default 500),}
#' \item{print}{ a boolean indicating if the estimate should be printed at each step.}
#' \item{S}{ the number of simulation performed use to compute integral in the covariance by important sampling.}
#' }
#' @references
#' Houndetoungan, E. A. (2022). Count Data Models with Social Interactions under Rational Expectations. Available at SSRN 3721250, \doi{10.2139/ssrn.3721250}.
#' @seealso \code{\link{sart}}, \code{\link{sar}}, \code{\link{simcdnet}}.
#' @examples
#' \donttest{
#' set.seed(123)
#' # Groups' size
#' nvec <- rep(100, 2)
#' M <- length(nvec)
#' n <- sum(nvec)
#'
#' # Parameters
#' lambda <- 0.4
#' beta <- c(1.5, 2.2, -0.9)
#' gamma <- c(1.5, -1.2)
#' delta <- c(1, 0.87, 0.75, 0.6)
#' delbar <- 0.05
#' rho <- 0.5
#' theta <- c(lambda, beta, gamma)
#'
#' # X
#' X <- cbind(rnorm(n, 1, 1), rexp(n, 0.4))
#'
#' # Network
#' Glist <- list()
#'
#' for (m in 1:M) {
#' nm <- nvec[m]
#' Gm <- matrix(0, nm, nm)
#' max_d <- 30
#' for (i in 1:nm) {
#' tmp <- sample((1:nm)[-i], sample(0:max_d, 1))
#' Gm[i, tmp] <- 1
#' }
#' rs <- rowSums(Gm); rs[rs == 0] <- 1
#' Gm <- Gm/rs
#' Glist[[m]] <- Gm
#' }
#'
#'
#' # data
#' data <- data.frame(x1 = X[,1], x2 = X[,2])
#'
#' ytmp <- simcdnet(formula = ~ x1 + x2 | x1 + x2, Glist = Glist, theta = theta,
#' deltabar = delbar, delta = delta, rho = rho, data = data)
#'
#' y <- ytmp$y
#'
#' # plot histogram
#' hist(y, breaks = max(y))
#'
#' data <- data.frame(yt = y, x1 = data$x1, x2 = data$x2)
#' rm(list = ls()[!(ls() %in% c("Glist", "data"))])
#'
#' out <- cdnet(formula = yt ~ x1 + x2, contextual = TRUE, Glist = Glist,
#' data = data, Rbar = 5, estim.rho = TRUE, optimizer = "nlm")
#' summary(out)}
#' @importFrom stats quantile
#' @importFrom utils head
#' @importFrom utils tail
#' @export
cdnet <- function(formula,
contextual,
Glist,
Rbar = NULL,
estim.rho = FALSE,
starting = list(theta = NULL, deltabar = NULL, delta = NULL, rho = NULL),
yb0 = NULL,
optimizer = "fastlbfgs",
npl.ctr = list(),
opt.ctr = list(),
cov = TRUE,
data) {
stopifnot(optimizer %in% c("fastlbfgs", "optim", "nlm"))
env.formula <- environment(formula)
# controls
npl.print <- npl.ctr$print
npl.tol <- npl.ctr$tol
npl.maxit <- npl.ctr$maxit
npl.S <- npl.ctr$S
npl.incdit <- npl.ctr$incdit
if (is.null(npl.print)) {
npl.print <- TRUE
}
if (is.null(npl.tol)) {
npl.tol <- 1e-4
}
if (is.null(npl.maxit)) {
npl.maxit <- 500L
}
if (is.null(npl.S)) {
npl.S <- 1e3L
}
if (is.null(npl.incdit)) {
npl.incdit<- 30L
}
# starting values
thetat <- starting$theta
Deltat <- starting$delta
Deltabart <- starting$deltabar
rhot <- starting$rho
tmp <- c(is.null(thetat), is.null(Deltat), is.null(Deltabart))
if(estim.rho){
tmp <- c(tmp, is.null(rhot))
}
if(all(tmp) != any(tmp)){
stop("All the parameters in starting should be either NULL or defined. Set delta = numeric() if Rbar = 1.")
}
# data
if (missing(contextual)) {
contextual <- FALSE
}
if (!is.list(Glist)) {
Glist <- list(Glist)
}
M <- length(Glist)
nvec <- unlist(lapply(Glist, nrow))
n <- sum(nvec)
igr <- matrix(c(cumsum(c(0, nvec[-M])), cumsum(nvec) - 1), ncol = 2)
f.t.data <- formula.to.data(formula, contextual, Glist, M, igr, data, theta0 = thetat)
formula <- f.t.data$formula
y <- f.t.data$y
X <- f.t.data$X
coln <- c("lambda", colnames(X))
maxy <- max(y)
K <- ncol(X)
# # simulations
# Ksimu <- 0
# nsimu <- 0
# if(!is.null(simu1)) {
# Ksimu <- 1
# nsimu <- ncol(simu1)
# stopifnot(nrow(simu1) == n)
# coln <- c(coln, "simu1")
# } else{
# if(!is.null(simu2)){
# stop("simu1 = NULL whereas simu2 != NULL")
# }
# }
# if(!is.null(simu2)) {
# Ksimu <- 2
# stopifnot(nrow(simu1) == n)
# coln <- c(coln, "simu2")
# }
# yb
ybt <- rep(0, n)
if (!is.null(yb0)) {
ybt <- yb0
}
# Gybt
Gybt <- unlist(lapply(1:M, function(x) Glist[[x]] %*% ybt[(igr[x,1]:igr[x,2])+1]))
# Rbar and delta
if(is.null(Rbar)){
if(is.null(Deltat)){
Rbar <- max(quantile(y, 0.95), 1)
} else {
Rbar <- length(Deltat) + 1
}
} else {
if(!is.null(Deltat)){
if(Rbar != (length(Deltat) + 1)) stop("Rbar != length(delta) + 1")
}
}
stopifnot(Rbar <= max(y))
stopifnot(Rbar >= 1)
# check starting and computed them if not defined
lmbd0 <- NULL
if (!is.null(thetat)) {
if(length(thetat) != (K + 1)) {
stop("Length of theta is not suited.")
}
if(thetat[1] <= 0) stop("Starting lambda is not positive.")
lmbd0 <- thetat[1]
thetat <- c(log(thetat[1]/(1 -thetat[1])), thetat[-1])
if(Rbar > 1) if(any(Deltat < lmbd0)) stop("Multiple equilibrium issue: Starting lambda is greater than starting delta.")
if(Deltabart <= 0) stop("Starting deltabar is not positive.")
if(estim.rho){
if(rhot <= 0) stop("Starting rho is not positive.")
}
} else {
Xtmp <- cbind(f.t.data$Gy, X)
b <- solve(t(Xtmp)%*%Xtmp, t(Xtmp)%*%y)
b <- b/sqrt(sum((y - Xtmp%*%b)^2)/n)
b[1] <- max(b[1], 0.01)
lmbd0 <- b[1]
Deltat <- rep(b[1] + 0.01, Rbar - 1)
thetat <- c(log(b[1]/(1 - b[1])), b[-1])
Deltabart <- 0.01
rhot <- 0.01
}
lDelta <- log(Deltat - lmbd0)
lDeltabart <- log(Deltabart)
thetat <- c(thetat, lDelta, lDeltabart)
if(estim.rho){
lrho <- log(rhot)
thetat <- c(thetat, lrho)
}
# other variables
cont <- TRUE
t <- 0
Rmax <- NULL
theta <- NULL
covout <- NULL
REt <- NULL
llht <- 0
par0 <- NULL
par1 <- NULL
like <- NULL
var.comp <- NULL
steps <- list()
# arguments
ctr <- c(list(Gyb = Gybt, X = X, Rbar = Rbar, maxy = maxy, K = K, n = n, y = y), opt.ctr)
# if(Ksimu == 1) {
# ctr <- c(ctr, list(Simu1 = simu1, nsimu = nsimu))
# }
# if(Ksimu == 2) {
# ctr <- c(ctr, list(Simu1 = simu1, Simu2 = simu2, nsimu = nsimu))
# }
# Arguments used in the optimizer
if (optimizer == "fastlbfgs"){
ctr <- c(ctr, list(par = thetat));
if(estim.rho){optimizer <- "cdnetLBFGSrho"} else{optimizer <- "cdnetLBFGS"}
par0 <- "par"
par1 <- "par"
like <- "value"
} else if (optimizer == "optim") {
# ctr <- c(ctr, list(fn = ifelse(Ksimu == 0, foptimREM_NPL,
# ifelse(Ksimu == 1, foptimREM_NPLncond1,
# foptimREM_NPLncond2)), par = thetat))
if(estim.rho){
ctr <- c(ctr, list(fn = foptimREM_NPL2, par = thetat))
} else{
ctr <- c(ctr, list(fn = foptimREM_NPL, par = thetat))
}
par0 <- "par"
par1 <- "par"
like <- "value"
} else {
# ctr <- c(ctr, list(f = ifelse(Ksimu == 0, foptimREM_NPL,
# ifelse(Ksimu == 1, foptimREM_NPLncond1,
# foptimREM_NPLncond2)), p = thetat))
if(estim.rho){
ctr <- c(ctr, list(f = foptimREM_NPL2, p = thetat))
} else{
ctr <- c(ctr, list(f = foptimREM_NPL, p = thetat))
}
par0 <- "p"
par1 <- "estimate"
like <- "minimum"
}
# Fonction used to compute new candidate for E(y)
fL_NPLR <- NULL
fcovCDIR <- NULL
if(estim.rho){
fL_NPLR <- fL_NPL2
fcovCDI <- fcovCDI2
} else{
fL_NPLR <- fL_NPL
fcovCDI <- fcovCDI
}
# NPL algorithm
ninc.d <- 0
dist0 <- Inf
if(npl.print) {
while(cont) {
# tryCatch({
ybt0 <- ybt + 0 #copy in different memory
# compute theta
REt <- do.call(get(optimizer), ctr)
thetat <- REt[[par1]]
llht <- -REt[[like]]
theta <- c(1/(1 + exp(-thetat[1])), thetat[2:(K +1)], exp(thetat[-(1:(K +1))]))
if(Rbar > 1){
theta[(K+2):(K+Rbar)] <- theta[(K+2):(K+Rbar)] + theta[1]
}
# compute y
fL_NPLR(ybt, Gybt, Glist, igr, M, X, thetat, Rbar, K, n)
# distance
# print(theta)
var.eps <- abs(c((ctr[[par0]] - thetat)/thetat, (ybt0 - ybt)/ybt))
if(all(!is.finite(var.eps))){
var.eps <- abs(c(ctr[[par0]] - thetat, ybt0 - ybt))
}
dist <- max(var.eps, na.rm = TRUE)
ninc.d <- (ninc.d + 1)*(dist > dist0) #counts the successive number of times distance increases
dist0 <- dist
cont <- (dist > npl.tol & t < (npl.maxit - 1))
t <- t + 1
REt$dist <- dist
ctr[[par0]] <- thetat
steps[[t]] <- REt
cat("---------------\n")
cat(paste0("Step : ", t), "\n")
cat(paste0("Distance : ", round(dist,3)), "\n")
cat(paste0("Likelihood : ", round(llht,3)), "\n")
cat("Estimate:", "\n")
print(theta)
if((ninc.d > npl.incdit) | (llht < -1e293)) {
cat("** Non-convergence ** Redefining theta and computing a new yb\n")
thetat[1] <- -4.5
dist0 <- Inf
if(estim.rho){
fnewyb2(ybt, Gybt, Glist, igr, M, X, thetat, Rbar, K, n, npl.tol, npl.maxit)
} else{
fnewyb(ybt, Gybt, Glist, igr, M, X, thetat, Rbar, K, n, npl.tol, npl.maxit)
}
ctr[[par0]] <- thetat
}
# },
# error = function(e){
# cat("** Non-convergence ** Redefining theta and computing a new yb\n")
# thetat[1] <- -4.5
# dist0 <- Inf
# if(estim.rho){
# fnewyb2(ybt, Gybt, Glist, igr, M, X, thetat, Rbar, K, n, npl.tol, npl.maxit)
# } else{
# fnewyb(ybt, Gybt, Glist, igr, M, X, thetat, Rbar, K, n, npl.tol, npl.maxit)
# }
# cont <- TRUE
# t <- t + 1
# ctr[[par0]] <- thetat
# steps[[t]] <- NULL
# })
}
} else {
while(cont) {
# tryCatch({
ybt0 <- ybt + 0 #copy in different memory
# compute theta
REt <- do.call(get(optimizer), ctr)
thetat <- REt[[par1]]
llht <- -REt[[like]]
# compute y
fL_NPLR(ybt, Gybt, Glist, igr, M, X, thetat, Rbar, K, n)
# distance
var.eps <- abs(c((ctr[[par0]] - thetat)/thetat, (ybt0 - ybt)/ybt))
if(all(!is.finite(var.eps))){
var.eps <- abs(c(ctr[[par0]] - thetat, ybt0 - ybt))
}
dist <- max(var.eps, na.rm = TRUE)
ninc.d <- (ninc.d + 1)*(dist > dist0) #counts the successive number of times distance increases
dist0 <- dist
cont <- (dist > npl.tol & t < (npl.maxit - 1))
t <- t + 1
REt$dist <- dist
ctr[[par0]] <- thetat
steps[[t]] <- REt
if((ninc.d > npl.incdit) | (llht < -1e293)) {
thetat[1] <- -4.5
dist0 <- Inf
if(estim.rho){
fnewyb2(ybt, Gybt, Glist, igr, M, X, thetat, Rbar, K, n, npl.tol, npl.maxit)
} else{
fnewyb(ybt, Gybt, Glist, igr, M, X, thetat, Rbar, K, n, npl.tol, npl.maxit)
}
ctr[[par0]] <- thetat
}
# },
# error = function(e){
# thetat[1] <- -4.5
# dist0 <- Inf
# if(estim.rho){
# fnewyb2(ybt, Gybt, Glist, igr, M, X, thetat, Rbar, K, n, npl.tol, npl.maxit)
# } else{
# fnewyb(ybt, Gybt, Glist, igr, M, X, thetat, Rbar, K, n, npl.tol, npl.maxit)
# }
# cont <- TRUE
# t <- t + 1
# ctr[[par0]] <- thetat
# steps[[t]] <- NULL
# })
}
theta <- c(1/(1 + exp(-thetat[1])), thetat[2:(K +1)], exp(thetat[-(1:(K +1))]))
if(Rbar > 1){
theta[(K+2):(K+Rbar)] <- theta[(K+2):(K+Rbar)] + theta[1]
}
}
# name theta
namtheta <- coln
if(Rbar > 1){
namtheta <- c(namtheta, paste0("delta", 2:(Rbar)))
}
namtheta <- c(namtheta, "deltabar")
if(estim.rho){
namtheta <- c(namtheta, "rho")
}
names(theta) <- namtheta
environment(formula) <- env.formula
sdata <- list(
"formula" = formula,
"Glist" = deparse(substitute(Glist)),
"nfriends" = unlist(lapply(Glist, function(u) sum(u > 0)))
)
if (!missing(data)) {
sdata <- c(sdata, list("data" = deparse(substitute(data))))
}
if (npl.maxit == t) {
warning("The maximum number of iterations of the NPL algorithm has been reached.")
}
# covariance and ME
tmp <- fcovCDI(n, Gybt, thetat, X, Rbar, K, npl.S, Glist, igr, M, cov)
meffects <- c(tmp$meffects)
var.comp <- tmp$var.comp
covt <- tmp$covt
covm <- tmp$covm
Rmax <- tmp$Rmax
colnME <- coln
if("(Intercept)" %in% coln) {
colnME <- coln[-2]
meffects <- meffects[-2]
covm <- covm[-2, -2]
}
names(meffects) <- colnME
if(!is.null(covt)) {
namecovt <- coln
if(Rbar > 1){
namecovt <- c(namecovt, paste0("log(delta", 2:(Rbar), ")"))
}
namecovt <- c(namecovt, "log(deltabar)")
if(estim.rho){
namecovt <- c(namecovt, "log(rho)")
}
colnames(covt) <- namecovt
rownames(covt) <- namecovt
colnames(covm) <- colnME
rownames(covm) <- colnME
namecovt[1] <- "log(lambda)"
rownames(var.comp$Sigma) <- namecovt
rownames(var.comp$Omega) <- namecovt
colnames(var.comp$Sigma) <- namecovt
colnames(var.comp$Omega) <- namecovt
}
AIC <- 2*length(theta) - 2*llht
BIC <- length(theta)*log(n) - 2*llht
INFO <- list("M" = M,
"n" = n,
"Kz" = K,
"nlinks" = unlist(lapply(Glist, function(u) sum(u > 0))),
"formula" = formula,
"estim.rho" = estim.rho,
"Rbar" = Rbar,
"Rmax" = Rmax,
"log.like" = llht,
"npl.iter" = t,
"AIC" = AIC,
"BIC" = BIC)
out <- list("info" = INFO,
"estimate" = list(parms = theta, marg.effects = meffects),
"yb" = ybt,
"Gyb" = Gybt,
"cov" = list(parms = covt, marg.effects = covm, var.comp = var.comp),
"details" = steps)
class(out) <- "cdnet"
out
}
#' @title Summarize Count Data Model with Social Interactions
#' @description Summary and print methods for the class `cdnet` as returned by the function \link{cdnet}.
#' @param object an object of class `cdnet`, output of the function \code{\link{cdnet}}.
#' @param x an object of class `summary.cdnet`, output of the function \code{\link{summary.cdnet}},
#' class `summary.cdnets`, list of outputs of the function \code{\link{summary.cdnet}}
#' (when the model is estimated many times to control for the endogeneity)
#' or class `cdnet` of the function \code{\link{cdnet}}.
#' @param Glist adjacency matrix or list sub-adjacency matrix. This is not necessary if the covariance method was computed in \link{cdnet}.
#' @param data a `dataframe` containing the explanatory variables. This is not necessary if the covariance method was computed in \link{cdnet}.
#' @param S number of simulation to be used to compute integral in the covariance by important sampling.
#' @param ... further arguments passed to or from other methods.
#' @return A list of the same objects in `object`.
#' @export
"summary.cdnet" <- function(object,
Glist,
data,
S = 1e3L,
...) {
stopifnot(class(object) == "cdnet")
out <- c(object, list("..." = ...))
if(is.null(object$cov$parms)){
env.formula <- environment(object$info$formula)
Rbar <- object$info$Rbar
Kz <- object$info$Kz
estim.rho <- object$info$estim.rho
parms <- object$estimate$parms
thetat <- c(log(parms[1]/(1 -parms[1])), parms[-1])
if(Rbar > 1){
thetat[(2 + Kz):(Kz + Rbar)] <- thetat[(2 + Kz):(Kz + Rbar)] - parms[1]
}
thetat[(2 + Kz):(Kz + Rbar + 1 + estim.rho)] <- log(tail(thetat, Rbar + estim.rho))
if(estim.rho){
fcovCDI <- fcovCDI2
} else{
fcovCDI <- fcovCDI
}
Gybt <- object$Gyb
npl.S <- S
if (is.null(npl.S)) {
npl.S <- 1e3L
}
contextual <- FALSE
formula <- object$info$formula
if (!is.list(Glist)) {
Glist <- list(Glist)
}
M <- length(Glist)
nvec <- unlist(lapply(Glist, nrow))
n <- sum(nvec)
igr <- matrix(c(cumsum(c(0, nvec[-M])), cumsum(nvec) - 1), ncol = 2)
f.t.data <- formula.to.data(formula, contextual, Glist, M, igr, data, theta0 = thetat)
formula <- f.t.data$formula
y <- f.t.data$y
X <- f.t.data$X
K <- ncol(X)
coln <- c("lambda", colnames(X))
tmp <- fcovCDI(n, Gybt, thetat, X, Rbar, K, npl.S, Glist, igr, M, TRUE)
meffects <- c(tmp$meffects)
var.comp <- tmp$var.comp
covt <- tmp$covt
covm <- tmp$covm
Rmax <- tmp$Rmax
colnME <- coln
if("(Intercept)" %in% coln) {
colnME <- coln[-2]
meffects <- meffects[-2]
covm <- covm[-2, -2]
}
names(meffects) <- colnME
if(!is.null(covt)) {
namecovt <- coln
if(Rbar > 1){
namecovt <- c(namecovt, paste0("log(delta", 2:(Rbar), ")"))
}
namecovt <- c(namecovt, "log(deltabar)")
if(estim.rho){
namecovt <- c(namecovt, "log(rho)")
}
colnames(covt) <- namecovt
rownames(covt) <- namecovt
colnames(covm) <- colnME
rownames(covm) <- colnME
namecovt[1] <- "log(lambda)"
rownames(var.comp$Sigma) <- namecovt
rownames(var.comp$Omega) <- namecovt
colnames(var.comp$Sigma) <- namecovt
colnames(var.comp$Omega) <- namecovt
}
out$cov <- list(parms = covt, marg.effects = covm, var.comp = var.comp)
}
class(out) <- "summary.cdnet"
out
}
#' @rdname summary.cdnet
#' @importFrom stats pchisq
#' @export
"print.summary.cdnet" <- function(x, ...) {
stopifnot(class(x) == "summary.cdnet")
M <- x$info$M
n <- x$info$n
iteration <- x$info$npl.iter
Rbar <- x$info$Rbar
formula <- x$info$formula
Kz <- x$info$Kz
estim.rho <- x$info$estim.rho
AIC <- x$info$AIC
BIC <- x$info$BIC
parms <- x$estimate$parms
coef <- parms[1:(1 + Kz)]
K <- length(coef)
meff <- x$estimate$marg.effects
std <- sqrt(head(diag(x$cov$parms), K))
std.meff <- sqrt(diag(x$cov$marg.effects))
delta <- tail(parms, Rbar + estim.rho)
llh <- x$info$log.like
tmp <- fcoefficients(coef, std)
out_print <- tmp$out_print
out <- tmp$out
out_print <- c(list(out_print), x[-(1:6)], list(...))
tmp.meff <- fcoefficients(meff, std.meff)
out_print.meff <- tmp.meff$out_print
out.meff <- tmp.meff$out
out_print.meff <- c(list(out_print.meff), x[-(1:6)], list(...))
nfr <- x$info$nlinks
cat("Count data Model with Social Interactions\n\n")
cat("Call:\n")
print(formula)
cat("\nMethod: Nested pseudo-likelihood (NPL) \nIteration: ", iteration, sep = "", "\n\n")
cat("Network:\n")
cat("Number of groups : ", M, sep = "", "\n")
cat("Sample size : ", n, sep = "", "\n")
cat("Average number of friends: ", sum(nfr)/n, sep = "", "\n\n")
cat("Coefficients:\n")
do.call("print", out_print)
cat("\nMarginal Effects:\n")
do.call("print", out_print.meff)
cat("---\nSignif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n\n")
cat("Rbar: ", Rbar, sep = "", "\n")
if(Rbar > 1){
cat("delta:", delta[1:(Rbar - 1)], sep = " ", "\n")
}
cat("deltabar: ", delta[Rbar], sep = "")
if(estim.rho){
cat(" -- rho: ", delta[Rbar+1], sep = "", "\n")
cat("Wald test H0: rho = 0, Prob = ", 1 - pchisq(1/x$cov$parms[Kz+Rbar+2,Kz+Rbar+2], df = 1), sep = "", "\n")
} else{
cat("\n")
}
cat("log pseudo-likelihood: ", llh, sep = "", "\n")
cat("AIC: ", AIC, " -- BIC: ", BIC, sep = "", "\n")
invisible(x)
}
#' @rdname summary.cdnet
#' @export
"print.cdnet" <- function(x, ...) {
stopifnot(class(x) == "cdnet")
print(summary(x, ...))
}
#' @rdname summary.cdnet
#' @importFrom stats cov
#' @export
"summary.cdnets" <- function(object, ...) {
stopifnot(class(object) %in% c("list", "cdnets", "summary.cdnets"))
lclass <- unique(unlist(lapply(object, class)))
if (!all(lclass %in%c("summary.cdnet"))) {
stop("All the components in `object` should be from `summary.cdnet` class")
}
nsim <- length(object)
coef <- do.call("rbind", lapply(object, function(z) t(z$estimate$parms)))
meff <- do.call("rbind", lapply(object, function(z) t(z$estimate$marg.effects)))
estimate <- colSums(coef)/nsim
meffects <- colSums(meff)/nsim
vcoef2 <- Reduce("+", lapply(object, function(z) z$cov$parms))/nsim
vmeff2 <- Reduce("+", lapply(object, function(z) z$cov$marg.effects))/nsim
vcoef1 <- cov(coef)
vmeff1 <- cov(meff)
vcoef <- vcoef1 + vcoef2
vmeff <- vmeff1 + vmeff2
llh <- unlist(lapply(object, function(z) z$info$log.like))
llh <- c("min" = min(llh), "mean" = mean(llh), "max" = max(llh))
M <- object[[1]]$info$M
n <- object[[1]]$info$n
Rbar <- object[[1]]$info$Rbar
Kz <- object[[1]]$info$Kz
estim.rho <- object[[1]]$info$estim.rho
INFO <- list("M" = M,
"n" = n,
"Kz" = Kz,
"estim.rho" = estim.rho,
"log.like" = llh,
"Rbar" = Rbar,
"simulation" = nsim)
out <- list("info" = INFO,
"estimate" = list(parms = coef, marg.effects = meffects),
"cov" = list(parms = vcoef, marg.effects = vmeff),
... = ...)
class(out) <- "summary.cdnets"
out
}
#' @rdname summary.cdnet
#' @importFrom stats cov
#' @export
"print.summary.cdnets" <- function(x, ...) {
stopifnot(class(x) %in% c("summary.cdnets"))
nsim <- x$info$simulation
Kz <- x$info$Kz
estim.rho <- x$info$estim.rho
parms <- x$estimate$parms
coef <- parms[1:(1 + Kz)]
K <- length(coef)
meff <- x$estimate$marg.effects
std <- sqrt(head(diag(x$cov$parms), K))
std.meff <- sqrt(diag(x$cov$marg.effects))
delta <- tail(parms, Rbar + estim.rho)
vcoef <- x$cov$parms
vmeff <- x$cov$marg.effects
llh <- x$info$log.like
M <- x$info$M
n <- x$info$n
Rbar <- x$info$Rbar
std <- sqrt(head(diag(vcoef), K))
std.meff <- sqrt(diag(vmeff))
tmp <- fcoefficients(coef, std)
out_print <- tmp$out_print
out <- tmp$out
tmp.meff <- fcoefficients(meff, std.meff)
out_print.meff <- tmp.meff$out_print
out.meff <- tmp.meff$out
out_print <- c(list(out_print), x[-c(1:3)], list(...))
out_print.meff <- c(list(out_print.meff), x[-c(1:3)], list(...))
cat("Count data Model with Social Interactions\n\n")
cat("Method: Replication of Nested pseudo-likelihood (NPL) \nReplication: ", nsim, "\n\n")
cat("Coefficients:\n")
do.call("print", out_print)
cat("\nMarginal Effects:\n")
do.call("print", out_print.meff)
cat("---\nSignif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n\n")
cat("Rbar: ", Rbar, "\n")
if(Rbar > 1){
cat("delta: ", delta[1:(Rbar - 1)], "\n")
}
cat("deltabar: ", delta[Rbar])
if(estim.rho){
cat(" -- rho: ", delta[Rbar+1], "\n")
cat("Wald test H0: rho = 0, Prob = ", 1 - pchisq(1/x$cov$parms[Kz+Rbar+2,Kz+Rbar+2], df = 1), "\n")
} else{
cat("\n")
}
cat("log pseudo-likelihood: ", "\n")
print(llh)
invisible(x)
}
#' @rdname summary.cdnet
#' @importFrom stats cov
#' @export
"print.cdnets" <- function(x, ...) {
stopifnot(class(x) %in% c("cdnets", "list"))
print(summary.cdnets(x, ...))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.