library(readxl)
rmark_ktm <- read_excel("rmark_ktm.xls") #Loading excel datasheet
View(rmark_ktm)

#New Stats
library(Hmisc)
describe(rmark_ktm) 
cor(rmark_ktm, method = c("spearman"),use = "complete.obs")

#Linear regression
Trendcheck <- lm(formula = dist ~ speed, data = cars)


#Kathmandu Post

modelKPost1 <- glm(kathmandu_post ~Red_Panda+ Tiger+ Musk_Deer+ Rhino+ Bear+ Others+ Owl+ Leopard+ Pangolin ,family=binomial(link='logit'),data= rmark_ktm)
summary(modelKPost1)

modelKPost2 <- glm(kathmandu_post ~arrestsize+ plevel+ year,family=binomial(link='logit'),data= rmark_ktm)
summary(modelKPost2)

UnimodelKpost3 <- glm(kathmandu_post ~arrestees+ plevel+ year+ Red_Panda+ Tiger+ Musk_Deer+ Rhino+ Bear+ Others+ Owl+ Leopard+ Pangolin ,family=binomial(link='logit'),data= rmark_ktm)
summary(UnimodelKpost3)

#Kantipur

modelKant1 <- glm(kantipur ~Red_Panda+ Tiger+ Musk_Deer+ Rhino+ Bear+ Others+ Owl+ Leopard+ Elephant,family=binomial(link='logit'),data= rmark_ktm)
summary(modelKant1)

modelKant2 <- glm(kantipur ~arrestsize+ plevel+ year,family=binomial(link='logit'),data= rmark_ktm)
summary(modelKant2)

UnimodelKant2 <- glm(kantipur ~arrestees+ plevel+ year+ Red_Panda+ Tiger+ Musk_Deer+ Rhino+ Bear+ Others+ Owl+ Leopard+ Elephant,family=binomial(link='logit'),data= rmark_ktm)
summary(UnimodelKant2)



#Gorkahapatra

modelgorkha1 <- glm(gorakhapatra ~Red_Panda+ Tiger+ Musk_Deer+ Others+ Owl+  Pangolin ,family=binomial(link='logit'),data= rmark_ktm)
summary(modelgorkha1)

modelgorkha2 <- glm(gorakhapatra ~arrestsize+ plevel+ year,family=binomial(link='logit'),data= rmark_ktm)
summary(modelgorkha2)

Unimodelgorkha <- glm(gorakhapatra ~arrestees+ plevel+ year+ Red_Panda+ Tiger+ Musk_Deer+ Others+ Owl+  Pangolin ,family=binomial(link='logit'),data= rmark_ktm)
summary(Unimodelgorkha)

#General

library (glmm)

modelglobal1 <- glm(AllSources ~Red_Panda+ Tiger+ Musk_Deer+ Rhino+ Bear+ Others+ Owl+ Leopard+ Pangolin+ Elephant ,family=binomial(link='logit'),data= rmark_ktm)
summary(modelglobal1)


modelglobal <- glm(AllSources ~arrestsize+ plevel+ year,family=binomial(link='logit'),data= rmark_ktm)
summary(modelglobal)

Unimodelglobal <- glm(AllSources ~arrestees + plevel+ year+ Red_Panda+ Tiger+ Musk_Deer+ Rhino+ Bear+ Others+ Owl+ Leopard+ Pangolin+ Elephant ,family=binomial(link='logit'),data= rmark_ktm)
summary(Unimodelglobal)


#total Seizures by species

TotalSpecies <- c(13, 14, 18, 15, 43, 24, 71, 23, 24, 50)
Totalchi <- chisq.test(TotalSpecies, p = c(1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10))
Totalchi

TotalSpeciesres <- Totalchi$stdres
TotalSpeciesres

TotalSpeciessig<-.05
TotalSpeciessigAdjust<-TotalSpeciessig/10
TotalSpeciessigAdjust
qnorm(TotalSpeciessigAdjust/2)

#Reported Seizures by species

RepSpecies <- c(2, 8, 6, 3, 10, 3, 16, 2, 5, 8)
RepChi <- chisq.test(RepSpecies, p = c(1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10))
RepChi

RepChires <- RepChi$stdres
RepChires

RepChisig<-.05
RepChisigAdjust<-RepChisig/10
RepChisigAdjust
qnorm(RepChisigAdjust/2)


#UnReported Seizures by species

UNRepSpecies <- c(11, 6, 12, 12, 33, 21, 55, 21, 19, 42)
UNRep <- chisq.test(UNRepSpecies , p = c(1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10, 1/10))
UNRep

UNRepSpecies <- c(2, 8, 6, 3, 10, 3, 16, 2, 5, 8)
UNRepres <- UNRep$stdres
UNRepres

UNRepsig4<-.05
UNRepsigAdjust4<-UNRepsig4/10
UNRepsigAdjust4
qnorm(UNRepsigAdjust4/2)


#Total Seizures by newspaper

Species <- c(24, 15, 24)
Test3 <- chisq.test(Species, p = c(1/3, 1/3, 1/3))
Test3

Species <- c(24, 15, 24)
Test3res <- chisq.test(Species, p = c(1/3, 1/3, 1/3))$stdres
Test3res

sig3<-.05
sigAdjust3<-sig3/50
sigAdjust3
qnorm(sigAdjust3/2)



#piecewise regression

library(dpseg)
segs <- dpseg(x=Seizure_Data_monthwise$Time, y=Seizure_Data_monthwise$Total, minl=6, type = "cor")
print(segs)
plot(segs)

LinearModel <- lm(Seizure_Data_monthwise$Total ~ Seizure_Data_monthwise$Time, data = Seizure_Data_monthwise)
summary(LinearModel)


#correlation between variables
cor(rmark_ktm, method="pearson")


#correlation between species
library(readxl)
species_rktm <- read_excel("species_rktm") #Loading excel datasheet
View(rmark_ktm)

cor(species_rktm, method="pearson")


#Reported vs Actual seizures comparison


#three Newspapers

Species <- c(2, 8, 6, 3, 10, 3, 16, 2, 5, 8)
ThreeNews <- chisq.test(Species, p = c(0.04, 0.05, 0.06, 0.05, 0.15, 0.08, 0.24, 0.08, 0.08, 0.17), simulate.p.value=T, B = 10000)
ThreeNews

Species <- c(2, 8, 6, 3, 10, 3, 16, 2, 5, 8)
ThreeNewsres <- chisq.test(Species, p = c(0.04, 0.05, 0.06, 0.05, 0.15, 0.08, 0.24, 0.08, 0.08, 0.17))$stdres
ThreeNewsres

sig3<-.05
sigAdjust3<-sig3/10
sigAdjust3
qnorm(sigAdjust3/2)

##Observed <- c(2, 8, 6, 3, 10, 3, 16, 2, 5, 8)
##Real <- c(0.04, 0.05, 0.06, 0.05, 0.15, 0.08, 0.24, 0.08, 0.08, 0.17)
##Test3 <- chisq.test(x= Observed, p= Real, simulate.p.value=T, B = 10000)
##Test3

#single Newspapers

#KTMP

ObsKTIMES <- c(0, 4, 4, 1, 3, 1, 9, 1, 1, 1)
KTIMES <- chisq.test(ObsKTIMES, p = c(4/100, 5/100, 6/100, 5/100, 15/100, 8/100, 24/100, 8/100, 8/100, 17/100), simulate.p.value=T, B = 10000)
KTIMES

##ObservedKTMP1 <- c(0, 4, 4, 1, 3, 1, 9, 1, 1, 1)
##ExpectedKTMP1 <- c(4/100, 5/100, 6/100, 5/100, 15/100, 8/100, 24/100, 8/100, 8/100, 17/100)
##ChiKTMP1 <- chisq.test(x= ObservedKTMP1, p= ExpectedKTMP1, simulate.p.value=T, B = 10000)

KTIMESres <- chisq.test(x= ObservedKTMP1 , p= ExpectedKTMP1, simulate.p.value=T, B = 200000)$stdres
KTIMES

sigKTIMES<-.05
sigAdjustKTIMES<-sigKTIMES/10
sigAdjustKTIMES
qnorm(sigAdjustKTIMES/2)

#GRPT

ObsGRPT <- c(0,1,0,1,5,2,0,0,2,4)
GRPT <- chisq.test(ObsGRPT, p = c(0.04, 0.05, 0.06, 0.05, 0.15, 0.08, 0.24, 0.08, 0.08, 0.17), simulate.p.value=T, B = 10000)
GRPT


##ObservedGRPT <- c(0,1,0,1,5,2,0,0,2,4)
##RealGRPT <- c(0.04, 0.05, 0.06, 0.05, 0.15, 0.08, 0.24, 0.08, 0.08, 0.17)
##ChiGRPT <- chisq.test(x= ObservedGRPT, p=RealGRPT, simulate.p.value=T, B= 200000)
##ChiGRPT

#KNTP

ObsKNTP <- c(2,3,2,1,4,0,7,1,2,3)
KNTP <- chisq.test(ObsKNTP, p = c(0.04, 0.05, 0.06, 0.05, 0.15, 0.08, 0.24, 0.08, 0.08, 0.17), simulate.p.value=T, B = 10000)
KNTP

##ObservedKNTP <- c(2,3,2,1,4,0,7,1,2,3)
##RealKNTP <- c(0.04, 0.05, 0.06, 0.05, 0.15, 0.08, 0.24, 0.08, 0.08, 0.17)
##ChiKNTP <- chisq.test(x= ObservedKNTP, p= RealKNTP, simulate.p.value=T, B= 200000)
##ChiKNTP

#KTP+GRPT
ObservedDuo1 <- c(0,5,4,2,8,3,9,1,3,5)
RealDuo1 <- c(0.04, 0.05, 0.06, 0.05, 0.15, 0.08, 0.24, 0.08, 0.08, 0.17)

ChiDuo1 <- chisq.test(x= ObservedDuo1, p= RealDuo1, simulate.p.value=T, B= 200000)
ChiDuo1


#KTP+KNTP  <----THIS ONE

KTPKNTPObs <- c(2,7,6,2,7,1,16,2,3,4)
KTPKNTP <- chisq.test(KTPKNTPObs, p = c(0.04, 0.05, 0.06, 0.05, 0.15, 0.08, 0.24, 0.08, 0.08, 0.17), simulate.p.value=T, B = 10000)
KTPKNTP

KTPKNTPres <- chisq.test(KTPKNTPObs, p = c(0.04, 0.05, 0.06, 0.05, 0.15, 0.08, 0.24, 0.08, 0.08, 0.17), simulate.p.value=T, B = 10000)$stdres
KTPKNTPres

sigKTPKNTPres<-.05
sigAdjustKTPKNTPres<-sigKTPKNTPres/10
sigAdjustKTPKNTPres
qnorm(sigAdjustKTPKNTPres/2)

##ObservedDuo2 <- c(2,7,6,2,7,1,16,2,3,4)
##RealDuo2 <- c(0.04, 0.05, 0.06, 0.05, 0.15, 0.08, 0.24, 0.08, 0.08, 0.17)
##ChiDuo2 <- chisq.test(x= ObservedDuo2, p= RealDuo2 , simulate.p.value=T, B = 200000)
##ChiDuo2

##ChiDuo2res2 <- chisq.test(x= ObservedDuo2, p = RealDuo2, simulate.p.value=T, B = 200000)$stdres
##ChiDuo2res2

##sigChiDuo2<-.05
##sigAdjustChiDuo2<-sigChiDuo2/10
##sigAdjustChiDuo2
##qnorm(sigAdjustChiDuo2/2)


#GRPT+KNTP
ObservedDuo3 <- c(2,4,2,2,9,2,7,1,4,7)
RealDuo3 <- c(0.04, 0.05, 0.06, 0.05, 0.15, 0.08, 0.24, 0.08, 0.08, 0.17)

ChiDuo3 <- chisq.test(x= ObservedDuo3, p= RealDuo3, simulate.p.value=T, B= 200000)
ChiDuo3


#test
ObservedKTMP <- c(0, 4, 4, 1, 3, 1, 9, 1, 1, 1)
ChiKTMP <- chisq.test(ObservedKTMP, p = c(4/100, 5/100, 6/100, 5/100, 15/100, 8/100, 24/100, 8/100, 8/100, 17/100), simulate.p.value=T, B = 200000)
ChiKTMP

ObservedKTMP <- c(0, 4, 4, 1, 3, 1, 9, 1, 1, 1)
ChiKTMP <- chisq.test(ObservedKTMP, p = c(4/100, 5/100, 6/100, 5/100, 15/100, 8/100, 24/100, 8/100, 8/100, 17/100), simulate.p.value=T, B = 200000)$stdres
ChiKTMP