library(lme4)
covs<-read.csv("LazegaLawyers/ELattr.csv")
Y<-as.matrix(read.table("LazegaLawyers/ELfriend.dat"))
advice<-as.matrix(read.table("LazegaLawyers/ELadv.dat"))
cowork<-as.matrix(read.table("LazegaLawyers/ELwork.dat"))
# only look at associates in the firm
Y<-Y[covs$status==2,covs$status==2]
advice<-advice[covs$status==2,covs$status==2]
cowork<-cowork[covs$status==2,covs$status==2]
covs<-covs[covs$status==2,]
N=nrow(Y)
# sender random effects
i<-as.factor(c(t(matrix(rep(1:N,N),N))))
# receiver random effects
j<-as.factor(rep(1:N,N))
diag(Y)<-NaN
M<-Y # using ergmMPLE(network(Y) ~ mutual) and checking M (predictor) and Y (response) values shows this is right
df<-data.frame(y=c(t(Y)),M=c(M),i=as.factor(i),j=as.factor(j))
# we want to predict Y[j,i] from Y[i,j]
df$location<-covs$office[i]!=covs$office[j]
df$aseniority<-abs(covs$seniority[i]-covs$seniority[j])
df$seniority<-covs$seniority[i]-covs$seniority[j]
df$gender<-covs$gender[i]!=covs$gender[j]
df$speciality<-covs$practice[i]!=covs$practice[j]
df$status<-covs$status[i]+covs$status[j]
df$years<-covs$years[i]+covs$years[j]
df$age<-covs$age[i]+covs$age[j]
df$school<-covs$school[i]+covs$school[j]

# same models as original p2 paper
m0<-lmer(data=df, y~1+M+(1|i)+(1|j),family=binomial)
m1<-lmer(data=df, y~1+location+aseniority+seniority+gender+speciality+M+(0+M|location)+(1|i)+(1|j),family=binomial)
m2<-lmer(data=df, y~1+location+aseniority+gender+c(advice)+c(cowork)+M+(1|i)+(1|j),family=binomial)

# results from original p2 paper
tmp1<-c(-0.64, -2.33, -0.58, 0.18, -0.55, -0.51, 2.21, 1.72, 1.19, 0.63, -0.01)
cor(tmp1,c(fixef(m1),sapply(ranef(m1),var),cor(ranef(m1)[[1]],ranef(m1)[[2]])))
tmp2<-c(-1.62,-1.45,-0.49,-0.58,1.5,0.53,2.22,1.36,0.65, 0.16)
cor(tmp2,c(fixef(m2),sapply(ranef(m2),var),cor(ranef(m2)[[1]],ranef(m2)[[2]])))