sse1=sse2=sse3=rep(NA,nsims)

for(i in 1:nsims)

{

reorder=sample(n)

train=data1[reorder[1:2500],]

test=data1[reorder[2500:n],]

###quasi_pois: result of backward model on Quasi Poisson

quasi_pois=glm(visits ~ stem + word_count + user_generated, family =

"quasipoisson",data=train)

###neg_binom: result of step function on Negative Binomial

neg_binom=glm.nb(visits ~ stem + upload_year_2013 + word_count +

user_generated + diff_level_2 + diff_level_3, data = train,

init.theta = 0.263265724, link = log)

###zero_infl: result of step function on zero inflated

zero_infl=zeroinfl(visits~stem + other_dept+upload_year_2013 +

word_count+doc_title_length+course_title_length +

user_generated+diff_level_2+diff_level_3 | stem+other_dept +

upload_year_2013 + word_count+doc_title_length + course_title_length

+ user_generated+diff_level_2 + diff_level_3,data=train)

# % of zeros in model calculation

zpoi[i]=mean(exp(-exp(predict(quasi_pois,new=test))))

z_nb[i]=mean((neg_binom$theta/(exp(predict(neg_binom,new=test))

+ neg_binom$theta))^neg_binom$theta)

#zip[i]=mean(predict(zero_infl,type="zero",new=test) + (1-

predict( zero_infl,type="zero",new=test))*exp(-

predict(zero_infl,type="count",new=test))

# deviance calculation

devia_qpoi[i]=2*sum(exp(predict(quasi_pois,new=test)) *

log(exp(predict(quasi_pois,new=test))/

mean(exp(predict(quasi_pois,new=test)))) -

(exp(predict(quasi_pois,new=test))-

mean(exp(predict(quasi_pois,new=test)))))

devia_nb[i]=deviance(neg_binom) # the computation by formula is #

# beyond the time scope of this project for us. To show that we have

# researched it I am pasting the link of the formula page here

# http://www.mathematica-journal.com/2013/06/negative-binomial-