email3b<-function(n,lambda.WD,lambda.WE,b1,b2,n.days){
  #vertex set (n people)
  V=1:n
  #starting time of simulation (a Monday)
  tm<-0
  #initialize 
  tms<-NULL
  sender<-NULL
  receiver<-NULL
  non.sender<-NULL
  non.receiver<-NULL
  rep<-matrix(-Inf,ncol=n,nrow=n)
  rec<-matrix(-Inf,ncol=n,nrow=n)
  repetition<-NULL
  reciprocity<-NULL
  non.repetition<-NULL
  non.reciprocity<-NULL
  while (tm<n.days){
    # decide which lambda to use
    day<-tm%%7
    if (day<5){
      lambda0<-lambda.WD
      WD<-TRUE
    } else {
      lambda0<-lambda.WE
      WD<-FALSE
    }
    
    lambda<-lambda0*exp(b1*(tm-rec<2)+b2*(tm-rep<2))
    diag(lambda)<-0
    
    delta.tm<-rexp(1,sum(lambda))
    
    if ((WD & (day+delta.tm<5))|(!WD & (day+delta.tm<7))){
      #sample event
      event<-sample(1:n^2,1,prob=lambda/sum(lambda))
      snd<-(event-1)%%n+1
      rcv<-(event-1)%/%n+1
      #sample non-event
      non.event<-sample(setdiff(1:n^2,event),1)
      non.snd<-(non.event-1)%%n+1
      non.rcv<-(non.event-1)%/%n+1
      #determine whether event was repetition/reciprocal
      repetition<-c(repetition,1*(tm-rep[snd,rcv]<2))
      reciprocity<-c(reciprocity,1*(tm-rec[snd,rcv]<2))
      #determine whether event was repetition/reciprocal
      non.repetition<-c(non.repetition,1*(tm-rep[non.snd,non.rcv]<2))
      non.reciprocity<-c(non.reciprocity,1*(tm-rec[non.snd,non.rcv]<2))
      #save relational event/non-event
      tm<-tm+delta.tm
      tms<-c(tms,tm)
      sender<-c(sender,snd)
      receiver<-c(receiver,rcv)
      non.sender<-c(non.sender,non.snd)
      non.receiver<-c(non.receiver,non.rcv)
      #update reciprocity and repetition
      rep[snd,rcv]<-tm
      rec[rcv,snd]<-tm
    } else {
      if (WD){
        tm<-(tm%/%7)*7+5
      } else {
        tm<-(tm%/%7)*7+7
      }
    }
  }
  remdat<-data.frame(tms=tms,sender=sender,receiver=receiver,
      non.sender=non.sender,non.receiver=non.receiver,
      repetition=repetition,reciprocity=reciprocity,
      non.repetition=non.repetition,non.reciprocity=non.reciprocity)
  return(remdat)
}


set.seed(3)
remdat<-email3b(10,.05,0.02,1,1.5,365)


# Let's inspect the data
hist(remdat$tms,breaks=365)

# Fitting REM
delta.rec<-remdat$reciprocity-remdat$non.reciprocity
delta.rep<-remdat$repetition-remdat$non.repetition
ones<-rep(1,length(delta.rep))
email3b.gam<-gam(ones~-1+delta.rec+delta.rep,family = binomial)
summary(email3b.gam)




#####################################
# WTC preprocessing

wtc<-read.csv2("WTC.csv")
head(wtc)
attach(wtc)


# Sample non-events
set.seed(1)
V<-unique(c(s,r))
ICR<-NULL
for (v in V){
  v.ICR<-c(s.ICR[s==v],r.ICR[r==v])
  v.ICR<-v.ICR[!is.na(v.ICR)]
  ICR<-c(ICR,v.ICR[1])
}

n<-nrow(wtc)
non.s<-NULL
non.r<-NULL
non.s.ICR<-NULL
non.r.ICR<-NULL
for (i in 1:n){
  same=TRUE
  while (same){
    non.si<-sample(V,1)
    non.ri<-sample(setdiff(V,non.si),1)
    same <- (s[i]==non.si)&(r[i]==non.ri)
  }
  non.s<-c(non.s,non.si)
  non.r<-c(non.r,non.ri)
  non.s.ICR<-c(non.s.ICR,ICR[V==non.si])
  non.r.ICR<-c(non.r.ICR,ICR[V==non.ri])
}

# Determine repetition, reciprocity 
repetition<-0
reciprocity<-0
non.repetition<-0
non.reciprocity<-0

for (i in 2:n){
  #repetition
  repi<-sum((s[i]==s[1:(i-1)])&(r[i]==r[1:(i-1)]))
  repetition<-c(repetition,repi)
  non.repi<-sum((non.s[i]==s[1:(i-1)])&(non.r[i]==r[1:(i-1)]))
  non.repetition<-c(non.repetition,non.repi)
  #reciprocity
  reci<-sum((s[i]==r[1:(i-1)])&(r[i]==s[1:(i-1)]))
  reciprocity<-c(reciprocity,reci)
  non.reci<-sum((non.s[i]==r[1:(i-1)])&(non.r[i]==s[1:(i-1)]))
  non.reciprocity<-c(non.reciprocity,non.reci)
}

remdat<-cbind(wtc,non.s,non.r,non.s.ICR,non.r.ICR,repetition,reciprocity,non.repetition,non.reciprocity)
write.csv2(remdat,"WTC-rem.csv",row.names = FALSE)
remdat<-read.csv2("WTC-rem.csv")



#################################################
# WTC analysis

attach(remdat)
delta.s.ICR<-s.ICR-non.s.ICR
delta.r.ICR<-r.ICR-non.r.ICR
delta.rec<-reciprocity-non.reciprocity
delta.rep<-repetition-non.repetition
ones<-rep(1,length(delta.r.ICR))
wtc1.rem<-gam(ones~-1+delta.s.ICR+delta.r.ICR+delta.rec+delta.rep, family=binomial)
summary(wtc1.rem)
wtc1.rem$aic

wtc2.rem<-gam(ones~-1+delta.s.ICR+delta.r.ICR+delta.rec, family=binomial)
summary(wtc2.rem)
wtc2.rem$aic




##########################################################
# Class data preprocessing

cls<-read.csv2("classroom.csv")
head(cls)
attach(cls)


# Sample non-events
set.seed(1)
V<-unique(c(s,r))
teacher<-NULL
female<-NULL
for (v in V){
  v.teacher<-c(s.teacher[s==v],r.teacher[r==v])
  v.teacher<-v.teacher[!is.na(v.teacher)]
  teacher<-c(teacher,v.teacher[1])
  v.female<-c(s.female[s==v],r.female[r==v])
  v.female<-v.female[!is.na(v.female)]
  female<-c(female,v.female[1])
}

n<-nrow(cls)
non.s<-NULL
non.r<-NULL
non.s.teacher<-NULL
non.r.teacher<-NULL
non.s.female<-NULL
non.r.female<-NULL
for (i in 1:n){
  same=TRUE
  while (same){
    non.si<-sample(V,1)
    non.ri<-sample(setdiff(V,non.si),1)
    same <- (s[i]==non.si)&(r[i]==non.ri)
  }
  non.s<-c(non.s,non.si)
  non.r<-c(non.r,non.ri)
  non.s.teacher<-c(non.s.teacher,teacher[V==non.si])
  non.r.teacher<-c(non.r.teacher,teacher[V==non.ri])
  non.s.female<-c(non.s.female,female[V==non.si])
  non.r.female<-c(non.r.female,female[V==non.ri])
}

# Determine repetition, reciprocity 
repetition<-0
reciprocity<-0
non.repetition<-0
non.reciprocity<-0

for (i in 2:n){
  #repetition
  repi<-sum((s[i]==s[1:(i-1)])&(r[i]==r[1:(i-1)]))
  repetition<-c(repetition,repi)
  non.repi<-sum((non.s[i]==s[1:(i-1)])&(non.r[i]==r[1:(i-1)]))
  non.repetition<-c(non.repetition,non.repi)
  #reciprocity
  reci<-sum((s[i]==r[1:(i-1)])&(r[i]==s[1:(i-1)]))
  reciprocity<-c(reciprocity,reci)
  non.reci<-sum((non.s[i]==r[1:(i-1)])&(non.r[i]==s[1:(i-1)]))
  non.reciprocity<-c(non.reciprocity,non.reci)
}

remdat<-cbind(cls,non.s,non.r,non.s.teacher,non.r.teacher,non.s.female,non.r.female,repetition,reciprocity,non.repetition,non.reciprocity)
write.csv2(remdat,"classroom-rem.csv",row.names = FALSE)
remdat<-read.csv2("classroom-rem.csv")



#################################################
# classroom analysis

attach(remdat)
delta.s.teacher<-s.teacher-non.s.teacher
delta.r.teacher<-r.teacher-non.r.teacher
delta.samerole<- (s.teacher==r.teacher)-(non.s.teacher==non.r.teacher)
delta.s.female<-s.female-non.s.female
delta.r.female<-r.female-non.r.female
delta.samegender<- (s.female==r.female)-(non.s.female==non.r.female)
delta.rec<-reciprocity-non.reciprocity
delta.rep<-repetition-non.repetition
ones<-rep(1,length(delta.r.teacher))
class1.rem<-gam(ones~-1+delta.s.teacher+delta.r.teacher+delta.samerole+delta.s.female+delta.r.female+delta.samegender+delta.rec+delta.rep, family=binomial)
summary(class1.rem)
class1.rem$aic

class2.rem<-gam(ones~-1+delta.s.teacher+delta.r.teacher+delta.r.female+delta.rec+delta.rep, family=binomial)
summary(class2.rem)
class2.rem$aic


delta.dat<-data.frame(ones,delta.s.teacher,delta.r.teacher,delta.samerole,delta.s.female,delta.r.female,delta.samegender,delta.rep,delta.rec)
class.rem<-glm(ones~-1++delta.s.teacher+delta.r.teacher+delta.samerole+delta.s.female+delta.r.female+delta.samegender+delta.rec+delta.rep, family=binomial,dat=delta.dat)
library(MASS)
stepAIC(class.rem)