---
title: "activity context teaching practices"
author: "Matteo"
date: "19/01/2021"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)



library(stats)
library(readr)
library(tidyverse)
library(MASS)
library(sjPlot)
library(DCluster)
library(pscl)
library(rcompanion)

#dataset
setwd("C:/Users/Utente/Documents/Matteo PhD/study3")
data<-read_csv("C:/Users/Utente/Documents/Matteo PhD/study3/TeachingFINALdata.csv")
#to check variance
data1<-data%>%
  filter(group %in% c(1))
data2<-data%>%
  filter(group %in% c(2))
data3<-data%>%
  filter(group %in% c(3))

data$group<-as.factor(data$group)
data$group
data$teacher<-as.factor(data$teacher)
data$teacher

```




```{r cars, message=FALSE, warning=FALSE, paged.print=FALSE, include=FALSE}


#IndividualActivity Poisson / Negative Binomial
mean(data1$IndividualActivity)
var(data1$IndividualActivity)
mean(data2$IndividualActivity)
var(data2$IndividualActivity)
mean(data3$IndividualActivity)
var(data3$IndividualActivity)

IndividualActivity1<-glm(IndividualActivity~ group+offset(log(Countofnumber)), family= "poisson",  data= data)
IndividualActivity2<-glm.nb(IndividualActivity ~ group +offset(log(Countofnumber)), data= data)

summary(IndividualActivity1)
summary(IndividualActivity2)

#likelyhood ratio test and Dean's test
test.nb.pois(IndividualActivity2, IndividualActivity1)
DeanB(IndividualActivity1, alternative="greater")
DeanB2(IndividualActivity1, alternative="greater")

#exponentialised coefficients
x<-coef(IndividualActivity2)
exp(x)

tab_model(IndividualActivity1,show.ci = 0.95,string.ci=TRUE, show.se = TRUE)
tab_model(IndividualActivity2,show.ci = 0.95,string.ci=TRUE, show.se = TRUE)




#PartnerActivity Poisson / Negative Binomial
mean(data1$PartnerActivity)
var(data1$PartnerActivity)
mean(data2$PartnerActivity)
var(data2$PartnerActivity)
mean(data3$PartnerActivity)
var(data3$PartnerActivity)

PartnerActivity1<-glm(PartnerActivity~ group+offset(log(Countofnumber)), family= "poisson",  data= data)
PartnerActivity2<-glm.nb(PartnerActivity ~ group +offset(log(Countofnumber)), data= data)
PartnerActivity3<-zeroinfl(PartnerActivity~ group+offset(log(Countofnumber)), dist= "negbin",  data= data)
PartnerActivity4<-hurdle(PartnerActivity~ group+offset(log(Countofnumber)), dist= "negbin",  data= data)

summary(PartnerActivity1)
summary(PartnerActivity3)


#likelyhood ratio test and Dean's test
test.nb.pois(PartnerActivity2, PartnerActivity1)
DeanB(PartnerActivity1, alternative="greater")
DeanB2(PartnerActivity1, alternative="greater")

#exponentialised coefficients
x<-coef(PartnerActivity2)
exp(x)

tab_model(PartnerActivity1,show.ci = 0.95,string.ci=TRUE, show.se = TRUE)
tab_model(PartnerActivity2,show.ci = 0.95,string.ci=TRUE, show.se = TRUE)
tab_model(PartnerActivity3,show.ci = 0.95,string.ci=TRUE, show.se = TRUE)
tab_model(PartnerActivity4,show.ci = 0.95,string.ci=TRUE, show.se = TRUE)

InterceptModelPartnerActivity4 <- update(PartnerActivity4, . ~ 1) 
logLik(PartnerActivity4)
logLik(InterceptModelPartnerActivity4)

1-logLik(PartnerActivity4)/(logLik(InterceptModelPartnerActivity4))

pchisq(2 * (logLik(PartnerActivity4) - logLik(PartnerActivity1)), df = 1, lower.tail = FALSE)



#SmallSidedActivity Poisson / Negative Binomial
mean(data1$SmallSidedActivity)
var(data1$SmallSidedActivity)
mean(data2$SmallSidedActivity)
var(data2$SmallSidedActivity)
mean(data3$SmallSidedActivity)
var(data3$SmallSidedActivity)

SmallSidedActivity1<-glm(SmallSidedActivity~ group+offset(log(Countofnumber)), family= "poisson",  data= data)
SmallSidedActivity2<-glm.nb(SmallSidedActivity ~ group +offset(log(Countofnumber)), data= data)
SmallSidedActivity3<-zeroinfl(SmallSidedActivity ~ group+offset(log(Countofnumber)), dist= "negbin",  data= data)
SmallSidedActivity4<-hurdle(SmallSidedActivity ~ group+offset(log(Countofnumber)), dist= "negbin",  data= data)


summary(SmallSidedActivity1)
summary(SmallSidedActivity3)

#likelyhood ratio test and Dean's test
test.nb.pois(SmallSidedActivity2, SmallSidedActivity1)
DeanB(SmallSidedActivity1, alternative="greater")
DeanB2(SmallSidedActivity1, alternative="greater")

#exponentialised coefficients
x<-coef(SmallSidedActivity2)
exp(x)

tab_model(SmallSidedActivity1,show.ci = 0.95,string.ci=TRUE, show.se = TRUE)
tab_model(SmallSidedActivity2,show.ci = 0.95,string.ci=TRUE, show.se = TRUE)
tab_model(SmallSidedActivity3,show.ci = 0.95,string.ci=TRUE, show.se = TRUE)
tab_model(SmallSidedActivity4,show.ci = 0.95,string.ci=TRUE, show.se = TRUE)


InterceptModelSmallSidedActivity4 <- update(SmallSidedActivity4, . ~ 1) 
logLik(SmallSidedActivity4)
logLik(InterceptModelSmallSidedActivity4)

1-logLik(SmallSidedActivity4)/(logLik(InterceptModelSmallSidedActivity4))


pchisq(2 * (logLik(SmallSidedActivity4) - logLik(PartnerActivity1)), df = 1, lower.tail = FALSE)





#WholeClassActivity Poisson / Negative Binomial
mean(data1$WholeClassActivity)
var(data1$WholeClassActivity)
mean(data2$WholeClassActivity)
var(data2$WholeClassActivity)
mean(data3$WholeClassActivity)
var(data3$WholeClassActivity)

WholeClassActivity1<-glm(WholeClassActivity~ group+offset(log(Countofnumber)), family= "poisson",  data= data)
WholeClassActivity2<-glm.nb(WholeClassActivity ~ group +offset(log(Countofnumber)), data= data)

summary(WholeClassActivity1)
summary(WholeClassActivity2)

#likelyhood ratio test and Dean's test
test.nb.pois(WholeClassActivity2, WholeClassActivity1)
DeanB(WholeClassActivity1, alternative="greater")
DeanB2(WholeClassActivity1, alternative="greater")

#exponentialised coefficients
x<-coef(WholeClassActivity2)
exp(x)

tab_model(WholeClassActivity1,show.ci = 0.95,string.ci=TRUE, show.se = TRUE)
tab_model(WholeClassActivity2,show.ci = 0.95,string.ci=TRUE, show.se = TRUE)








#WaitingActivity Poisson / Negative Binomial
mean(data1$WaitingActivity)
var(data1$WaitingActivity)
mean(data2$WaitingActivity)
var(data2$WaitingActivity)
mean(data3$WaitingActivity)
var(data3$WaitingActivity)

WaitingActivity1<-glm(WaitingActivity~ group+offset(log(Countofnumber)), family= "poisson",  data= data)
WaitingActivity2<-glm.nb(WaitingActivity ~ group +offset(log(Countofnumber)), data= data)

summary(WaitingActivity1)
summary(WaitingActivity2)


#likelyhood ratio test and Dean's test
test.nb.pois(WaitingActivity2, WaitingActivity1)
DeanB(WaitingActivity1, alternative="greater")
DeanB2(WaitingActivity1, alternative="greater")

#exponentialised coefficients
x<-coef(WaitingActivity2)
exp(x)

tab_model(WaitingActivity1,show.ci = 0.95,string.ci=TRUE, show.se = TRUE)
tab_model(WaitingActivity2,show.ci = 0.95,string.ci=TRUE, show.se = TRUE)



#ChildrenOffTask Poisson / Negative Binomial
mean(data1$ChildrenOffTask)
var(data1$ChildrenOffTask)
mean(data2$ChildrenOffTask)
var(data2$ChildrenOffTask)
mean(data3$ChildrenOffTask)
var(data3$ChildrenOffTask)

ChildrenOffTask1<-glm(ChildrenOffTask~ group+offset(log(Countofnumber)), family= "poisson",  data= data)
ChildrenOffTask2<-glm.nb(ChildrenOffTask ~ group +offset(log(Countofnumber)), data= data)

summary(ChildrenOffTask1)
summary(ChildrenOffTask2)

#likelyhood ratio test and Dean's test
test.nb.pois(ChildrenOffTask2, ChildrenOffTask1)
DeanB(ChildrenOffTask1, alternative="greater")
DeanB2(ChildrenOffTask1, alternative="greater")

#exponentialised coefficients
x<-coef(ChildrenOffTask2)
exp(x)

tab_model(ChildrenOffTask1,show.ci = 0.95,string.ci=TRUE, show.se = TRUE)
tab_model(ChildrenOffTask2,show.ci = 0.95,string.ci=TRUE, show.se = TRUE)







```

## Results


```{r pressure, echo=FALSE, warning=FALSE}


summary(IndividualActivity2)
tab_model(IndividualActivity2,show.ci = 0.95,string.ci=TRUE, show.se = TRUE)
nagelkerke(IndividualActivity2)


summary(PartnerActivity4)
tab_model(PartnerActivity4,show.ci = 0.95,string.ci=TRUE, show.se = TRUE)
1-logLik(PartnerActivity4)/(logLik(InterceptModelPartnerActivity4))


summary(SmallSidedActivity4)
tab_model(SmallSidedActivity4,show.ci = 0.95,string.ci=TRUE, show.se = TRUE)
1-logLik(SmallSidedActivity4)/(logLik(InterceptModelSmallSidedActivity4))


summary(WholeClassActivity2)
tab_model(WholeClassActivity2,show.ci = 0.95,string.ci=TRUE, show.se = TRUE)
nagelkerke(WholeClassActivity2)


summary(WaitingActivity2)
tab_model(WaitingActivity2,show.ci = 0.95,string.ci=TRUE, show.se = TRUE)
nagelkerke(WaitingActivity2)


summary(ChildrenOffTask2)
tab_model(ChildrenOffTask2,show.ci = 0.95,string.ci=TRUE, show.se = TRUE)
nagelkerke(ChildrenOffTask2)







```


