library(lme4)
options(rgl.useNULL=TRUE)
library(LMERConvenienceFunctions)

library(ggplot2)
options(jupyter.plot_mimetypes = 'image/png')
library(repr)
options(repr.plot.width=6, repr.plot.height=6)

#no zeros in the outcome variable
data = read.csv("/home/psinger/tmp/data2.csv", header=TRUE)

nrow(data)

m_lmer = lmer(body_length~1+index+(1|author)+(1|length), data=data, REML=FALSE)

summary(m_lmer)

mcp.fnc(m_lmer)

data_limit <- romr.fnc(m_lmer, data, trim = 2.5)$data

m_lmer_limit = lmer(body_length~1+index+(1|author)+(1|length), data=data_limit, REML=FALSE)

summary(m_lmer_limit)

mcp.fnc(m_lmer_limit)

m_lmer_log = lmer(log(body_length)~1+index+(1|author)+(1|length), data=data, REML=FALSE)

summary(m_lmer_log)

mcp.fnc(m_lmer_log)

m_lognormal = glmer(body_length~1+index+(1|author)+(1|length),data=data,family=gaussian(link=log))

summary(m_lognormal)

plot(resid(m_lognormal,type="deviance") ~ fitted(m_lognormal), pch = ".")

qqnorm(resid(m_lognormal,type="deviance"), pch = ".", main = "")
qqline(resid(m_lognormal,type="deviance"))

m_poisson = glmer(body_length~1+index+(1|author)+(1|length),data=data,family=poisson)

summary(m_poisson)

overdisp_fun <- function(model) {
  ## number of variance parameters in 
  ##   an n-by-n variance-covariance matrix
  vpars <- function(m) {
    nrow(m)*(nrow(m)+1)/2
  }
  model.df <- sum(sapply(VarCorr(model),vpars))+length(fixef(model))
  rdf <- nrow(model.frame(model))-model.df
  rp <- residuals(model,type="pearson")
  Pearson.chisq <- sum(rp^2)
  prat <- Pearson.chisq/rdf
  pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
  c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)
}

overdisp_fun(m_poisson)

plot(scale(resid(m_poisson,type="deviance")) ~ fitted(m_poisson), pch = ".")

data$obs<-seq.int(nrow(data))

m_poisson2 = glmer(body_length~1+index+(1|author)+(1|length)+(1|obs),data=data,family=poisson)

m_poisson2 = glmer(body_length~1+index+length+(1|author)+(1|obs),data=data,family=poisson)

m_poisson2 = glmer(body_length~1+index+length+(1|author)+(1|obs),data=data[1:10000,],family=poisson)

summary(m_poisson2)

overdisp_fun(m_poisson2)

plot(resid(m_poisson2,type="deviance") ~ fitted(m_poisson2), pch = ".")

qqnorm(resid(m_poisson2,type="deviance"), pch = ".", main = "")
qqline(resid(m_poisson2,type="deviance"))

m_nbinom = glmer.nb(body_length~1+index+length+(1|author),data=data[1:10000,])

summary(m_nbinom)

plot(resid(m_nbinom,type="deviance") ~ fitted(m_nbinom), pch = ".")

qqnorm(resid(m_nbinom,type="deviance"), pch = ".", main = "")
qqline(resid(m_nbinom,type="deviance"))