# Assess temporary effect of treatment

Imagine that I have a treatment that reduces the likelihood of response to a stimulus. This could be anything you like, but the simplest example is of a treatment (e.g., hand washing, mask wearing, etc.) that prevents disease when exposed.

For simplicity, I built an (overly complicated) model of the risk of responding to the stimulus over time under the two treatments (all of this is done in R using the tidyverse):

day_prob <-
list(
A = c(0.06, 0.06, pmax(dchisq(1:13, 3)*1, 0.06))
, B = c(0.06, 0.06, pmax(dchisq(seq(1, 25, 2), 5)*6, 0.06))
)


This can be plotted to show the risk:

day_prob %>%
lapply(function(x){
tibble(
Day = 1:length(x)
, Prob = x
)
}) %>%
bind_rows(.id = "Group") %>%
ggplot(aes(x = Day
, y = Prob
, col = Group)) +
geom_line() +
geom_point() +
scale_color_brewer(palette = "Dark2")


Plot:

Note that there is a similar baseline risk for both group (0.06), some latency to develop after exposure (risk rises on day 3), and that risk falls over time back to baseline.

Now, assuming that I randomize individuals into the two treatments, what is my best approach to identify this effect? I can just analyze each day separately, though that raises some repeated-measures questions since (unlike the sample data below) it is likely that a positive individual is more (or even less) likely to test positive on subsequent days.

I've searched through a number of alternative questions, but nothing seems to quite capture this issue. I could try repeated-measures analyses, but the latency and return-to-baseline should both swamp my effect if they are sufficiently long. Further, it would be ideal to actually know when the effect is observed.

A sample data set:

make_obs <- function(group, day){
sapply(1:length(group), function(idx){
rbinom(1, 1, day_prob[[group[idx]]][day[idx]] )
})
}

set.seed(12345)
example_data <-
tibble(
Group = rep(c("A", "B"), each = 50)
, Ind = 1:100
, Day = 1
) %>%
complete(nesting(Group, Ind), Day = 1:15) %>%
mutate(
Obs = make_obs(Group, Day)
)


Plotting to show the observed data:

example_data %>%
group_by(Group, Day) %>%
summarise(Prop_pos = mean(Obs)) %>%
ungroup() %>%
ggplot(aes(x = Day
, y = Prop_pos
, col = Group)) +
geom_line() +
scale_color_brewer(palette = "Dark2")


• The chi-squared call you use in the simulated data - does that represent an actual model of the outcome, or is it just a function that produces good-looking results? If the latter is the case, would you be able to describe the response to treatment as a parametric curve and/or differential equation? Vaguely reminds me of PK/PD (pharmacokinetic/pharmacodynamic) models - is there a connection? Jun 3 '19 at 19:58
• The chi-square model was simply handy to produce roughly the pattern that I was looking for (after lots of tinkering with df etc.) -- it does not describe any particular underlying result. I am fine using any underlying model; the output data are probably the more relevant, but I thought that it would be instructive to include the generative model. There is no direct connection of PK/PD to the model I have in mind, though it could explain the pattern in similar cases. Jun 3 '19 at 20:52
• Let me understand one thing: in the model you wrote in R, it seems to me that the same individual has probability p(day, group) of responding to stimulus at each day, independent on other days events. This is not what happens on disease onset models (survival models), where every individual who shows the stimulus is then removed from the sample. Explain better the ideas behind your model, please. Nov 8 '19 at 12:23
• @carlo, In this model, the probability of testing positive (presenting with symptoms) is independent of preceding days of testing. I understand that this may not match expectations in all cases, but in the abstract, it is plausible (recall that the example applies to any case of positive-negative testing over time, not just to disease scenarios). A survival model is an intriguing possibility, though I am not sure that it captures the relative risk that I am looking for here. Nov 8 '19 at 16:02
• For a plausible scenario that matches this independence, we can consider a priming event (e.g., advertisement) that increases the probability of clicking on a banner ad in the future. We are comparing the effect of two different priming events (e.g., email and phone call). While there is an argument that a single click on a banner ad may indicate an individual that is more (or less) likely to click on a subsequent day, it is not unreasonable to treat the days as independent (especially if the banner ads change each day). Survival analysis (to first click) is still an interesting possibility. Nov 8 '19 at 16:09

In biostatistics, time-dependent responses are often inspected in terms of the area under the curve. This is not to be confused to be the AUC of the ROC which is a classification metric. For instance, in an oral glucose challenge, blood is drawn hourly to see the glucose concentration over a period of time. Because diabetes and its complications are caused by cumulative concentration of blood sugar, the AUC of 2hOGTT blood sugar describes how effectively the body produces insulin.

https://www.researchgate.net/figure/Comparison-of-the-area-under-the-curve-AUC-values-of-blood-glucose-levels-obtained-from_fig1_273472980

Here the AUC reduces an overly complex question into a simple one: does the intervention reduce the time-averaged risk in one group versus another? In doing this, we are able to exonerate ourselves from making strong and possibly incorrect assumptions about the dependence structure. With imbalanced data, the curves can be estimated using splines over a reasonable range (restricted mean). For instance, if the curve in Group A does not "peak" as largely as that in Group B, but remains elevated over a long period of time, none of the modeling approaches proposed above will correctly estimate the standard deviation of the conditional response, and yet AUC simply addresses why Group B might be preferred (a less sustained response).

• +1 Very often, simpler is better.
– EdM
Nov 14 '19 at 19:47

To tackle your problem you could fit a logistic regression playing with your covariates to capture the time and group effects conjointly.

library(Hmisc)
#>
#> Attaching package: 'Hmisc'
#> The following objects are masked from 'package:base':
#>
#>     format.pval, units
library(rms)
#>
#> Attaching package: 'SparseM'
#> The following object is masked from 'package:base':
#>
#>     backsolve
library(broom)
library(modelr)
#>
#> Attaching package: 'modelr'
#> The following object is masked from 'package:broom':
#>
#>     bootstrap
library(yardstick)
#> For binary classification, the first factor level is assumed to be the event.
#> Set the global option yardstick.event_first to FALSE to change this.
#>
#> Attaching package: 'yardstick'
#> The following objects are masked from 'package:modelr':
#>
#>     mae, mape, rmse
library(tidyverse)

make_obs <- function(group, day){
sapply(1:length(group), function(idx){
rbinom(1, 1, day_prob[[group[idx]]][day[idx]] )
})
}

day_prob <-
list(
A = c(0.06, 0.06, pmax(dchisq(1:13, 3)*1, 0.06))
, B = c(0.06, 0.06, pmax(dchisq(seq(1, 25, 2), 5)*6, 0.06))
)

set.seed(12345)
example_data <-
tibble(
Group = rep(c("A", "B"), each = 50)
, Ind = 1:100
, Day = 1
) %>%
complete(nesting(Group, Ind), Day = 1:15) %>%
mutate(
Obs = make_obs(Group, Day)
)

example_data %>%
group_by(Group, Day) %>%
summarise(prop = mean(Obs)) %>%
ggplot(aes(Day, prop, group=Group)) +
geom_line()


By examining the grouped plot by day and group form the data you generated, you can see that there is a highly non linear effect from day and group interactions.

I propose the following 4 model to study how the effect of Day and Group onto Obs. In models 3 and 4 I fit a non-linear transformation to Days, a restricted cubic spline to attempt to capture the effect mentioned in the previous plot.

df <- example_data

md_log <- glm(Obs ~ Day + Group, data = df, family = 'binomial')
md_log2 <- glm(Obs ~ Day * Group, data = df, family = 'binomial')
md_log3 <- glm(Obs ~ rcs(Day) + Group, data = df, family = 'binomial')
md_log4 <- glm(Obs ~ rcs(Day) * Group, data = df, family = 'binomial')

gather_predictions(df, md_log, md_log2, md_log3, md_log4, type='response') %>%
group_by(model, Group, Day) %>%
summarise(prop = mean(Obs),
fitted = mean(pred)
) %>%
ggplot(aes(x = Day, group=Group)) +
geom_line(aes(y = prop)) +
geom_line(aes(y = fitted, group=model, color = model), lty='dashed') +
facet_wrap(~Group, ncol =1)


I the gather the predictions from each model and aggregate them to have a visual check of fit. As you can see, models 3 and 4 capture the effect you are interested in.

md_frame <-
tibble(number= seq(4), model = list(md_log, md_log2, md_log3, md_log4))

get_auc <- . %>% augment( . , df, type.predict='response') %>%
mutate(Obs = as.factor(Obs)) %>%
yardstick::roc_auc(., truth=Obs, .fitted) %>%
pull(.estimate)

md_frame %>%
mutate(auc = map_dbl(model, get_auc)) %>%
mutate(other_metrics = map(model, glance)) %>%
unnest(other_metrics) %>%
select(model_number = number, auc, AIC, deviance, BIC, logLik) -> metrics

metrics %>%
pivot_longer(-model_number, names_to = 'metric_name', values_to='values') %>%
mutate(model_number = paste0('model ', model_number)) %>%
ggplot(aes(model_number, values)) +
geom_point(aes(fill=model_number), size = 4, pch=21, show.legend = F) +
facet_wrap(~metric_name, scales = 'free') +
labs(title = 'Model metrics', x = NULL, y=NULL)


Now, to check the metrics for accessing the quality of fit from each of these models also indicate that models 3 and 4 capture the Days effect. By inspecting the models you can check the coefficients.

Now, if you intend to perform some statistical inferences with these models, there are better approaches when it comes estimating nested variances and time dependent effects, so a longitudinal approach could be advantageous here.

Created on 2019-11-14 by the reprex package (v0.3.0)

• A decent first pass would be a generalized ARMA model followed by the analysis of impulse response functions (IRFs).
• Interesting IRF measurements: peak, delay until peak, halftime of peak decay, cumulative response (AUC, area under the IRF's curve, as mentioned by @AdamO).
• In R VGAM package seems to do the job.
• Construction of confidence intervals for chosen statistics might need a moderate amount of additional coding.

Details:

If I understand correctly, each period each individual in you sample is exposed to some random shock (you call it stimulus, e.g. seeing a banner of random appearance, shaking hands with the contaminated person etc.). Upon a shock the individual may react (click on the banner, contract a desease) or not to the shock. It is reasonable to assume that continuous exposure to shock for several period increases the likelihood of reaction (the degree of this relation is to be inferred from the data). And you are seeking to gauge how this sensitivity to shocks changes with some treatment.

Assuming you have binary outcome data, i.e. $$y_{i,t} \in \{0,1\}$$.

1. The model would look like: $$\begin{cases} p_{y^k_{i,t}=1}& = & f(z^k_{i,t})\\ z^k_{i,t} &= & \alpha_0 + \sum_{s=1}^m \alpha^k_s z_{i,t-s} + \epsilon_{i,t} + \sum_{i=1}^{n} \beta^k_s \epsilon_{i,t-s} \\ \epsilon_{i,t} &\sim_{iid} & N(0,\sigma^2) \quad \forall i,t, \end{cases}$$

where $$f: \mathbb{R} \to [0,1]$$ is some function mapping the hidden variable $$z_{i,t}$$ into the probability of reacting to stimulus (e.g. logit), and $$k\in\{\text{treated},\text{nontreated}\}$$ is the index of the group and $$i$$ is the index of the individual. Less explicit assumption is that the shock $$\epsilon$$ hits everyone and at all periods evenly (i.e. it's variance does not depend on treatment nor on the identity $$i$$), which may be relaxed if data allows and/or the initial formulation is deemed insufficent.

Estimation may be performed by likelihood maximization or Bayesian methods.

2. Upon estimation, two impulse response functions, derived from estimates $$\{a^k_s, b^k_s\}_s$$ of $$\{\alpha^k_s, \beta^k_s\}_s$$ can be obtained, one per $$k$$. You can compare their peaks, delay of the peak, cumulative response, half-time after the peak etc (for more ideas one may check the litterature on linear time-invariant systems). Judgements about statistical significance can be based on the Bayesian or bootstrapped confindence intervals.

3. Quick search showed some examples of packages able to handle this kind of model gave VGAM package, wich seems to be well maintained and even accompanied by a book.