Plottar besök per årskull + Poissonmodeller

install.packages("pscl")
library(dplyr)
library(MASS)
library(pscl)

Besokarskull <- Beso_k_per_a_rskull_Anpassad_fo_r_R_och_SAS
plot(Besokarskull$Ålder,Besokarskull$tot_total_tot)
BesokUV <- filter(Besokarskull, Kommun=="Upplands-Väsby")
BesokOST <- filter(Besokarskull, Kommun=="Östersund")
BesokTOM <- filter(Besokarskull, Kommun=="Tomelilla")
BesokMLM <- filter(Besokarskull, Kommun=="Malmö")

plot(BesokMLM$Ålder, BesokMLM$tot_total_tot, type ="n")
lines(BesokMLM$Ålder, BesokMLM$tot_total_tot, col="black")
lines(BesokUV$Ålder, BesokUV$tot_total_tot, col="blue")
lines(BesokOST$Ålder, BesokOST$tot_total_tot, col="green")
lines(BesokTOM$Ålder, BesokTOM$tot_total_tot, col="red")

plot(BesokUV$Ålder, BesokUV$tot_total_tot, type ="n")
lines(BesokUV$Ålder, BesokUV$tot_total_tot, col="blue")
lines(BesokOST$Ålder, BesokOST$tot_total_tot, col="green")
lines(BesokTOM$Ålder, BesokTOM$tot_total_tot, col="red")

plot(BesokTOM$Ålder, BesokTOM$tot_total_tot, type ="n")
lines(BesokTOM$Ålder, BesokTOM$tot_total_tot, col="red")

#Import and rename dataframe.
Poissonbesok <- Sammansta_llning_2
#Creating expected count to use as offset
Poissonbesok$E <- Poissonbesok$`Befolkning 2019` * (sum(Poissonbesok$`total 2019`)/sum(Poissonbesok$`Befolkning 2019`))
Poissonbesok$ptatortsgrad <- Poissonbesok$`Tätortsgrad, 2019`/100
#Poisson model
fit.poi.besok<-glm(Poissonbesok$`total 2019`~offset(log(Poissonbesok$E))+scale(Poissonbesok$`Andel inv 65+, 2019`)+
scale(Poissonbesok$`Andel inv 0-18 år, 2019`)+scale(Poissonbesok$ptatortsgrad),
family=poisson, data=Poissonbesok)
summary(fit.poi.besok)
#Negative Binomial.
fit.nb.besok<-glm.nb(Poissonbesok$`total 2019`~offset(log(Poissonbesok$E))+scale(Poissonbesok$`Andel inv 65+, 2019`)+
                     scale(Poissonbesok$`Andel inv 0-18 år, 2019`)+scale(Poissonbesok$ptatortsgrad), data=Poissonbesok)
summary(fit.nb.besok)

#Testing model fit.
pchisq(fit.nb.besok$deviance, df = fit.nb.besok$df.residual)

hist(fitted(fit.nb.besok)/Poissonbesok$E,main="Sweden SMR Distribution from Negative Binomial Model" )
hist(Poissonbesok$`total 2019`/Poissonbesok$E)

hist(fitted(fit.poi.besok)/Poissonbesok$E,main="Sweden SMR Distribution from Poisson Model" )
#QuasiPoisson.
fit.quasi.besok<-glm(Poissonbesok$`total 2019`~offset(log(Poissonbesok$E))+scale(Poissonbesok$`Andel inv 65+, 2019`)+
                       scale(Poissonbesok$`Andel inv 0-18 år, 2019`)+scale(Poissonbesok$ptatortsgrad),
                     family=quasipoisson, data=Poissonbesok)
summary(fit.quasi.besok)
#Negative binomial with befolkning(n) as off set.
fit.nb.besok2<-glm.nb(Poissonbesok$`total 2019`~offset(log(Poissonbesok$`Befolkning 2019`))+scale(Poissonbesok$`Andel inv 65+, 2019`)+
                       scale(Poissonbesok$`Andel inv 0-18 år, 2019`)+ scale(Poissonbesok$`Mediannettoinkomst, kr/inv 20+, 2018`)+ factor(Poissonbesok$`Kommunindelning i tre typer, Namn`), data=Poissonbesok)

summary(fit.nb.besok2)

#Test fit Neg Binomial vs.Poisson, this is also a test of equidispersion.
odTest(fit.nb.besok)