GAM2year456factorweek

ggplot(y,aes(x=y$Start.Date,y=y$freq))+geom_line()

gampoisson<- mgcv ::gam(freq ~ s(day,bs = "cc", k = 350) +s(day,bs = "cc",by=weekdays)+
                      weekdays,
                    data = y,
                    family = "poisson")

dispersiontest(gampoisson)

## 
##  Overdispersion test
## 
## data:  gampoisson
## z = 15.03, p-value < 2.2e-16
## alternative hypothesis: true dispersion is greater than 1
## sample estimates:
## dispersion 
##   670.6098

gamnoar<- mgcv ::gam(freq ~ s(day,bs = "cc", k = 350) +s(day,bs = "cc", k = 5,by=weekdays)+
                      weekdays,
                    data = y,
                    family = "quasipoisson")

acf(gamnoar$residuals)

pacf(gamnoar$residuals)

#####################################################

gammar0<- mgcv ::gamm(freq ~ s(day,bs = "cc", k = 350) +s(day,bs = "cc", k = 5,by=weekdays)+
                       weekdays,
                     data = y,
                     family = "quasipoisson")

## 
##  Maximum number of PQL iterations:  20

## iteration 1

## iteration 2

## iteration 3

gamar1<- mgcv ::gamm(freq ~ s(day,bs = "cc", k = 350) +s(day,bs = "cc", k = 5,by=weekdays)+
                         weekdays,
                       data = y,
                       family = "quasipoisson",full=T,
                       correlation = corARMA(form = ~ 1|day, q = 1))

## 
##  Maximum number of PQL iterations:  20

## iteration 1

## iteration 2

## iteration 3

gamgu<-  mgcv ::gamm(freq ~ s(day,bs = "cc", k = 350) +s(day,bs = "cc", k = 5,by=weekdays)+
                       weekdays,
                     data = y,
                     family = "gaussian")


AIC(gamar1$lme,gammar0$lme,gamnoar,gampoisson)

##                    df         AIC
## gamar1$lme   17.00000   -339.2336
## gammar0$lme  16.00000   -340.6318
## gamnoar      58.01489          NA
## gampoisson  402.40771 771784.2964

qqnorm(gamgu$lme)

qqnorm(gamgu$lme)

pgo<- predict(gampoisson,newdata = data.frame(day=re2017$day, weekdays = re2017$weekdays))

pnoar<- exp(predict(gamnoar,newdata = data.frame(day=re2017$day, weekdays = re2017$weekdays)))

par0<- exp(predict(gammar0$gam,newdata = data.frame(day=re2017$day, weekdays = re2017$weekdays)))

par1<- exp(predict(gamar1$gam,newdata = data.frame(day=re2017$day, weekdays = re2017$weekdays)))

pgamgu<- predict(gamgu$gam,newdata = data.frame(day=re2017$day, weekdays = re2017$weekdays))

mean((re2017$freq-pgo)^2)

## [1] 921221504

mean((re2017$freq-pnoar)^2)

## [1] 34506195

mean((re2017$freq-par0)^2)

## [1] 33616281

mean((re2017$freq-par1)^2)

## [1] 33796118

mean((re2017$freq-pgamgu)^2)

## [1] 33182208

datas <- data.frame(re2017$freq,re2017$day,par0,pnoar,par1,pgamgu)
datas <- data.frame(re2017$freq,re2017$day,pgamgu)

ggplot(data=datas,aes(x=datas$re2017.day,y=datas$re2017.freq))+geom_line()+
  geom_line(data=datas,aes(x=datas$re2017.day,y=datas$pgamgu),col="red")+ggtitle("GAMM no p/q")