rm(list=ls())
setwd("D:/Research/PR 11")
data = read.csv("jakselchjam.csv", sep=";")
#View(datajakpus)

library(dlnm)
library(splines)
summary(data$rain)
crosshujan = crossbasis(data$rain, lag=12, argvar = list(fun="poly"), arglag = list(fun="ns"))

#varknots <- equalknots(data$suhu,fun="bs",df=6,degree=2)
#lagknots <- logknots(6, 1)
#crosssuhu = crossbasis(data$suhu, lag=6, 
                       #argvar=list(fun="bs",knots=varknots),
                       #arglag=list(knots=lagknots))
summary(crosshujan)
#crossangin = crossbasis(datajakpus$wind, lag =14, argvar = list(fun="lin"), arglag = list(fun="strata"))
#summary(crossangin)

#rumus = glm(PM ~ crosshujan + datajakpus$wind +ns(month, 2), family = quasipoisson(), datajakpus)
data$rain = as.numeric(data$rain)
nzmean <- sum(data$rain)/sum(!!data$rain)
#datajakpus$wind = as.numeric(datajakpus$wind)
data$PM = as.numeric(data$PM)
data$suhu = as.numeric(data$suhu)

rumus = glm(PM ~ crosshujan + data$suhu + ns(time,2), family = quasipoisson(), data)

predhujan = crosspred(crosshujan, rumus, cen=nzmean, by=1, at=0:20, ci.level=0.95, cumul=T)
plot(predhujan, xlab="Hourly Precipitation (mm)", zlab="Relative Effect",
     theta=380, phi=20, 
     lphi=500, main="3D graph of rain effect")
ouput<-predhujan[["matRRfit"]]