Modeling Target Variable based on Days Since Last Engagement I have a data frame of engagement data for different sites for every day like this:




SITE_ID
Date
Engagement A Count
Engagement B Count
Target Variable




1
1/1/2022
0
0
1


1
1/2/2022
0
0
0


1
1/3/2022
1
0
0


1
1/4/2022
0
0
0


1
1/5/2022
0
1
0


1
1/6/2022
0
0
1


...
...
...
...
...


2
1/1/2022
0
0
0


2
1/2/2022
0
0
0


2
1/3/2022
1
0
0


2
1/4/2022
0
1
2


2
1/5/2022
0
0
0


2
1/6/2022
0
0
0


...
...
...
...
...




I was wondering what would be the best way to longitudinally model the effect of engagements on the target variable (count data, assume Poisson). I have tried calculating indicator columns of whether each type engagement occurred in the last N days and throwing those into a Poisson Linear Mixed Model which achieved significant results. The drawback with this approach I see is that only the most recent engagement  gets "credit" if the target variable equals 1 when I think all engagements up until that target variable equaling 1 should be credited (with a decaying effect based on how long it's been since the engagement). Another drawback is that choosing N can be quite arbitrary and "hand-wavey" to explain.
I also considered throwing instead into the model the cumulative total of each engagements up until each day, but I don't like this approach because it doesn't take into account of how long ago those engagements took place. But I am willing to listen to other opinions if this is the best approach.
What I am thinking is having columns like "Days since kth most recent A Engagement", and similarly for B engagements. But then I get NA values for sites on days that have not yet had an engagement. So I was thinking of transforming this variable via $e^{-x}$ or with a base that appropriately quantifies the decaying correlation of time since engagements with the target variable. That way engagements that occurred a long time ago have a near zero value and recent engagements are closer to 1. Also I can then impute NAs with 0 which makes sense. But now I have a bunch of columns like ["Days since kth most recent A Engagement (k=1)", "Days since kth most recent A Engagement (k=2)", ..., "Days since kth most recent A Engagement (k=14)"] that are bounded between 0 and 1. Do I just throw all of these into a linear model now? If this is a valid approach, could I make any conclusions about the correlation of types of engagements with the target variable besides directionally?
I'm kind of stuck on what to do because it seems every approach seems to have a serious drawback or unnatural way of capturing the true essence of the data. I am curious what the minds of stats Stack Exchange have to say. And also, this is truly observational data with no random assignment of engagements, so it is very possible that more engagements tended to be in sites with a higher amount of target variable equaling 1 by correlation and selection bias. There are other site level features that I omitted for simplicity.

There are lots of ways you could model this type of data, but I will confine my attention here to models that are roughly in accordance with what you are proposing with your idea of having the engagements given an exponentially decaying effect on the target variable.
Your proposed method of assuming exponential decay in relevance might be a reasonable way to proceed, but your present form assumes a fixed rate of decay rather than estimating this from the data.  To be effective, your model would need to estimate the rate of decay from the data, which would yield a nonlinear count model.  Suppose we let $\{ A_t \}$, $\{ B_t \}$, $\{ S_t \}$ and $\{ Y_t \}$ denote the time-series of engagement values for A and B, sites (categorical variable), and the target values respectively.  Then you would want to use some kind of GLM equation like:
$$\mathbb{E}(g(Y_t)) \equiv \mu_t 
= \beta_A \sum_{k=0}^\infty \exp(-\theta_A A_{t-k}) + \beta_B \sum_{k=0}^\infty \exp(-\theta_B B_{t-k}) + \sum_{s=1}^{K-1} \beta_{k} \mathbb{I} (S_t = s),$$
where in practice we would cut off the sums at some sufficiently large value $N$.  (The function $g$ is the link function in the GLM.)  You are correct that there is some arbitrariness in cutting this off at a finite value $N$, but if this value is sufficiently large then the effect of engagements further back in time should be minimal.  (You can check the upper bound on earlier effects as a post-hoc inference once you've estimated your coefficients, so you should be able to satisfy yourself that you have enough model terms included.)  There is also a drawback in assuming that unobserved earlier values are no-engagement, but again, this effect can be limited by noting the upper bounds on the exponential decay.
In any case, if you decide to go with something like this, it is a nonlinear GLM equation where the rate parameters in the exponentials cannot be "linearised".  Consequently, you would need to use a function like the gnlr function in the gnlm package to fit the model.  You should use a negative binomial count model rather than a Poisson count model, because the negative binomial is a two-parameter model that allows you to correctly estimate the dispersion.
If you were to successfully implement a model like this, it would allow you to estimate the effects $\beta_A$ and $\beta_B$ from the two types of engagements to the target variable and the rates of decay $\theta_A$ and $\theta_B$ that determine the degree to which older engagements lose relevance to the target variable.  You would also be able to estimate the effect of the site (categorical variable) on the target variable.  You would need to augment your analysis with appropriate diagnostic plots (e.g., added variable plots) to see if your posited model looks reasonable from the data.  Ultimately, you would be looking to get a model where the fit is reasonable and the diagnostic plots do not falsify your model assumptions.
