# code to fit the GLMM model p2 to the Caltech data # Note that this is an inappropriate model as the network is undirected so we drop the reciprocity term library(lme4) Y<-as.matrix(read.table("CaltechY.dat")) N<-nrow(Y) i<-c(t(matrix(rep(1:N,N),N))) j<-rep(1:N,N) df<-data.frame(y=c(t(Y)),i=as.factor(i),j=as.factor(j)) # read in the covariate information localinfo<-read.csv("Caltech_localinfo.csv") df$status<-localinfo$status[i]!=localinfo$status[j] df$gender<-localinfo$gender[i]!=localinfo$gender[j] df$major<-localinfo$major[i]!=localinfo$major[j] df$minor<-localinfo$minor[i]!=localinfo$minor[j] df$dorm<-localinfo$dorm[i]!=localinfo$dorm[j] df$year<-localinfo$year[i]!=localinfo$year[j] df$highschool<-localinfo$highschool[i]!=localinfo$highschool[j] m<-lmer(data=df, y~1+status+gender+major+minor+dorm+year+highschool+(1|i)+(1|j),family=binomial) # now look at model fit using a ROC curve library(boot) probs<-matrix(inv.logit(fixef(m)[[1]]+fixef(m)[[2]]*df$status+fixef(m)[[3]]*df$gender+fixef(m)[[4]]*df$major+ fixef(m)[[5]]*df$minor+fixef(m)[[6]]*df$dorm+fixef(m)[[7]]*df$year+fixef(m)[[8]]*df$highschool+ ranef(m)[[1]][i,1]+ranef(m)[[2]][j,1]),N) source("roc.R") roc(probs,Y)