getwd()
install.packages("readxl")
library(readxl)

getwd()

data<-read_excel("Article Misuse_Results.xlsx")

head(data,n=10)

class(data)


data<-as.data.frame(data)


data$Essay<-substring(data$Essay, 5, 8) #extract four digit numbers

data$CorrectnessOmission<-factor(data$CorrectnessOmission,levels=c(0,1),
                                 labels = c("No","Yes"))
data$PreMod<-factor(data$PreMod,levels=c(0,1),labels=c("Non","Pre"))

data$PreModPOS<-factor(data$PreModPOS,levels=c(1,2),labels=c("Adj","Noun"))

data$Animacy<-factor(data$Animacy,levels=c(2,1),labels=c("Inanimate","Animate"))

data$Number<-factor(data$Number,levels=c(2,1),labels=c("Plural","Singular"))

str(data)

summary(data)

data2<-data[complete.cases(data[,c("CorrectnessOmission", "PreMod", "Animacy")]),] #optional


head(data2) #optional

str(data2) #optional

install.packages("lme4")

library(lme4)

# fitting generalized linear mixed effects models with logit link


with(data,table(CorrectnessOmission,PreMod))

with(data,table(CorrectnessOmission,PreModPOS))

with(data,table(CorrectnessOmission,Animacy))

with(data,table(CorrectnessOmission,Number))


#Univariate
model1<-glmer(CorrectnessOmission ~ PreMod + (1|Essay),data=data,
             family = binomial(link="logit"),
             control=glmerControl(optimizer = "Nelder_Mead"))

summary(model1)


model2<-glmer(CorrectnessOmission ~ PreModPOS + (1|Essay),data=data,
              family = binomial(link="logit"),
              control=glmerControl(optimizer = "Nelder_Mead"))

summary(model2)


model3<-glmer(CorrectnessOmission ~ Animacy + (1|Essay),data=data,
              family = binomial(link="logit"),
              control=glmerControl(optimizer = "Nelder_Mead"))

summary(model3)

model4<-glmer(CorrectnessOmission ~ Number + (1|Essay),data=data,
              family = binomial(link="logit"),
              control=glmerControl(optimizer = "Nelder_Mead"))

summary(model4)

#Multiple predictors

model5<-glmer(CorrectnessOmission ~ PreMod + Animacy + Number +(1|Essay),data=data,
              family = binomial(link="logit"),
              control=glmerControl(optimizer = "Nelder_Mead"))

summary(model5)

formula(model1)
formula(model2)
formula(model3)
formula(model4)
formula(model5)

model6<-glmer(CorrectnessOmission ~ Animacy + Number + Animacy*Number+(1|Essay),data=data,
              family = binomial(link="logit"),
              control=glmerControl(optimizer = "Nelder_Mead"))

summary(model5)

#optional

anova(model1, model3, type="LRT") #log likelyhood ratio test)

#p-value is less than alpha = 0.05, then choose the model with small AIC
#p-value is greater than alpha = 0.05, then choose the model with small df.

#AIC, BIC< deviance ->goodness of fit test statistics->smaller is better

####################################################
install.packages("sjPlot")

library(sjPlot)

tab_model(model1, show.r2=F, show.aic=T,show.dev=T)

tab_model(model2, show.r2=F, show.aic=T,show.dev=T)

tab_model(model3, show.r2=F, show.aic=T,show.dev=T)

tab_model(model4, show.r2=F, show.aic=T,show.dev=T)

tab_model(model5, show.r2=F, show.aic=T,show.dev=T)

tab_model(model6, show.r2=F, show.aic=T,show.dev=T)