#############
#PROGRAM INFO
#############

# Noel and Nyhan simulation code for "The 'Unfriending' Problem: The Consequences of Homophily in Friendship Retention for Causal Estimates of Social Influence," which is available at http://www-personal.umich.edu/~bnyhan/unfriending.pdf

# Based on Christakis and Fowler's simulation code for "Estimating Peer Effects on Health in Social Networks" (JHE, 2008), which is available at http://jhfowler.ucsd.edu/homophily_and_influence_monte_carlo.R

# Note: We have added the prefix "CF:" to comments from Christakis and Fowler's original code. All other comments are our own (we have added "NN:" if there is ambiguity).

##############
#BEGIN PROGRAM
##############

library(geepack) # CF: load library to do GEE estimation

nsim <- 1000 # number of simulations for each unique combination of parameters

# make array of unique combinations of simulation parameters

betacoef0seq<-c(0,.0125,.025,.0375,.05) #homophily coefficient at time 0
betacoef1seq<-c(0,.025,.05) #homophily coefficient at time 1
constant1seq<-c(0,.5,1,1.85) #baseline attrition at time 1
sdseq<-c(5) #we hold s.d. of shock fixed at 5
truebseq<-c(0) #CF: true influence coefficient (NN: fixed at 0)
nseq<-c(1000) #number of individuals

Param<-list(betacoef0=betacoef0seq,betacoef1=betacoef1seq,constant1=constant1seq,shocksd=sdseq,b1=truebseq,n=nseq)
CombParam <- expand.grid(Param)

trueb<-.0 #true influence coefficient

# create make/break function

makebreak<-function(betacoef0,betacoef1,constant1,shocksd,b1,n) {

   # CF: fix vector of a random normal feeling data  
   # CF: for time period 0
   y0<-rnorm(n=n,mean=50,sd=10)
  
   # diff in y0 (abs. value positive)
   distance0 <- -abs(outer(y0,y0,"-"))
   # make mean 0
   distance0 <- distance0 - mean(as.vector(distance0))
   # eps~N(0,1) for latent variable probit
   epsilon0<-rnorm(n^2,0,1)
   constant0 <- -2.5
   ystar0 <- constant0 + (betacoef0*distance0)        
   # CF: realize a random draw for each probability of a tie
   tie<-ystar0>epsilon0
   # CF: set probability of self-tie to 0
   diag(tie)<-0
   tie<-tie>0
 
   # CF: create a list of all ties
   tie_index<-which(tie)
   # CF: set non-ties to missing values
   tie[!tie]<-NA
   # CF: create index values for the senders of the ties (the 'namers')
   namer<-trunc((tie_index-1)/n)+1
   # CF: create index values for the receivers of the ties (the 'named')
   named<-tie_index%%n 
   named[named==0]<-n

   # correlation in y0 among friends at t0 
   r.form <- cor(y0[namer],y0[named])
   
   # average number of friends per individual
   numt0 <- length(named)/n

   # CF: let new feeling equal old feeling plus a shock that is a random normal    
   #     change in feeling
   y1<-y0+rnorm(n,0,shocksd)
   # CF: find mean feeling of all alters for each ego
   y1_alters<-rowMeans(t(t(tie)*y1),na.rm=T)
   # CF: replace mean feeling of egos without alters so that it equals the ego 
   #     feeling
   y1_alters[is.na(y1_alters)]<-y1[is.na(y1_alters)]
   # CF: let ego feeling be influenced by own feeling and by alter feeling
   y1<-b1*y1_alters+(1-b1)*y1

   # prob of staying friends
   distance1 <- -abs(y1[named]-y1[namer])
   #make mean 0
   distance1 <- distance1 - mean(as.vector(distance1))
   # eps~N(0,1) for latent variable probit
   ystar1 <- constant1 + (betacoef1*distance1)        
   epsilon1<-rnorm(length(ystar1),0,1)
   # realize a random draw for each probability of a tie
   observed<-ystar1>epsilon1
     
   # create new vectors of namers and named
   named1 <- named[observed==1]
   namer1 <- namer[observed==1]

   # correlation in y0 among friends at t1
   r.break <- cor(y1[namer],y1[named])
   
   # average number of friends per individual
   numt1<-length(named1)/n
        
   # GEE model of influence per CF
   g<-summary(geeglm(y1[namer1]~y0[namer1]+y1[named1]+y0[named1],id=namer1))

   # return estimate of influence coefficients, s.e., p-values and other 
   # quantities of interest
   cbind(g$coefficients[3,][1],g$coefficients[3,][2],g$coefficients[3,][4]<.05,
   r.form,r.break,numt0,numt1)
}

#function that repeats makebreak nsim times

makebreakrep<-function(betacoef0,betacoef1,constant1,shocksd,b1,n){
	
sims<-replicate(nsim,makebreak (betacoef0=betacoef0,betacoef1=betacoef1,constant1=constant1,n=n,b1=b1,shocksd=shocksd),simplify=TRUE)
 
   #record mean values for parameters of interest over 500 simulations      
   meanbias<-mean(unlist(sims[1,]))
   meanse<-mean(unlist(sims[2,]))
   rejection<-mean(unlist(as.numeric(sims[3,])))
   coverage<-1-rejection
   meanrform <- mean(unlist(sims[4,]))
   meanrbreak <- mean(unlist(sims[5,]))
   meannumt0 <- mean(unlist(sims[6,]))
   meannumt1 <- mean(unlist(sims[7,]))
   retention <- meannumt1/meannumt0
	
   rbind(meanbias,meanse,coverage,meanrform,meanrbreak,meannumt0,meannumt1,retention)
}

#use apply to run makebreakrep over CombParam array
MonteCarlo <- apply(CombParam, MARGIN = 1,
function(x) makebreakrep(x[1],x[2],x[3],x[4],x[5],x[6]))

#create final simulation results
options(digits=3)
options(scipen=999)
MCdata<-cbind(CombParam[,3],CombParam[,1],CombParam[,2],t(MonteCarlo)[,1],t(MonteCarlo)[,3:8])
MCdata<-MCdata[order(MCdata[,1],MCdata[,2]),]
MCdata