## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup, include=FALSE-----------------------------------------------------
library(LongitudinalEvalue)
library(yuima)
library(ggplot2)
library(np)
library(matrixStats)

data("Covidexample", package="LongitudinalEvalue")

## ----TATE, fig.cap="The temporal average treatment effect (TATE) over time between the exposure regional lockdown and outcome incidence of SARS Covid-19", echo=FALSE----
ggplot(data = Covidexample) +
  geom_line(aes(x = time, y = TATE)) +
    ylab("Absolute difference in incidence pr. 100.000") +
  xlab("Date") +
  ggtitle("Temporal average treatment effect")+
  theme_minimal()

## ----difplot, out.width="33%", include=TRUE, fig.align="center", fig.cap=c("Difference in difference plots, stating the differencing order such that the TATE process becomes stationary. From the plots stationarity is achieved when the difference in difference operator is used twice, as the spline (blue line) reflects no temporal dependence, when the order is 2"),fig.show="hold",warning=FALSE,echo=FALSE----
knitr::include_graphics(c("DiD_order_0.pdf","DiD_order_1.pdf","DiD_order_2.pdf"))

## ----senscurves, out.width="80%", include=TRUE, fig.align="center", fig.cap=c("The confounder plot, stating the hypothetical association population density needs to have with the regional lockdown and incidence to explain the TATE between the regional lock down and incidence."), echo=FALSE,warning=FALSE----
knitr::include_graphics("Confounder_order_2.pdf")

## ----qmlefitting, eval=TRUE---------------------------------------------------
Delta<-1/3
Covidexample$t<-1:30
mod<-setModel(drift="a/TATE+c*TATE", diffusion="sigma", state.var="TATE", time.var = "t", 
              solve.var = "TATE",xinit = 5.425966)
model<-setYuima(model=mod, data=setData(zoo(Covidexample$TATE, order.by = Covidexample$t),
                                        delta=Delta))
fit <- qmle(model,start = list(a=-5,c = 5,sigma = 0.2),
            lower = list(a=-10,c = -10,sigma = 0), 
            upper = list(a=0,c = 10,sigma =1))


## ----simulated values of the confounder, eval=TRUE, warning=FALSE-------------
U_t<-SDEconfounderYuima1d(fit, EFXU=1.2, EFUY=1.5,time=Covidexample$t,
                          initial.value = 5.425966/(1.2*1.5))

## ----paraplot, eval=TRUE, warning=FALSE, fig.cap="The trajectory based on SDEconfounderYuima1d of an unmeasured confounder capable of explaining the association between regional lockdown and incidence. The association between the exposure and confounder was set to 1.2 on the risk ratio scale and the association between the outcome and confounder was set to 1.5 on the absolute scale."----
plot(U_t)

## ----plot_unpara visible, eval=FALSE------------------------------------------
# U_t<-SDEconfounderNonparametric(Covidexample$TATE, ATE=Covidexample$ATE,Covidexample$t,
#                                 EFXU=1.2, EFUY=1.5, scheme = "Euler-Maruyama")
# plot(U_t)

## ----plotunpara creation, eval=TRUE, include=FALSE----------------------------
U_t<-SDEconfounderNonparametric(Observed=Covidexample$TATE,
ATE=Covidexample$ATE,time=Covidexample$t, EFXU=1.2,EFUY=1.5, scheme = "Euler-Maruyama")

## ----plotunpara, eval=TRUE, echo=FALSE,fig.cap="The trajectory based on SDEconfounderNonparametric of an unmeasured confounder capable of explaining the association between regional lockdown and incidence. The association between the exposure and confounder was set to 1.2 on the risk ratio scale and the association between the outcome and confounder was set to 1.5 on the absolute scale."----
plot(U_t)

## ----unparatimevarying, eval=FALSE--------------------------------------------
# U_t<-SDEconfounderNonparametric(Covidexample$TATE, ATE=Covidexample$ATE,Covidexample$time,
#                                 EFXU=1, EFUY=1, scheme = "Euler-Maruyama")

## ----unparatimevarying hidden, eval=TRUE, include=FALSE-----------------------
U_t<-SDEconfounderNonparametric(Covidexample$TATE, ATE=Covidexample$ATE,Covidexample$time, 
                                EFXU=1, EFUY=1, scheme = "Euler-Maruyama")

## ----timevaryingplot, fig.cap="The trajectory based of an unmeasured confounder capable of explaining the association between regional lockdown and incidence, where the association from the confounder to the exposure and outcome is time-varying."----

beta<-1.2
alpha<-0.3
gamma<-0.5

U_mat <- matrix(NA_real_, nrow = 30, ncol = 100)


for (k in 2:100){


 E <- U_t@trajectories[, k]


U <- numeric(30)

U[1] <- E[1] * 0.27

for (t in 1:29) {
  denom <- (beta * U[t]^(beta - 1)) *
    ((1/t) * (U[t]^beta + alpha*t) *
       exp(U[t]^beta + (alpha - gamma)*t))
  
  U[t + 1] <- (E[t + 1] - E[t]) / denom + U[t]
}
  U_mat[, k] <- U


}

alpha_level<-0.05

U_mat<-t(as.matrix(U_mat))
U_mat<-U_mat[2:100,]
Mean <- colMeans2(U_mat, na.rm = TRUE)
CIll <- colQuantiles(U_mat, probs = alpha_level/2, na.rm = TRUE)
CIul <- colQuantiles(U_mat, probs = 1 - alpha_level/2, na.rm = TRUE)



  Confounderdataset<-data.frame(
    Confounder =  Mean,
    time = Covidexample$t,
    ConfounderCIll =  CIll,
    ConfounderCIul = CIul
  )



names(Confounderdataset)<-c("mean_Confounder","time", "CIll", "CIul")

ggplot(Confounderdataset, aes(x = time, y = mean_Confounder)) +
  geom_line(color = "blue", linewidth = 1) +
  geom_ribbon(
    aes(ymin = CIll, ymax = CIul),
    fill = "lightblue",
    alpha = 0.5
  ) +
  labs(
    x = "Time",
    y = "Confounder Value",
    title = "Confounder with Credibility Intervals"
  ) +
  theme_minimal()