################################################################################
# Updated version of the R code for the analysis in:
#
#   "Attributable risk from distributed lag models"
#   Antonio Gasparrini & Michela Leone
#   BMC Medical Research Methodology 2014
#   http://www.ag-myresearch.com/2014_gasparrini_bmcmrm.html
#
# Update: 15 January 2017
# * an updated version of this code, compatible with future versions of the
#   software, is available at:
#   https://github.com/gasparrini/2014_gasparrini_BMCmrm_Rcodedata
#
#   
#
################################################################################

#modified by Kaisa Lakkala and Reija Ruuhela
#Latest modification 14 January 2026

library(dlnm) ; library(splines) ; library(foreign) ; library(tsModel)
library(ggplot2); library(dplyr); library(lubridate) ; library(mgcv); library(splines);library(nlme) 

# CHECK VERSION OF THE PACKAGE
if(packageVersion("dlnm")<"2.2.0")
  stop("update dlnm package to version >= 2.2.0")

# Load the suicide deaths and radiation data

lndn2 <- read.csv2("../itsemurhadata/suicidedata.csv",sep=",", head=TRUE, na.strings = -999) #suicides
absdata <- read.csv2("../itsemurhadata/radiation_data.csv",sep=",", head=TRUE, na.strings = -999)#dailydose

lndn2$Dailydose <- as.numeric(absdata$Dailydose) 
lndn2$date <- make_date(lndn2$Year, lndn2$month, lndn2$day) 


#choose the time period
jakso2 <- subset(lndn2, date>=as.Date("1979-01-01") & date<=as.Date("2019-12-31"))
lndn3 <- jakso2
 
#difference between two consecutive days
uusi <- (lndn3$Dailydose - dplyr::lag(lndn3$Dailydose, n = 1))
lndn3$Dailydose <- uusi

#choose the months March-June
lndn <- lndn3[lndn3$month%in%3:6, ]

 
###################################
# DERIVE THE CROSS-BASIS


# KNOTS FOR EXPOSURE-RESPONSE FUNCTION
vk <- equalknots(lndn$Dailydose,nk=1) #ns on default

# KNOTS FOR THE LAG-RESPONSE FUNCTION
maxlag <- 7 
lk <- logknots(maxlag,nk=1) 

# CENTERING VALUE (AND PERCENTILE)
cen <- 0

cb <- crossbasis(lndn$Dailydose, lag=maxlag, argvar=list(fun="ns",knots=vk), arglag=list(fun="ns",knots=lk)) 

# SUMMARY
summary(cb)


####################################
# RUN THE MODEL AND OBTAIN PREDICTIONS

model <- glm(dead_all~cb+ns(number,4*41)+dow,family=quasipoisson(),lndn, na.action="na.exclude")

####################################
# https://github.com/gasparrini/2010_gasparrini_StatMed_Rcode
# following Gasparini et al. 2010 update for 

cp <- crosspred(cb,model,at=seq(min(lndn$Dailydose,na.rm=T),max(lndn$Dailydose,na.rm=T),length=50),cen=cen)

##############################
# RESULTS AND PLOTS
##############################


#Contour-plot
cairo_pdf("./kuvat/UV_contour_1979_2019_ALL.pdf")
cp <- crosspred(cb,model,at=seq(min(lndn$Dailydose,na.rm=T),max(lndn$Dailydose,na.rm=T),length=50),cen=cen)
plot(cp, "contour", xlab="Difference in UV dose UVR doses[J/m²]", ylab="Lag [day]", key.title=title("RR"), cex.lab=1.5, cex.axis=1.5)
dev.off()

# SLICES
percentiles <- round(quantile(lndn$Dailydose,c(0.01,0.05,0.95,0.99),na.rm = TRUE),1)
cp <- crosspred(cb,model,at=(min(lndn$Dailydose,na.rm=T)*10):(max(lndn$Dailydose,na.rm=T)*10)/10,cen=cen)
cairo_pdf("./kuvat/UV_slices_1979_2019_ALL.pdf", 
           width = 6,        # Width in inches
           height = 8,       # Height in inches
           pointsize = 12) 
plot(cp,var=percentiles,lag=c(0,1,2,3), ylab="RR",ylim=c(0.8,1.5),cex.lab=1.5, cex.axis=1.5)
mtext("Var=Difference in UVR doses [J/m²]", side = 1, line = 4.3)
dev.off()


# OVERALL EFFECT
percentiles <- round(quantile(lndn$Dailydose,c(0.01,0.05,0.95,0.99),na.rm = TRUE),1)
cp <- crosspred(cb,model,at=seq(min(lndn$Dailydose,na.rm=T),max(lndn$Dailydose,na.rm=T),length=50),cen=cen)

cairo_pdf("./kuvat/UV_overall_1979_2019_ALL.pdf")
plot(cp,"overall",xlab="Difference in UV dose UVR doses[J/m²]",ylim=c(-2,4),ylab="RR",
     main="Overall effect, 1979-2019")
abline(v=quantile(lndn$Dailydose,c(0.01,0.99),na.rm = TRUE ),lty=2)
abline(v=quantile(lndn$Dailydose,c(0.05,0.95),na.rm = TRUE ),lty=4)
abline(v=min(lndn$Dailydose,na.rm=T))
abline(v=max(lndn$Dailydose,na.rm=T))
abline(v=0)
par(new=T)

hist(lndn$Dailydose,xlim=c(min(lndn$Dailydose,na.rm=T),max(lndn$Dailydose,na.rm=T)),ylim=c(0,1500),axes=F,ann=F,col=grey(0.95),breaks=100)
abline(v=quantile(lndn$Dailydose,c(0.01,0.99),na.rm = TRUE ),lty=2)
abline(v=quantile(lndn$Dailydose,c(0.05,0.95),na.rm = TRUE ),lty=4)
abline(v=min(lndn$Dailydose,na.rm=T))
abline(v=max(lndn$Dailydose,na.rm=T))
abline(v=0)
dev.off()

# RR AT CHOSEN PERCENTILES VERSUS 0 radiation difference (AND 95%CI)
percentiles <- round(quantile(lndn$Dailydose,c(0.01,0.025,0.05,0.95,0.975,0.99),na.rm = TRUE),1)
c(min(lndn$Dailydose,na.rm=T),max(lndn$Dailydose,na.rm=T))
cp <- crosspred(cb,model,at=-40000:40000/10,cen=cen)
lag <- 0
cp$matRRfit[as.character(percentiles), lag + 1] 
cp$matRRlow[as.character(percentiles),lag + 1]
cp$matRRhigh[as.character(percentiles),lag + 1]
M=max(lndn$Dailydose,na.rm=T)
cp$matRRfit[as.character(M),lag + 1] 
cp$matRRlow[as.character(M),lag + 1]
cp$matRRhigh[as.character(M),lag + 1]


#THIS IS TO CALCULATE THE OVERALL INCREASE:
cp$allRRfit[as.character(percentiles)]
cbind(cp$allRRlow,cp$allRRhigh)[as.character(percentiles),]
#at max diff
cp$allRRfit[as.character(max(lndn$Dailydose,na.rm=T))]
cbind(cp$allRRlow,cp$allRRhigh)[as.character(max(lndn$Dailydose,na.rm=T)),]



####################################################################
# FUNCTION TO COMPUTE THE Q-AIC IN QUASI-POISSON MODELS
fqaic <- function(model) {
  loglik <- sum(dpois(model$y,model$fitted.values,log=TRUE))
  phi <- summary(model)$dispersion
  qaic <- -2*loglik + 2*summary(model)$df[3]*phi
  return(qaic)
}
# KNOTS GRID 
grid <- as.matrix(expand.grid(var=1:6,lag=1:6))
# SEARCH (TAKES )
system.time({
  aicval <- sapply(seq(nrow(grid)), function(i) {
    vk <- equalknots(lndn3$Dailydose,nk=grid[i,1])

    lk <- logknots(7,nk=grid[i,2])

    cb <- crossbasis(lndn3$Dailydose, lag=7, argvar=list(fun="ns",knots=vk), arglag=list(fun="ns",knots=lk)) 
    m <- glm(dead_all~cb+ns(number,4*41)+dow,family=quasipoisson(),lndn3, na.action="na.exclude") #1979-1990->12,1979-2019->41, 1991-2019->29
    return(fqaic(m))
  })
})
# BEST FITTING MODEL
(best <- grid[which.min(aicval),])
plot(aicval,col=3,pch=19)
best