Last active
September 12, 2017 18:39
-
-
Save nickeubank/5970facdf644639d44f9cde7116d6fd7 to your computer and use it in GitHub Desktop.
modifications to functions.R to get reciprocal graph
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| stripAttributes <- function(g) { | |
| for (a in list.graph.attributes(g)) { | |
| g <- delete_graph_attr(g, a) | |
| } | |
| for (a in list.edge.attributes(g)) { | |
| g <- delete_edge_attr(g, a) | |
| } | |
| for (a in list.vertex.attributes(g)) { | |
| if (a != "name") { | |
| g <- delete_vertex_attr(g, a) | |
| } | |
| } | |
| return(g) | |
| } | |
| makeNetwork <- function(nodes, ties, type, village) { | |
| g <- ties[ties$type == type,c("i","j")] | |
| n_no_interview <- length(unique(g$j[!g$j %in% nodes$i[nodes$interviewed]])) | |
| reciprocity_construct <- NA | |
| if (type %in% c("friend", "family", "lender", "solver", "speak")) { | |
| g <- graph.data.frame(g, vertices = nodes) | |
| if (type %in% c("friend", "family")) { | |
| reciprocity_construct <- reciprocity(g) | |
| g <- as.undirected(g, "mutual", edge.attr.comb = list("type" = "first")) | |
| } | |
| # Add friends and family that does not require | |
| # reciprocation (nick, sept 12, 2017) | |
| if (type %in% c("friends_nonrecip", "family_nonrecip")) { | |
| reciprocity_construct <- reciprocity(g) | |
| g <- as.undirected(g, "collapse", edge.attr.comb = list("type" = "first")) | |
| } | |
| g <- induced_subgraph(g, V(g)[V(g)$interviewed]) | |
| } else if (type %in% c("Leader")) { | |
| leaders <- unique(g$j) | |
| g <- g[!g$i %in% leaders,] | |
| g <- graph.data.frame(g, vertices = nodes) | |
| g <- induced_subgraph(g, V(g)[V(g)$interviewed]) | |
| low <- V(g)[degree(g,mode = "in") < 5] | |
| g <- g - E(g)[V(g) %->% low] | |
| } else { | |
| stop("Unsupported type") | |
| } | |
| g$etype <- type | |
| g$village <- village | |
| stats_construct <- data.frame(village = village, | |
| etype = type, | |
| reciprocity_construct = reciprocity_construct, | |
| n_no_interview = n_no_interview, | |
| stringsAsFactors = F) | |
| return(list(g = g, | |
| stats_construct = stats_construct)) | |
| } | |
| makeNetworks <- function(nodes, ties, village) { | |
| types <- c("family","friend","solver","lender","Leader", "speak", "friends_nonrecip", "family_nonrecip") | |
| networks <- lapply(types, makeNetwork, nodes = nodes, ties = ties, village = village) | |
| stats_construct <- lapply(networks, function(l) l$stats_construct) | |
| networks <- lapply(networks, function(l) l$g) | |
| names(networks) <- types | |
| stats_construct <- do.call(rbind, stats_construct) | |
| # distances network | |
| p <- as_data_frame(networks$family, what = "vertices")[,c("name","lon","lat","nodist")] | |
| colnames(p)[1] <- "i" | |
| sto <- as.data.frame(t(combn(p$i, 2))) | |
| colnames(sto) <- c("i","j") | |
| sto <- merge(sto, p) | |
| sto <- merge(sto, p, by.x = "j", by.y = "i") | |
| sto$weight <- sapply(1:nrow(sto), function(i) distGeo(sto[i,c("lon.x","lat.x")],sto[i,c("lon.y","lat.y")])) | |
| sto$weight[sto$nodist.x | sto$nodist.y] <- NA | |
| sto$weight[is.na(sto$weight)] <- mean(sto$weight, na.rm = T) # those that don't have coordinates get assigned mean distance | |
| sto <- sto[,c("i","j","weight")] | |
| networks$geo <- graph.data.frame(sto, | |
| directed = F, | |
| vertices = as_data_frame(networks$family, what = "vertices")) | |
| networks$geo$etype <- "geo" | |
| networks$geo$village <- village | |
| # union of all networks | |
| networks$union <- networks$family %u% | |
| stripAttributes(networks$friend) %u% | |
| as.undirected(stripAttributes(networks$solver), mode = "collapse") %u% | |
| as.undirected(stripAttributes(networks$lender), mode = "collapse") | |
| networks$union <- simplify(networks$union) | |
| networks$union$etype <- "union" | |
| networks$union$village <- village | |
| # Union without requirement of reciprocity | |
| networks$union_nonrecip <- networks$family_nonrecip %u% | |
| stripAttributes(networks$friend_nonrecip) %u% | |
| as.undirected(stripAttributes(networks$solver), mode = "collapse") %u% | |
| as.undirected(stripAttributes(networks$lender), mode = "collapse") | |
| networks$union_nonrecip <- simplify(networks$union_nonrecip) | |
| networks$union_nonrecip$etype <- "union_nonrecrip" | |
| networks$union_nonrecip$village <- village | |
| networks$strong_ties <- networks$family %u% stripAttributes(networks$friend) | |
| networks$strong_ties <- simplify(networks$strong_ties) | |
| networks$strong_ties$etype <- "strong_ties" | |
| networks$weak_ties <- networks$solver %u% stripAttributes(networks$lender) | |
| networks$weak_ties <- simplify(networks$weak_ties) | |
| networks$weak_ties$etype <- "weak_ties" | |
| networks$leader_obj <- networks$union | |
| networks$leader_obj <- networks$leader_obj - | |
| difference(E(networks$leader_obj), | |
| E(networks$leader_obj)[inc(V(networks$leader_obj)[is_leader_obj == 1])]) | |
| networks$leader_obj$etype <- "leader_obj" | |
| n_no_interview <- nrow(nodes) - length(unique(ties$j[ties$j %in% nodes$i[nodes$interviewed]])) | |
| stats_construct <- rbind(stats_construct, | |
| data.frame(village = village, | |
| etype = "union", | |
| reciprocity_construct = NA, | |
| n_no_interview = n_no_interview, | |
| stringsAsFactors = F)) | |
| return(list(networks = networks, | |
| stats_construct = stats_construct)) | |
| } | |
| cleanDf <- function(df) { | |
| cols <- c("IDLINE", "key", "connection", | |
| "exist", "id", "confirmentryid", | |
| "village", "family_name", "first_name", | |
| "alias_name", "gender", "gender_rec", | |
| "hh_quadrant", "hh_number", "confirm", | |
| "newfamily", "newfirst", "newalias", | |
| "newgender", "qa13relation", "lineperperson", | |
| "connectionid", "fewconnections") | |
| nodes <- df[,which(!colnames(df) %in% cols)] | |
| nodes <- unique(nodes) | |
| if (length(unique(nodes$resid)) != nrow(nodes)) { | |
| dup <- table(nodes$resid) | |
| dup <- names(dup)[dup>1] | |
| stop("# unique ids: ", length(unique(nodes$resid)), | |
| ", # unique individual rows: ", nrow(nodes), "\n", | |
| "IDs: ", paste(dup, collapse = ", ")) | |
| } | |
| pol_index <- as.matrix(nodes[,c("qe2b_Attended", "qe2b_Contributesproject", | |
| "qe2b_Contributedmember", "qe2b_Contributedlabour", "qe2b_Reportedlead", | |
| "qe2b_Reportedgov")]) | |
| nodes <- nodes[,c("i", "interviewed", "heard", "adopt", "age", "female", "income", "edu", "phone", "immigrant", | |
| "is_leader_obj", | |
| "sent", "pborrow", "plend", "picontact", | |
| "no_dk", "no_dksms", "no_phone", "no_worth", "no_way", | |
| "lon", "lat", "nodist", | |
| "efficacy_int", "efficacy_ext", | |
| "prosoc_dic", "prosoc_pub", "attend", "satisfaction", "satisfaction2")] | |
| nomiss <- nodes[!nodes$nodist,c("lon","lat")] | |
| nomiss <- geomean(nomiss) | |
| nodes$lon[nodes$nodist] <- nomiss[1,1] | |
| nodes$lat[nodes$nodist] <- nomiss[1,2] | |
| nodes$adopt[is.na(nodes$adopt)] <- F | |
| pol_index <- scale(pol_index) | |
| weights <- cov(pol_index) | |
| weights <- solve(weights) | |
| weights <- colSums(weights) | |
| pol_index <- apply(pol_index, 1, weighted.mean, w = weights) | |
| nodes$pol_index <- pol_index | |
| ties <- df[,c("resid","id","connection")] | |
| number <- regexpr("[0-9]", ties$connection) | |
| # ties$order <- ifelse(number == -1, NA, as.numeric(substr(ties$connection, number, number))) | |
| ties$type <- ifelse(number == -1, ties$connection, substr(ties$connection, 1, number-1)) | |
| ties$connection <- NULL | |
| colnames(ties) <- c("i", "j", "type") | |
| ties$type[ties$type == "member"] <- "family" | |
| ties <- ties[ties$i != ties$j,] | |
| ties <- ties[!duplicated(ties),] | |
| # add in non-interviewed nodes | |
| missing <- unique(ties$j[!ties$j %in% nodes$i]) | |
| nodes <- rbind(nodes, | |
| data.frame(i = missing, | |
| interviewed = F, | |
| heard = NA, | |
| adopt = NA, | |
| age = NA, | |
| female = NA, | |
| income = NA, | |
| edu = NA, | |
| phone = NA, | |
| immigrant = NA, | |
| is_leader_obj = NA, | |
| lat = NA, | |
| lon = NA, | |
| nodist = NA, | |
| efficacy_int = NA, | |
| efficacy_ext = NA, | |
| pol_index = NA, | |
| prosoc_dic = NA, | |
| prosoc_pub = NA, | |
| attend = NA, | |
| satisfaction = NA, | |
| satisfaction2 = NA, | |
| sent = NA, pborrow = NA, plend = NA, picontact = NA, | |
| no_dk = NA, no_dksms = NA, no_phone = NA, no_worth = NA, no_way = NA), | |
| stringsAsFactors = F) | |
| leaders <- unique(ties$j[ties$type == "Leader"]) | |
| # nodes$is_leader <- nodes$i %in% leaders | |
| return(list(nodes = nodes, | |
| ties = ties)) | |
| } | |
| getStatistics <- function(g) { | |
| comp_g <- components(g) | |
| g_large <- groups(comp_g) | |
| # g_large <- induced.subgraph(g, g_large[[which.max(comp_g$csize)]]) | |
| if(g$etype %in% c("family","friend","solver","lender","union")) { | |
| graph_level <- data.frame(village = g$village, | |
| etype = g$etype, | |
| n_nodes = vcount(g), | |
| n_ties = ecount(g), | |
| diameter = diameter(g), | |
| mean_deg = mean(degree(g)), | |
| mean_deg_in = ifelse(is.directed(g), mean(degree(g, mode = "in")), NA), | |
| mean_deg_out = ifelse(is.directed(g), mean(degree(g, mode = "out")), NA), | |
| density = graph.density(g), | |
| mean_dist = mean_distance(g), | |
| clustering_coef = transitivity(g, type = "global"), | |
| n_ntrivial_comp = length(comp_g$csize[comp_g$csize > 1]), | |
| n_lcomp = max(comp_g$csize), | |
| pct_lcomp = max(comp_g$csize)/vcount(g), | |
| n_isolates = length(degree(g)[degree(g)==0]), | |
| pct_isolates = length(degree(g)[degree(g)==0])/vcount(g), | |
| reciprocity = ifelse(is.directed(g), reciprocity(g), NA), | |
| stringsAsFactors = F) | |
| node_level <- data.frame(village = g$village, | |
| etype = g$etype, | |
| i = V(g)$name, | |
| degree = degree(g), | |
| degree_in = ifelse(is.directed(g), degree(g, mode = "in"), NA), | |
| degree_out = ifelse(is.directed(g), degree(g, mode = "out"), NA), | |
| clustering_coef = transitivity(g, type = "local"), | |
| between = centr_betw(g)$res, | |
| close = centr_clo(g)$res, | |
| eigen = centr_eigen(g)$vector, | |
| stringsAsFactors = F) | |
| } else if(g$etype == "Leader") { | |
| graph_level <- data.frame(village = g$village, | |
| etype = g$etype, | |
| n_nodes = vcount(g), | |
| n_ties = ecount(g), | |
| n_ntrivial_comp = length(comp_g$csize[comp_g$csize > 1]), | |
| n_lcomp = max(comp_g$csize), | |
| pct_lcomp = max(comp_g$csize)/vcount(g), | |
| n_isolates = length(degree(g)[degree(g,mode="all")==0]), | |
| pct_isolates = length(degree(g)[degree(g,mode="all")==0])/(vcount(g)-length(comp_g$csize[comp_g$csize > 1]))) | |
| node_level <- NULL | |
| } else { | |
| graph_level <- NULL | |
| node_level <- NULL | |
| } | |
| return(list(graph_level = graph_level, | |
| node_level = node_level)) | |
| } | |
| processVillage <- function(df, village) { | |
| df <- df[df$village_rec == village,] | |
| df <- cleanDf(df) | |
| networks <- makeNetworks(df$nodes, df$ties, village) | |
| stats_construct <- networks$stats_construct | |
| networks <- networks$networks | |
| etypes <- names(networks) | |
| stats <- lapply(networks, getStatistics) | |
| stats_lead <- stats[[which(etypes == "Leader")]]$graph_level | |
| rownames(stats_lead) <- NULL | |
| stats_graph <- lapply(stats[etypes != "Leader"], function(l) l$graph_level) | |
| stats_graph <- do.call(rbind, stats_graph) | |
| stats_graph <- merge(stats_graph, stats_construct) | |
| stats_graph$pct_no_interview <- stats_graph$n_no_interview / (stats_graph$n_no_interview + stats_graph$n_nodes) | |
| rownames(stats_graph) <- NULL | |
| stats_node <- lapply(stats[etypes != "Leader"], function(l) l$node_level) | |
| stats_node <- do.call(rbind, stats_node) | |
| rownames(stats_node) <- NULL | |
| return(list(village = village, | |
| networks = networks, | |
| stats_graph = stats_graph, | |
| stats_node = stats_node, | |
| stats_lead = stats_lead)) | |
| } | |
| processVillages <- function(df) { | |
| villages <- unique(d$village_rec) | |
| v <- lapply(villages, function(v) tryCatch(processVillage(df,v), | |
| error = function(m) return(m))) | |
| names(v) <- villages | |
| return(v) | |
| } | |
| prepareVillage <- function(village) { | |
| if(is.null(village)) return(NULL) | |
| n <- list(union = list(as_adjacency_matrix(village$networks$union, sparse = F)), | |
| solver = list(as_adjacency_matrix(village$networks$solver, sparse = F)), | |
| lender = list(as_adjacency_matrix(village$networks$lender, sparse = F)), | |
| family = list(as_adjacency_matrix(village$networks$family, sparse = F)), | |
| friend = list(as_adjacency_matrix(village$networks$friend, sparse = F)), | |
| leader = list(as_adjacency_matrix(village$networks$Leader, sparse = F)), | |
| speak = list(as_adjacency_matrix(village$networks$speak, sparse = F)), | |
| strong_ties = list(as_adjacency_matrix(village$networks$strong_ties, sparse = F)), | |
| weak_ties = list(as_adjacency_matrix(village$networks$weak_ties, sparse = F)), | |
| leader_obj = list(as_adjacency_matrix(village$networks$leader_obj, sparse = F)), | |
| geo = list(as_adjacency_matrix(village$networks$geo, sparse = F, attr = "weight")), | |
| df = as_data_frame(village$networks$union, what = "vertices")) | |
| n$df$heard <- as.numeric(n$df$heard) | |
| n$df$adopt <- as.numeric(n$df$adopt) | |
| n$df$is_leader <- degree(village$networks$Leader, mode = "in") > 0 | |
| n$heard <- data.frame(t1 = n$df$heard, row.names = rownames(n$df)) | |
| n$adopt <- data.frame(t1 = n$df$adopt, row.names = rownames(n$df)) | |
| return(n) | |
| } | |
| prepare <- function(villages) { | |
| out <- lapply(villages, prepareVillage) | |
| names(out) <- names(villages) | |
| return(out) | |
| } | |
| fit <- function(y, formula.glm, formula.tnam, data, sub = NULL, rename = NULL, ...) { | |
| if(!is.null(sub)) data <- data[sub] | |
| df.glm <- lapply(data, function(l) { | |
| df.glm <- model.frame(formula.glm, l$df) | |
| for (i in 1:length(colnames(df.glm))) { | |
| if(is.logical(df.glm[,i])) df.glm[,i] <- as.numeric(df.glm[,i]) | |
| } | |
| df.glm$node <- rownames(df.glm) | |
| return(df.glm) | |
| }) | |
| df.glm <- do.call(rbind, df.glm) | |
| if(!is.null(formula.tnam)) { | |
| df.tnam <- lapply(data, function(n) { | |
| makeCol <- function(v) { | |
| out <- data.frame(t1 = v) | |
| rownames(out) <- rownames(n$df) | |
| out | |
| } | |
| env <- new.env() | |
| isActive <- as.integer(n$df$pol_index > 0.1625646) | |
| assign("uni", n$union, envir = env) | |
| assign("solver", n$solver, envir = env) | |
| assign("lender", n$lender, envir = env) | |
| assign("family", n$family, envir = env) | |
| assign("friend", n$friend, envir = env) | |
| assign("leader", n$leader, envir = env) | |
| assign("leader_obj", n$leader_obj, envir = env) | |
| assign("is_leader", makeCol(n$df$is_leader), envir = env) | |
| assign("is_leader_obj", makeCol(n$df$is_leader_obj), envir = env) | |
| assign("is_peer_obj", makeCol(1-n$df$is_leader_obj), envir = env) | |
| assign("is_active", makeCol(isActive), envir = env) | |
| assign("is_inactive", makeCol(1-isActive), envir = env) | |
| assign("geo", n$geo, envir = env) | |
| assign("heard", n$heard, envir = env) | |
| assign("adopt", n$adopt, envir = env) | |
| assign("speak", n$speak, envir = env) | |
| satisfaction <- data.frame(t1 = n$df$satisfaction, row.names = n$df$name) | |
| satisfaction2 <- data.frame(t1 = n$df$satisfaction2, row.names = n$df$name) | |
| assign("satisfaction", satisfaction, envir = env) | |
| assign("satisfaction2", satisfaction2, envir = env) | |
| assign("weak_ties", n$weak_ties, envir = env) | |
| assign("strong_ties", n$strong_ties, envir = env) | |
| assign("leaderAdopt", makeCol(n$df$is_leader_obj * n$df$adopt)) | |
| assign("peerAdopt", makeCol((1-n$df$is_leader_obj) * n$df$adopt)) | |
| assign("activeAdopt", makeCol(isActive * n$df$adopt)) | |
| assign("inactiveAdopt", makeCol((1-isActive) * n$df$adopt)) | |
| environment(tnamdata) <- env | |
| df.tnam <- tnamdata(formula.tnam) | |
| df.tnam <- df.tnam[,-(1:2)] | |
| return(df.tnam) | |
| }) | |
| cols <- sapply(df.tnam, colnames) | |
| cols <- apply(cols, 1, function(i) rev(sort(unique(i)))[1]) | |
| df.tnam <- lapply(df.tnam, function(df) { | |
| colnames(df) <- cols | |
| return(df) | |
| }) | |
| df.tnam <- do.call(rbind, df.tnam) | |
| if(!is.null(rename)) colnames(df.tnam)[-1] <- rename | |
| df <- merge(df.glm, df.tnam) | |
| } else { | |
| df <- df.glm | |
| } | |
| nodes <- df$node | |
| df$node <- NULL | |
| formula <- as.formula(paste0(y, "~.")) | |
| mod <- glm(formula, data = df, ...) | |
| # mod <- lm(formula, data = df, ...) | |
| return(list(mod = mod, | |
| nodes = nodes, | |
| formula.glm = formula.glm, | |
| formula.tnam = formula.tnam)) | |
| } | |
| simIter <- function(mod, y, offset, coef, mnet, mdist) { | |
| adopt <- data.frame(t1 = y) | |
| env <- new.env() | |
| assign("uni", mnet, envir = env) | |
| assign("geo", mdist, envir = env) | |
| assign("adopt", adopt, envir = env) | |
| environment(tnamdata) <- env | |
| df.tnam <- tnamdata(mod$formula.tnam) | |
| df.tnam <- df.tnam[,-(1:2)] | |
| df.tnam$node <- NULL | |
| df.tnam <- as.matrix(df.tnam) | |
| pr <- as.numeric(1/(1+exp(-(offset + df.tnam %*% coef)))) | |
| y.new <- sapply(pr[!as.logical(y)], rbinom, size = 1, n = 1) | |
| y[!as.logical(y)] <- y.new | |
| return(list(y = y, pr = pr)) | |
| } | |
| simulate <- function(mod, offset, coef, mnet, mdist, i = NULL, maxit = 150, verbose = F, seeds = NULL) { | |
| nobs <- nrow(mnet) | |
| y <- rep(0,nobs) | |
| if(!is.null(seeds)) y <- seeds | |
| if(is.null(i)) i <- 1:nobs | |
| rownames(mnet) <- 1:nobs | |
| rownames(mdist) <- 1:nobs | |
| output <- data.frame(t = 0, | |
| i = i, | |
| y = y, | |
| pr = 0) | |
| t <- 0 | |
| while(sum(y) < nobs & t < maxit) { | |
| t <- t+1 | |
| it <- simIter(mod, y, offset, coef, mnet, mdist) | |
| output <- rbind(output, | |
| data.frame(t = t, | |
| i = i, | |
| y = it$y, | |
| pr = it$pr)) | |
| y <- it$y | |
| if (verbose) message("t = ", t, | |
| "; pct = ", round(sum(y)/nobs), | |
| "; Pr(adopt) = ", round(mean(it$pr[it$y == 0]), 4)) | |
| } | |
| return(output) | |
| } | |
| bootPredict <- function(formula, coefs, vcov, newdata, n) { | |
| boot <- mvrnorm(n, coefs, vcov) | |
| df <- model.matrix(formula, data = newdata) | |
| predict <- df %*% t(boot) | |
| predict <- 1/(1+exp(-predict)) | |
| return(predict) | |
| } | |
| toOffset <- function(x) -log(1/x - 1) | |
| adjust <- function(g1, g2) { | |
| if(vcount(g1) == vcount(g2)) { | |
| return(NULL) | |
| } else if (vcount(g1) < vcount(g2)) { | |
| gsmall <- g1 | |
| glarge <- g2 | |
| } else { | |
| gsmall <- g2 | |
| glarge <- g1 | |
| } | |
| diff <- vcount(glarge) - vcount(gsmall) | |
| deg <- table(degree(glarge)) | |
| min.sat <- as.numeric(names(deg)[max(which(cumsum(deg) <= diff))]) | |
| max.sat <- as.numeric(names(deg)[min(which(cumsum(deg) > diff))]) | |
| select <- V(glarge)[degree(glarge) <= min.sat] | |
| if (max.sat > min.sat) { | |
| select <- c(select, | |
| sample(V(glarge)[degree(glarge) == max.sat], diff - length(select))) | |
| } | |
| glarge <- glarge - select | |
| return(glarge) | |
| } | |
| getSeeds <- function(g, n, strat) { | |
| if (strat == "base") { | |
| y <- as.numeric(V(g) %in% sample(V(g),n)) | |
| } else if (strat == "cty") { | |
| y0 <- sample(V(g),1) | |
| y <- y0 | |
| sizes <- sapply(1:5, neighborhood.size, graph = g, nodes = y0) | |
| max.sat <- min(which(sizes > n)) | |
| min.sat <- max(which(sizes <= n)) | |
| if(min.sat != -Inf) y <- unique(c(y, ego(g,min.sat,y0)[[1]])) | |
| if (max.sat > min.sat) { | |
| target <- ego(g, max.sat, y0)[[1]] | |
| target <- target[!target %in% y] | |
| y <- c(y, sample(target, n-length(y))) | |
| y <- as.numeric(V(g) %in% y) | |
| } | |
| } else { | |
| weight <- degree(g, V(g))/sum(degree(g, V(g))) | |
| y <- as.numeric(V(g) %in% sample(V(g),n,prob = weight)) | |
| } | |
| return(y) | |
| } | |
| aggregateSims <- function(sims) { | |
| maxt <- max(sapply(sims, function(sim) sapply(sim, function(df) max(df$t)))) | |
| sims <- lapply(sims, function(sim) lapply(sim, function(df) { | |
| df <- aggregate(y ~ t + sim + type, df, mean) | |
| sim <- df$sim[1] | |
| type <- df$type[1] | |
| df <- merge(df, data.frame(t=1:maxt), all = T) | |
| df[is.na(df)] <- 1 | |
| df$sim <- sim | |
| df$type <- type | |
| return(df) | |
| })) | |
| sims <- lapply(sims, function(sim) do.call(rbind, sim)) | |
| sims <- do.call(rbind, sims) | |
| sims$sim <- paste0(sims$sim, sims$type) | |
| sims$sim <- match(sims$sim, unique(sims$sim)) | |
| return(sims) | |
| } | |
| topUp <- function(g, ediff, strat, cut = .1) { | |
| if (strat == "random") { | |
| select <- length(g[(!g) & lower.tri(g)]) | |
| select <- sample.int(n = select, size = ediff) | |
| g[(!g) & lower.tri(g)][select] <- 1 | |
| g <- t(g) | |
| g[(!g) & lower.tri(g)][select] <- 1 | |
| } else { | |
| diag(g) <- 1 | |
| n <- nrow(g) | |
| for (k in 1:ediff) { | |
| deg <- colSums(g) | |
| if (strat == "low") { | |
| ccut <- quantile(deg, cut) | |
| select <- which(deg <= ccut & deg < n) | |
| } else { | |
| ccut <- quantile(deg, 1-cut) | |
| select <- which(deg >= ccut & deg < n) | |
| } | |
| if (length(select) == 0) { | |
| stop("target is fully connected") | |
| } else if (length(select) == 1) { | |
| i <- select | |
| } else { | |
| i <- i <- sample(select, 1) | |
| } | |
| j <- sample(which(!g[i,]), 1) | |
| if (length(j) > 1) j <- sample(j, 1) | |
| g[i,j] <- g[j,i] <- 1 | |
| } | |
| diag(g) <- 0 | |
| } | |
| return(g) | |
| } | |
| fit2 <- function(formula.glm.hear, formula.tnam.hear = NULL, | |
| formula.glm.adopt, formula.tnam.adopt = NULL, | |
| rename.hear = NULL, rename.adopt = NULL, | |
| data, bayes = F, clogit = F, ...) { | |
| # if(!is.null(sub)) data <- data[sub] | |
| df.glm.hear <- lapply(data, function(l) { | |
| df.glm <- model.frame(formula.glm.hear, l$df) | |
| for (i in 1:length(colnames(df.glm))) { | |
| if(is.logical(df.glm[,i])) df.glm[,i] <- as.numeric(df.glm[,i]) | |
| } | |
| df.glm$node <- rownames(df.glm) | |
| return(df.glm) | |
| }) | |
| df.glm.hear <- do.call(rbind, df.glm.hear) | |
| df.glm.adopt <- lapply(data, function(l) { | |
| df.glm <- model.frame(formula.glm.adopt, l$df) | |
| for (i in 1:length(colnames(df.glm))) { | |
| if(is.logical(df.glm[,i])) df.glm[,i] <- as.numeric(df.glm[,i]) | |
| } | |
| df.glm$node <- rownames(df.glm) | |
| return(df.glm) | |
| }) | |
| df.glm.adopt <- do.call(rbind, df.glm.adopt) | |
| if(!is.null(formula.tnam.hear)) { | |
| df.tnam.hear <- lapply(data, function(n) { | |
| env <- new.env() | |
| assign("uni", n$union, envir = env) | |
| assign("solver", n$solver, envir = env) | |
| assign("lender", n$lender, envir = env) | |
| assign("family", n$family, envir = env) | |
| assign("friend", n$friend, envir = env) | |
| assign("leader", n$leader, envir = env) | |
| assign("leader_obj", n$leader_obj, envir = env) | |
| assign("geo", n$geo, envir = env) | |
| assign("heard", n$heard, envir = env) | |
| assign("adopt", n$adopt, envir = env) | |
| assign("is_leader", n$df$is_leader, envir = env) | |
| assign("is_leader_obj", n$df$is_leader_obj, envir = env) | |
| assign("weak_ties", n$weak_ties, envir = env) | |
| assign("strong_ties", n$strong_ties, envir = env) | |
| assign("speak", n$speak, envir = env) | |
| satisfaction <- data.frame(t1 = n$df$satisfaction, row.names = n$df$name) | |
| satisfaction2 <- data.frame(t1 = n$df$satisfaction2, row.names = n$df$name) | |
| assign("satisfaction", satisfaction, envir = env) | |
| assign("satisfaction2", satisfaction2, envir = env) | |
| environment(tnamdata) <- env | |
| df.tnam <- tnamdata(formula.tnam.hear) | |
| df.tnam <- df.tnam[,-(1:2)] | |
| return(df.tnam) | |
| }) | |
| cols <- sapply(df.tnam.hear, colnames) | |
| cols <- apply(cols, 1, function(i) rev(sort(unique(i)))[1]) | |
| df.tnam.hear <- lapply(df.tnam.hear, function(df) { | |
| colnames(df) <- cols | |
| return(df) | |
| }) | |
| df.tnam.hear <- do.call(rbind, df.tnam.hear) | |
| if(!is.null(rename.hear)) colnames(df.tnam.hear)[-1] <- rename.hear | |
| df.hear <- merge(df.glm.hear, df.tnam.hear) | |
| df.hear <- df.hear[, | |
| c(colnames(df.tnam.hear)[-which(colnames(df.tnam.hear) == "node")], | |
| colnames(df.glm.hear))] | |
| } else { | |
| df.hear <- df.glm.hear | |
| } | |
| if(!is.null(formula.tnam.adopt)) { | |
| df.tnam.adopt <- lapply(data, function(n) { | |
| env <- new.env() | |
| assign("uni", n$union, envir = env) | |
| assign("solver", n$solver, envir = env) | |
| assign("lender", n$lender, envir = env) | |
| assign("family", n$family, envir = env) | |
| assign("friend", n$friend, envir = env) | |
| assign("leader", n$leader, envir = env) | |
| assign("leader_obj", n$leader_obj, envir = env) | |
| assign("is_leader", n$df$is_leader, envir = env) | |
| assign("is_leader_obj", n$df$is_leader_obj, envir = env) | |
| assign("geo", n$geo, envir = env) | |
| assign("heard", n$heard, envir = env) | |
| assign("adopt", n$adopt, envir = env) | |
| assign("speak", n$speak, envir = env) | |
| satisfaction <- data.frame(t1 = n$df$satisfaction, row.names = n$df$name) | |
| satisfaction2 <- data.frame(t1 = n$df$satisfaction2, row.names = n$df$name) | |
| assign("satisfaction", satisfaction, envir = env) | |
| assign("satisfaction2", satisfaction2, envir = env) | |
| assign("weak_ties", n$weak_ties, envir = env) | |
| assign("strong_ties", n$strong_ties, envir = env) | |
| environment(tnamdata) <- env | |
| df.tnam <- tnamdata(formula.tnam.adopt) | |
| df.tnam <- df.tnam[,-(1:2)] | |
| return(df.tnam) | |
| }) | |
| cols <- sapply(df.tnam.adopt, colnames) | |
| cols <- apply(cols, 1, function(i) rev(sort(unique(i)))[1]) | |
| df.tnam.adopt <- lapply(df.tnam.adopt, function(df) { | |
| colnames(df) <- cols | |
| return(df) | |
| }) | |
| df.tnam.adopt <- do.call(rbind, df.tnam.adopt) | |
| if(!is.null(rename.adopt)) colnames(df.tnam.adopt)[-1] <- rename.adopt | |
| df.adopt <- merge(df.glm.adopt, df.tnam.adopt) | |
| df.adopt <- df.adopt[, | |
| c(colnames(df.tnam.adopt)[-which(colnames(df.tnam.adopt) == "node")], | |
| colnames(df.glm.adopt))] | |
| } else { | |
| df.adopt <- df.glm.adopt | |
| } | |
| # nodes <- intersect(df.hear$node, df.adopt$node) | |
| # df.hear <- df.hear[df.hear$node %in% nodes,] | |
| # df.adopt <- df.adopt[df.adopt$node %in% df$nodes,] | |
| df.adopt.full <- df.adopt[df.adopt$node %in% df.hear$node,] | |
| df.adopt <- df.adopt[df.adopt$node %in% df.hear$node[df.hear$hear == 1],] | |
| nodes <- df.hear$node | |
| nodesa <- df.adopt$node | |
| df.hear$node <- NULL | |
| df.adopt$node <- NULL | |
| df.adopt.full$node <- NULL | |
| if (bayes) { | |
| mod.hear <- bayesglm(heard ~ ., data = df.hear, ...) | |
| mod.adopt <- bayesglm(adopt ~ ., data = df.adopt, ...) | |
| } else if (clogit) { | |
| mod.hear <- clogit(heard ~ . - village + strata(village), data = df.hear, model = TRUE, y = TRUE, ...) | |
| mod.adopt <- clogit(adopt ~ . - village + strata(village), data = df.adopt, model = TRUE, y = TRUE, ...) | |
| mod.hear$AIC <- AIC(mod.hear) | |
| mod.adopt$AIC <- AIC(mod.adopt) | |
| } else { | |
| mod.hear <- glm(heard ~ ., data = df.hear, ...) | |
| mod.adopt <- glm(adopt ~ ., data = df.adopt, ...) | |
| } | |
| obj <- list(nodes = nodes, | |
| nodesa = nodesa, | |
| mod.hear = mod.hear, | |
| mod.adopt = mod.adopt, | |
| df.adopt.full = df.adopt.full, | |
| formula.glm.hear = formula.glm.hear, | |
| formula.glm.adopt = formula.glm.adopt, | |
| formula.tnam.hear = formula.tnam.hear, | |
| formula.tnam.adopt = formula.tnam.adopt) | |
| if(clogit) obj$loglik <- mod.hear$loglik[2] + mod.adopt$loglik[2] | |
| obj$AIC <- 2 * (length(coef(obj$mod.hear)) + length(coef(obj$mod.adopt)) - obj$loglik) | |
| class(obj) <- "mfit2" | |
| return(obj) | |
| } | |
| predict.mfit2 <- function(obj, dfHear = NULL, dfAdopt = NULL) { | |
| if(is.null(dfHear) || is.null(dfAdopt)) { | |
| return(predict(obj$mod.hear, type = "response") * | |
| predict(obj$mod.adopt, newdata = obj$df.adopt.full, type = "response")) | |
| } | |
| return(predict(obj$mod.hear, newdata = dfHear, type = "response") * | |
| predict(obj$mod.adopt, newdata = dfAdopt, type = "response")) | |
| } | |
| summary.mfit2 <- function(obj) summary(obj$mod.adopt) | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment