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")


 A: 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).
A: To tackle your problem you could fit a logistic regression playing with your covariates to capture the time and group effects conjointly.   
library(Hmisc)
#> Loading required package: lattice
#> Loading required package: survival
#> Loading required package: Formula
#> Loading required package: ggplot2
#> 
#> Attaching package: 'Hmisc'
#> The following objects are masked from 'package:base':
#> 
#>     format.pval, units
library(rms)
#> Loading required package: SparseM
#> 
#> 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: Answer:


*

*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\}$.


*

*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.

*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.

*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.
