In the last two blogs, we went through linear regression and generalized linear regression in search of the best solution to the covid data by relaxing the assumptions of linear regression and OLS. Based on what have done so far, our conclusion is negative binomial regression fit the data best. However, still, we have this problem that for some states with high vaccination rates, the increasing cases are over-estimated, and a few state such as FL and KY, are more like outliers to the model. Today, we are going to fixed these issues with the mixed-effects model, specifically, the generalized linear mixed-effects model(GLMM), which is very popular in biostatistics and econometrics, especially useful for economists when modeling cross-country macro data.

A mixed-effects model contains both fixed effects and random effects. They are particularly useful in settings where measurements are made on clusters of related statistical units. There are several ways people usually choose to deal with this type of inherently “lumpy” data.

1. Pooled Data Regression

One can run a regular regression on $y$ versus $x$ ignoring the group information, like what we have done in the last two blogs.

It works but:

• It throws away useful information.
• Linear Regression assumes observations are independent here they are not.
• The $\alpha$ and $\beta$ coefficients are unbiased estimates (not bad).
• The confidence intervals will be too optimistic.

2. Unpooled Data Regression

Instead one can run a regression for each of the groups.

• One needs to estimate too many coefficients (51 $\alpha$s and 51 $\beta$s in our case).
• This results in unbiased, but very noisy estimates if some groups have few observations.

3. Pooled Slope, Unpooled Intercept

We can instead fit a regression with a dummy variable for each group.

$$
y_i = \alpha + \beta x + \rho^T Z_i
$$
where $\rho$ is now a vector with as many coefficients as groups, and $Z$ is a vector of dummies, one column per group (we must leave one group out to avoid co-linearity of Z).

This method:

• Pools information for the slope $\beta$.
• Treats each group intercept independently. If some group $g$ has few examples $\rho_g$ will be noisy.

4. Penalized Linear Regression

To pool over the different groups, we can introduce a regression penalty to penalize model coefficients to reduce the model degrees of freedom.

• Theoretically, the correct penalty is proportional to the ration of variances for the observations $\sigma_Y^2$ and the random effects $\sigma_\rho^2$.
• Because there is no penalty for $\alpha$, the mean over all the population will be matched exactly.

5. Linear Mixed-effect Model

For people who are familiar with penalized regressions and regularizations like lasso and rigid, the estimation of mixed models would be very straightforward as essentially, they are the same. Same result can be obtained with less work by fitting to a Linear Mixed Effect Model with statsmodels.mixed_linear module.

6. Generalized Linear Mixed-effects Model (GLMM)

Today, as we already know that negative binomial distribution fits the covid data better, we will skip the linear mixed model and start with the non-linear GLMM. Since the model is relatively complex, we will be using a much cooler statistical inference method – Bayesian.

Preliminaries

1
2
3
4
5
6
7
8
9
10

import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import seaborn as sns
import pymc3 as pm
import arviz as az
import sys
import covid_analysis as covid

As usual, let’s prepare the data. Hospital admissions and vaccinations are from the CDC website, and the most up-to-date populations by state are from the US Census Bureau.

1

data_dir="../data"

1
2
3
4
5
6
7
8
9
10
11

hospitalizations=pd.read_csv(f"{data_dir}/covid_hospitalizations.csv",parse_dates=["date"])
vaccinations=pd.read_csv(f"{data_dir}/covid_vaccinations.csv",parse_dates=["date"])
population=pd.read_csv(f"{data_dir}/population.csv")

event="admissions"
data=hospitalizations

data=data.merge(vaccinations,on=["date","state"])
data=data.merge(population,on="state")
data["vaccinated"]=data["vaccinated"]/data["population"]
data.head()

Split the train and test data,

1
2
3
4
5
6
7
8
9

# train dates
start_train="2021-07-15"
end_train="2021-08-15"
train_data=data[(data['date'] >= start_train) & (data['date'] <= end_train)].copy()
date0=train_data["date"].min()
# test dates
test_period=7 # days
test_end=pd.Timestamp(end_train)+pd.DateOffset(days=test_period)
test_data=data[(data["date"]>end_train) & (data["date"]<=test_end)].copy()

Here we write a function to generate our desired matrix. This time, we will need an extra step to label data by state.

1
2
3
4
5
6
7
8
9
10

def define_variables(data,date0):
        N=len(data)
        T=(data["date"]-date0).dt.days/30
        X=pd.DataFrame({"const":np.ones(N)})
        X["vaccinated"]=data["vaccinated"].values
        X["T"]=T.values
        X["T2"]=T.values**2
        P=data["population"].values 
        Z=pd.get_dummies(data["state"])        
        return X,Z,P

1
2
3

X_train,Z_train,P_train=define_variables(train_data,date0)
Y_train=train_data[event]
G_train=Z_train.values.argmax(axis=1)

We have 51 regions in total,

1
2
3

# number of states/regions
K = G_train.max()+1
K

51

1
2
3

# all states/regions
states=Z_train.columns
states

Index(['AK', 'AL', 'AR', 'AZ', 'CA', 'CO', 'CT', 'DC', 'DE', 'FL', 'GA', 'HI',
       'IA', 'ID', 'IL', 'IN', 'KS', 'KY', 'LA', 'MA', 'MD', 'ME', 'MI', 'MN',
       'MO', 'MS', 'MT', 'NC', 'ND', 'NE', 'NH', 'NJ', 'NM', 'NV', 'NY', 'OH',
       'OK', 'OR', 'PA', 'RI', 'SC', 'SD', 'TN', 'TX', 'UT', 'VA', 'VT', 'WA',
       'WI', 'WV', 'WY'],
      dtype='object')

Mixed Negative Binomial Regression Model

Model Definition

We have data from multiple dates, so it makes sense to incorporate Time as a extra variable.

\begin{eqnarray}
\eta_{i,j} =& \beta_0 + \beta_1 x_{i,j} + \beta_2 t_{j} + \beta_3 t_{j}^2 + \rho_i\newline
\rho_i \sim & N(0,\sigma_{\rho}^2) \newline
\lambda_{i,j} =& e^{\eta_{i,j}} \newline
y_{i,j} \sim & \text{NB}(y;\lambda_{i,j} P_i,\alpha)
\end{eqnarray}
with priors
\begin{eqnarray}
\beta_0 \sim & N( \bar{y},1) \newline
\beta_i \sim & N(0,1) \newline
\log \sigma_\rho \sim & N(0,1) \newline
\log \alpha \sim & N(0,1)\newline
\end{eqnarray}
where

• Each index $i$ represents a state in the US.
• Each index $j$ represents a particular date.
• $y_{i,j}$ is the number of hospital admissions on state $i$ at date $j$.
• the random effects (random intercept) $\rho_i$ for different states $i$ are independent and have a normal distribution. It will be fit to the data.
• $x_{i,j}$ is the percentage of the state population fully vaccinated on state $i$ at date $j$.
• $t_j$ is the number of days elapsed since the beginning of the training period.
• $P_i$ is the population of state $i$.
• $\bar{y}$ is log of the average per person incidence on the event in the US Population.
• the parameters and $\beta_1,\beta_2,\beta_3$ will be fitted to the observed data to maximize agreement with the model.
• the parameter $\sigma_\beta$ is our prior uncertainty about the value of $\beta_i$, it is set to 1.

To estimate the model coefficients, we will go for full Bayesian.

Bayesian Approach

Unlike the classical frequentist methods, in Bayesian analysis, a parameter is summarized by an entire distribution of values instead of one fixed value, and it provides a natural and principled way of combining prior information or expert knowledge with the data observed.

Both Bayesian methods and classical methods have advantages and disadvantages, and there are some similarities. When the sample size is large, Bayesian inference often provides results for parametric models that are very similar to the results produced by frequentist methods. Some advantages to using Bayesian analysis include the following:

1. It provides a natural and principled way of combining prior information with data, within a solid decision theoretical framework. You can incorporate past information about a parameter and form a prior distribution for future analysis.

2. It provides interpretable answers based on parameter distributions, such as “the true parameter has a probability of 0.95 of falling in a 95% credible interval.”

3. It provides a convenient setting for a wide range of complex models, such as the GLMM model we are trying to build but not easy to get an analytical solution. MCMC, along with other numerical methods, makes computations tractable for virtually all parametric models.

The disadvantages of Bayesian are quite obvious as well:

1. It does not tell you how to select a prior. There is no correct way to choose a prior. Bayesian inferences require skills to translate subjective prior beliefs into a mathematically formulated prior.

2. It often comes with a high computational cost, especially in models with a large number of parameters.

To solve our model with Bayesian approach, we have to specify the distributions of the parameters/variables. The good thing is that there are convenient functions provided by pymc3.

1
2
3

Y_bar=np.log(Y_train.sum()/P_train.sum())
beta0=np.array([Y_bar,0,0,0,])
sigma_beta=np.array([1,1,1,1])

I choose $N(0,1)$ for $\beta$s , and $lognormal(0,1)$ for $\sigma_\rho$s just to speed up the convergence.

1
2
3
4
5
6
7
8
9
10
11
12

with pm.Model() as mod:
    beta=pm.Normal("beta",mu=beta0,
                   sigma=sigma_beta,
                   shape=(len(beta0))
                  )
    lsigma_group=pm.Normal("lsigma_group",mu=0,sigma=1)
    sigma_group=pm.Deterministic("sigma_group",np.exp(lsigma_group))
    rho=pm.Normal("rho",mu=0.0,sigma=sigma_group,shape=(K)) # a vector of K values
    eta=pm.math.dot(X_train.values,beta)+rho[G_train]
    y_hat=pm.Deterministic("eta",pm.math.exp(eta))
    a=pm.Lognormal("alpha",mu=0,sigma=1)
    y=pm.NegativeBinomial("y",mu=y_hat*P_train,alpha=a,observed=Y_train)

Then we can start the Monte Carlo sampling process. 4 chains, of which 1000 samples are generated. That means there will be 4000 samples for each of the parameters. Be patient, the computation can take a while. (Honestly, sometimes HOURS if using personal laptops)

Monte Carlo Sampling

1
2

with mod:
    trace = pm.sample(1000,chains=4,cores=4, tune=500)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [alpha, rho, lsigma_group, beta]
Sampling 4 chains, 0 divergences: 100%|█████████████████████████████████████████| 6000/6000 [04:36<00:00, 21.72draws/s]
The acceptance probability does not match the target. It is 0.8915155058428972, but should be close to 0.8. Try to increase the number of tuning steps.
The rhat statistic is larger than 1.05 for some parameters. This indicates slight problems during sampling.
The estimated number of effective samples is smaller than 200 for some parameters.

Regression Coefficients

One of the advantages of Bayesian is that a parameter is summarized by a distribution of values instead of a fixed value, which definately helps us get a better understanding of the things going on, like how high the uncertainty is, how much confidence the model has in the estimations.

1
2
3
4

with mod:
    summary=az.summary(trace, var_names=["beta","sigma_group"], fmt="wide",round_to=2)
summary.index=["intercept","vaccinated","T","T2","sigma_group"]
summary

We can see the distribution of parameters very clearly,

1
2

with mod:
    az.plot_trace(trace,var_names=["beta","sigma_group","alpha"])

Random Effects

We also get a distribution for the random, state specific, effects.

1
2
3

with mod:
    random_effects=az.summary(trace, var_names=["rho"], fmt="wide",round_to=2)
random_effects.index=states

1

random_effects.head(9)

Comparing to the other states, states such as FL, KY which are seen as outliers in our previous models have much larger random effects.

1

random_effects.loc[['FL','KY']]

In Sample Predictions

Now, let’s make some predictions and see how the mixed model works.

1
2
3
4
5

def select_data(data,periods):
    dates=data["date"].unique()
    selected_dates=covid.select_dates(dates,periods)
    used_data=data.merge(selected_dates,on="date")
    return used_data.groupby("date")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

def plot_result(data,norm,ax):
    P=data["population"]
    m=ax.scatter(data["vaccinated"],data[event]/data["population"]*100_000,
                c=P,norm=norm,cmap="Greys",label=event)
    
    ax.errorbar(data["vaccinated"],data["y_pred"],yerr=data["y_std"],fmt="D",alpha=0.25,
               label="predicted")#,c=P,norm=norm,cmap="Blues",label="predicted")
    ax.set_xlabel("vaccinated")
    ax.set_ylabel(f"{event} per 100k")
    ax.legend()
    return m

def plot_time_facets(data,event):
    # Create two subplots and unpack the output array immediately
    fig, axes = plt.subplots(2, 2, figsize=(14,12),sharey=True,sharex=True)
    P=data["population"]
    norm=matplotlib.colors.LogNorm(vmin=100_000, vmax=P.max(), clip=False)
    count=0
    for date,group in select_data(data,4):
        col=count %2 
        row=count //2
        ax=axes[row][col]
        ax.set_title(date.strftime("%Y-%m-%d"))
        im=plot_result(group,norm,ax)
        count+=1
    cbar = fig.colorbar(im, ax=axes.ravel().tolist(), shrink=0.5,label="population")

1
2
3

with pm.Model() as mod:
    y_pred=pm.NegativeBinomial("y_pred",mu=y_hat*100_000,alpha=a,shape=(len(X_train)))
    posterior_predictive = pm.sample_posterior_predictive(trace, var_names=["y_pred"])

100%|█████████████████████████████████████████████████████████████████████████████| 4000/4000 [00:05<00:00, 798.51it/s]

Remember, there are 4000 sample predictions for each of the records, I won’t plot out all of them, but will do their means with one standard deviation of uncertainty. The error bars would make prefect graphical representations.

1
2

predicted=posterior_predictive["y_pred"]
predicted.shape,train_data.shape

((4000, 1632), (1632, 6))

1
2
3
4
5
6

y_pred_mean=predicted.mean(axis=0)
y_pred_std=predicted.std(axis=0)

data=train_data.copy()
data["y_pred"]=y_pred_mean
data["y_std"]=y_pred_std

Looks like the predictions are closer to the actuals since we incorporate the random effects. Hooray!

1

plot_time_facets(data,event)

Differences comparing to linear regression

Random Intercepts

Let’s take a further look by state to better understand what I meant by random effects/intercepts.

Below is all the observed data points marked by state. We can see that roughly the slopes of different groups are close, while the intercepts can vary.

In case you wonder why the slopes are upward, that’s due to the passage of time, remember time is our explanatory variable as well, but we can only plot two dimensions.

1
2
3
4
5
6
7
8

fig, ax = plt.subplots()
P=data["population"]
norm=matplotlib.colors.LogNorm(vmin=100_000, vmax=P.max(), clip=False)
m=ax.scatter(data["vaccinated"],data[event]/data["population"]*100_000,
            c=P,norm=norm,cmap="Greys",label=event, s = 5)
ax.set_xlabel("vaccinated")
ax.set_ylabel(f"{event} per 100k by state")
ax.legend()

Let’s select three states in the middle to plot their predictions against actuals.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

def plot_states(data, states, var):
    data = data[data['state'].isin(states)]
    data = data.sort_values('vaccinated')
    x = data[var]
    y = data["y_pred"]
    std = data["y_std"]
    cmp=cm.get_cmap("tab10")
    plt.xlim([0.4475,0.466])
    plt.scatter(x,y, c='black',s=15)
    colors = np.array(['tab:blue','tab:green','tab:orange','tab:red'])
    for i in range(len(states)):
        s = states[i]
        ds = data[data['state']==s]
        x = ds[var]
        y = ds["y_pred"]
        std = ds["y_std"]
        color = colors[i]
        plt.plot(x,y,'k-',color=color,linewidth=2.0, label=f"predicted {s}")
        plt.fill_between(x,y+std,y-std,alpha=0.15,color=color)
        plt.xlabel("vaccinated")
        plt.ylabel(f"{event} per 100k by state")
        plt.legend(loc="upper left")

This shows clearly how the mixed model leverages all the data points to estimate the correlation (slope) between admissions and vaccinations, but at the same time allows states to have different starts (intercept) due to fundamental differences which we consider are random.

1

plot_states(data,['KY','AK','KS'],"vaccinated")

Variance

The Bayesian approach gives the distributions of the predictions, so we don’t need to rack our brains for the analytical solution of the prediction band, instead, we can simply plot it out with the samples generated by pymc3. Here I used one unit of the sample standard deviation.

As we can see, the variance of prediction is larger with greater admissions, meaning they are positively correlated, which is not what we’ve seen with the linear regression that the variance of dependent variable is always constant. So using GLMM based on a negative binomial distribution allows the magnitude of the variance to change as a function of the predicted value.

1
2
3
4
5
6
7
8
9

def plot_state(data, state, color):
    data = data[data['state']==state]
    fig, ax = plt.subplots()
    ax.scatter(data["vaccinated"],data[event]/data["population"]*100_000,c='black',s=15,label=event)    
    ax.errorbar(data["vaccinated"],data["y_pred"],yerr=data["y_std"],fmt="o",alpha=0.25,c=color,
               label="predicted")
    ax.set_xlabel("vaccinated")
    ax.set_ylabel(f"{event} per 100k")
    ax.legend()

1

plot_state(data,'KY','tab:blue')

1

plot_state(data,'FL','tab:green')

Predicted Dependence on vaccination Rate

One of the things that we are interested is that how does the vaccination rate impact the admissions?

To see this better, I fixed the time at Aug.15, and let the vaccination rate to vary from 20% to 100%.

1
2
3
4
5
6

# given current time point, generate test set for vaccination rate from 20% to 100%
V_N=201
V=np.linspace(0.2,1,V_N)
T=np.ones(V_N)
T2=T**2
X2_test=np.c_[np.ones(T_N),V,T,T2]

1
2
3
4

with pm.Model() as mod:
    eta2_pred=pm.math.dot(X2_test,beta)
    y_pred2=pm.Deterministic("y_pred2",pm.math.exp(eta2_pred)*100_000)
    posterior_predictive = pm.sample_posterior_predictive(trace, var_names=["y_pred2"])

100%|████████████████████████████████████████████████████████████████████████████| 4000/4000 [00:01<00:00, 3891.10it/s]

1
2
3
4

# generate predicted mean and standard deviation
predicted=posterior_predictive["y_pred2"]
y_pred_mean=predicted.mean(axis=0)
y_pred_std=predicted.std(axis=0)

We can see that vaccinations and admissions are negatively correlated. And when vaccination rate is low, the prediction can be quite uncertain.

1
2
3
4
5

plt.plot(V,y_pred_mean,"k--",linewidth=1)
plt.fill_between(V,y_pred_mean+y_pred_std,y_pred_mean-y_pred_std,alpha=0.07,color="k")
plt.fill_between(V,y_pred_mean+2*y_pred_std,y_pred_mean-2*y_pred_std,alpha=0.07,color="k")
plt.xlabel("Vaccinated")
plt.ylabel(f"{event} per 100k")

Predicted Time Evolution

The other thing that we are most interested in is that how the third wave of covid19 is going to evolve. How many admissions are there going to be given the current vaccination rate? So I use the GLMM model to project the next two months from Aug.15th on.

1
2
3

# given the average vaccination rates
mean_vac=(train_data["vaccinated"]*train_data["population"]).sum()/train_data["population"].sum()
mean_vac

0.48852180987716365

1
2
3
4
5

# generate test set for the past month and the next 2 months
T_N=201
T=np.linspace(0,3,T_N)
T2=T**2
X1_test=np.c_[np.ones(T_N),mean_vac*np.ones(T_N),T,T2]

1
2
3
4

with pm.Model() as mod:
    eta1_pred=pm.math.dot(X1_test,beta)
    y_pred1=pm.Deterministic("y_pred1",pm.math.exp(eta1_pred)*100_000)
    posterior_predictive = pm.sample_posterior_predictive(trace, var_names=["y_pred1"])

100%|████████████████████████████████████████████████████████████████████████████| 4000/4000 [00:01<00:00, 3831.47it/s]

1
2
3
4

# generate predicted mean and standard deviation
predicted=posterior_predictive["y_pred1"]
y_pred_mean=predicted.mean(axis=0)
y_pred_std=predicted.std(axis=0)

The prediction tells us that the admission will peak in one and half month from Aug.15th, which is by the end of September. However, I would be careful as the model have low confidence when projecting beyond one month. Look at that wide confidential interval band after T=2. Basically, the reality can fall any where between the curve flattening out as early as the start of September to that we won’t be able to see a peak until mid October. But overall I am very optimistic as there’s at least 80% of chance that we’ll see it peaking before mid October. Eventually, we will get there. So hang in there, everyone!!

1
2
3
4
5
6
7
8

T0=T[T<1]
plt.plot(T0,y_pred_mean[T<1],"k",linewidth=4,label="fitted")
plt.plot(T,y_pred_mean,"k--",linewidth=1,label="extrapolated")
plt.fill_between(T,y_pred_mean+y_pred_std,y_pred_mean-y_pred_std,alpha=0.07,color="k")
plt.fill_between(T,y_pred_mean+2*y_pred_std,y_pred_mean-2*y_pred_std,alpha=0.07,color="k",)
plt.xlabel("T (months)")
plt.ylabel(f"{event} per 100k")
plt.legend(loc="lower right")

Officially we are done, cheers! I hope you find this project inspiring. And if you read it all the way through - congrats, you are an incredibly patient person, have a cookie!

In the last blog, we try to use linear regressions to model hospital admissions and vaccinations. The findings were the linear model tend to over-estimating how fast the admission rate is increasing on the higher range of vaccination rates, and residual variance is negatively correlated with vaccinations rates, suggesting heteroskedasticity. But it’s easy-peasy for Generalized linear models to fix those issues. GLM generalizes linear regression, 1. by allowing the linear model to be related to the response variable via a link function, 2. by allowing for the response variable to have an error distribution other than the normal distribution so that the magnitude of the variance of each measurement can change as a function of its predicted value.

Preliminaries

1
2
3
4
5
6
7
8
9
10

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.colors
import seaborn as sns
import statsmodels.api as sm
from statsmodels.discrete.discrete_model import Poisson
from statsmodels.discrete.discrete_model import NegativeBinomial
import sys
import covid_analysis as covid

The data I use includes hospital admissions and vaccinations from the CDC website, and the most up-to-date populations by state from the US Census Bureau.

1

data_dir="../data"

1
2
3
4
5
6
7
8
9
10
11

hospitalizations=pd.read_csv(f"{data_dir}/covid_hospitalizations.csv",parse_dates=["date"])
vaccinations=pd.read_csv(f"{data_dir}/covid_vaccinations.csv",parse_dates=["date"])
population=pd.read_csv(f"{data_dir}/population.csv")

event="admissions"
data=hospitalizations

data=data.merge(vaccinations,on=["date","state"])
data=data.merge(population,on="state")
data["vaccinated"]=data["vaccinated"]/data["population"]
data.head()

Poisson Regression

Like what we did with the linear regression, we can choose either to use the Poisson Regression model from statsmodels, or the GLM from the same lib. The two will give us the exact same estimations and statistics.

Model Definition

We have data from multiple dates, so it makes sense to incorporate Time as an extra variable.

\begin{eqnarray}
&\eta_i &=c + \beta_1 x_i + \beta_2 t_i + \beta_3 t_i^2 \newline
\mathbb{E}(y_i|x_i) =& \lambda_iP_i &= e^{\eta_i }P_i \newline
&p(y_i|x_i) & = \text{Poisson}(y;\lambda_i P_i)
\end{eqnarray}
where

• Each observation $i$ represents a single state at one particular date.
• $y_i$ is the number of hospital admissions on that state.
• $x_i$ is the percentage of the state population fully vaccinated that state.
• $t_i$ is the number of days elapsed since the beginning of the training period.
• $P_i$ is the population of the state.
• the parameters $c$ and $\beta_1,\beta_2,\beta_3$ will be fitted to the observed data to maximize agreement with the model

We including coefficients for $t$ and $t^2$ the dependence of the admission rate on time can have some curvature and does not need to be linear.

This can be fitter as a linear model where the inputs are $x_i$, $t_i$ and $t_i^2$

First we select the training period:

1
2
3
4

start_train="2021-07-15"
end_train="2021-08-15"
train_data=data[(data['date'] >= start_train) & (data['date'] <= end_train)].copy()
date0=train_data["date"].min()

Then a 7-day testing period:

1
2
3

test_period=7 # days
test_end=pd.Timestamp(end_train)+pd.DateOffset(days=test_period)
test_data=data[(data["date"]>end_train) & (data["date"]<=test_end)].copy()

1
2
3
4

X_train,P_train=covid.define_variables(train_data,date0)
Y_train=train_data[event]
X_test,P_test=covid.define_variables(test_data,date0)
Y_test=test_data[event]

Fit the model to the train data,

1
2
3

mod=Poisson(Y_train,X_train,exposure=P_train.values)
res=mod.fit()
res.summary()

Optimization terminated successfully.
         Current function value: 32.115992
         Iterations 6

1
2
3

mod=sm.GLM(Y_train,X_train,exposure=P_train,family=sm.families.Poisson(sm.families.links.log()))
glm_res=mod.fit()
glm_res.summary()

Poisson regression is fit by maximum likelihood, there are several choices of Pseudo R-square
Here, what statsmodels implements for poisson regression is McFadden $R_{McF}^2$ that is defined as ratio of log likelihood for the fitted model and log likelihood of a model fitted to just a constant.

1
2

# Pseudo R-squ McF
res.prsquared

0.4492953669633617

1

res.llf, res.llnull

(-52413.29834953645, -95174.97258108124)

1

1-res.llf/res.llnull

0.4492953669633617

Certainly, we can implement other measures of the Pseudo R-square ourselves. Here I did $R_{L}^2$, which is deviance-based. It’s the most analogous index to the squared multiple correlations in linear regression. We will stick to it in our analysis from now on.

1
2

def poisson_deviance(y,y_pred):
    return 2*np.sum(y*np.log(np.maximum(y,1e-12)/y_pred)-(y-y_pred))

1
2
3
4
5
6
7
8

mu_train=Y_train.sum()/P_train.sum()
Y_pred=glm_res.predict(X_train,exposure=P_train)

deviance=poisson_deviance(Y_train,Y_pred)
null_deviance=poisson_deviance(Y_train,mu_train*P_train)
R2=1-deviance/null_deviance
# Pseudo R-squ
R2

0.4725915412837366

If you don’t want bother yourself with the mathematical formulas, here’s a simpler way (my favorite as the laziest person ever). Just fit the model with a constant, and take the log-likelihood/deviance as your null case.

1
2
3
4

mod=sm.GLM(Y_train,X_train['const'],exposure=P_train,family=sm.families.Poisson(sm.families.links.log()))
res_null=mod.fit()
# deviance for a const fit model
res_null.deviance

180966.73552551522

1
2

# null deviance from our math formula
null_deviance

180966.7355255152

After all, the R-squared we have for now is 0.473.

Next let’s check on how the model predicts.

Extrapolation of Time trends

Remember what we saw last time with the linear model? Its projection was a straight line. But the poisson model is telling us a different story, that the admissions will flatten out assuming the current vaccination rate. Well, it is onto something, isn’t it?

1

covid.plot_time_extrapolation(train_data,glm_res,event,date0)

In-sample model Fit

To better evaluate the model fit, we visualize the real observations and model predictions. Here’s some snapshots of the fit for a few selected dates on the training period.

As we can see, curves from the poisson regression model fits the data much better than the linear model, indicating the correlation between admissions and vaccinations are not linear.

1

covid.plot_time_facets(train_data,glm_res,event,date0)

But is this enough? Certainly not, let’s verify the deviance residuals by visualizing them,

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

def plot_deviance_residuals(data,event,res):
    P=data["population"]
    # Create two subplots and unpack the output array immediately
    f, (ax1, ax2) = plt.subplots(1, 2,sharex=True,figsize=(16,4))
    norm=matplotlib.colors.LogNorm(vmin=100_000, vmax=P.max(), clip=False)
    m=ax1.scatter(data["vaccinated"],data[event]/data["population"]*100_000,
                c=P,norm=norm,cmap="Blues",label=event)
    plt.colorbar(m,label="State Population")
    x=np.linspace(0.3,0.75,201)
    x= pd.DataFrame({'vaccinated':x, 'T':31, 'T2':31*31})
    x["const"]=np.ones(201)
    y_pred=res.predict(x)
    ax1.plot(x.vaccinated, y_pred*100_000,"k--",label="predicted") 
    
    idx = data.index
    ax2.scatter(data["vaccinated"], res.resid_deviance[idx],
                c=P,norm=norm,cmap="Blues",label=event)
    ax2.set_title("Residuals")
    ax2.set_xlabel("vaccination")
    ax2.set_ylabel("residual")

    # select outlier state
    outlier_idx=np.argpartition(np.array(res.resid_deviance[idx]), -2)[-2:]
    
    for i in outlier_idx:
        outlier=data.iloc[i]
        ax1.annotate(outlier["state"],(outlier["vaccinated"]+0.01,outlier[event]/outlier["population"]*100_000))
        ax2.annotate(outlier["state"],(outlier["vaccinated"]+0.01,res.resid_deviance[idx].iloc[i]))

The date I choose is Aug.15,

1
2
3

offset=9 
interested_date=data["date"].max()-pd.offsets.Day(offset)
interested_data=data[data["date"]==interested_date]

1

plot_deviance_residuals(interested_data,event,glm_res)

The deviance residuals are defined as the signed square roots of the unit deviances, representing the contributions of individual samples to the deviance. The deviance indicates the extent to which the likelihood of the saturated model exceeds the likelihood of the proposed model. If the proposed model has a good fit, the deviance will be small. If the proposed model has a bad fit, the deviance will be high.

Thus, the deviance residuals are analogous to the conventional residuals: when they are squared, we obtain the sum of squares that we use for assessing the fit of the model. However, while the sum of squares is the residual sum of squares for linear models, for GLMs, this is the deviance.

1
2
3

deviance=poisson_deviance(Y_train,Y_pred)
sum_of_squared_deviance_residuals = sum(glm_res.resid_deviance**2)
deviance, sum_of_squared_deviance_residuals

(95443.38706242564, 95443.38706242583)

In a properly specified model, the deviance is approximately chi-square distributed with n-k-1 degrees of freedom, and the deviance residuals would be independent, standard normal random variables, i.e., ~ norm(0,1). So we would expect the residuals to be mostly in range of -1 to 1, and rarely fall outside the ± 3 limits.

However, the model residuals are 10 times larger than expected and we say there is overdispersion. If overdispersion is present in a dataset, the estimated standard errors and test statistics, the overall goodness-of-fit will be distorted and adjustments must be made.

A more appropriate model will be a Negative Binomial regression.

Model Prediction

We will now use the model without recalibrating to make predictions about admission rates on new data.

This is out of sample evaluation. It is the only way to make sure the model really works.

Out of Sample $R^2$

1
2
3
4
5
6
7

test_admission_mean=Y_test.sum()/P_test.sum()
Y_pred=glm_res.predict(X_test,exposure=P_test)

null_deviance_test=poisson_deviance(Y_test,P_test*test_admission_mean)
deviance_test=poisson_deviance(Y_test,Y_pred)
R2 = 1-deviance_test/null_deviance_test
R2

0.371917369051069

Out of sample $R^2$ is worse than for the in sample (approx. 0.473). I think we’ll all agree, predicting the future is hard.

Negative Binomial Regression

Model Definition

We have data from multiple dates, so it makes sense to incorporate Time as a extra variable.

\begin{eqnarray}
&\eta_i &=c + \beta_1 x_i + \beta_2 t_i + \beta_3 t_i^2 \newline
\mathbb{E}(y_i|x_i) =& \lambda_i P_i &= e^{\eta_i } P_i \newline
&p(y_i|x_i) & = \text{NB}(y_i; \lambda_i P_i, \alpha)
\end{eqnarray}
where

• Each observation $i$ represents a single state at one particular date.
• $y_i$ is the number of hospital admissions on that state.
• $x_i$ is the percentage of the state population fully vaccinated that state.
• $t_i$ is the number of days elapsed since the beginning of the training period.
• $P_i$ is the population of the state.
• the parameters $c$ and $\beta_1,\beta_2,\beta_3$ will be fitted to the observed data to maximize agreement with the model

We including coefficients for $t$ and $t^2$ the dependence of the admission rate on time can have some curvature and does not need to be linear.

This can be fitter as a linear model where the inputs are $x_i$, $t_i$ and $t_i^2$

Variance Comparison to Poisson Model

Introducing a free additional parameter $\alpha$ give more accurate models than simple parametric models like the Poisson distribution by allowing the mean and variance to be different, unlike the Poisson. The negative binomial distribution has a variance $\hat{y}+\hat{y}^2\alpha$ . This can make the distribution a useful overdispersed alternative to the Poisson distribution.

As we can see, with an extra parameter $\alpha$ to control the variance, the expected variance of observations is larger with the Negative Binomial Model than with the Poisson model.

1
2
3
4
5
6
7

alpha = 0.5
y_hat=np.linspace(0,10,201)
plt.plot(y_hat,y_hat,"k--",label="Poisson")
plt.plot(y_hat,y_hat+alpha*y_hat**2,"k-",linewidth=3,label="Negative Binomial")
plt.xlabel(r"$\hat{y}$")
plt.ylabel(r"Var($y|\hat{y}$)")
plt.legend()

Choosing $\alpha$

Negative Binomial is a GLM model with an extra parameter $\alpha$ that we must choose by maximizing the log likelihood.

The train/test data we fit to the Negative Binomial Regression is the same as what we create at start.

1
2
3
4
5
6
7
8
9
10
11
12

lls=[]
alphas=np.linspace(0.1,0.5,500)
for alpha in alphas:
    glm_model=sm.GLM(Y_train,X_train,exposure=P_train,family=sm.families.NegativeBinomial(sm.families.links.log(),alpha=alpha))
    glm_res=glm_model.fit()
    ll=glm_res.llf
    lls.append(ll)
plt.semilogx(alphas,lls) 
plt.xlabel(r"$\alpha$")
plt.ylabel("Log Likelihood")
alpha=alphas[np.argmax(lls)]
print("alpha",alpha)

alpha 0.3004008016032064

The best alpha is 0.3004 given a run of 500 times. Or easier, we can go for the Negative Regression model from statsmodels, which can do the dirty work, optimizing alpha for us.

1
2
3
4

#mod=sm.GLM(Y_train,X_train,exposure=P_train,family=nb)
mod=NegativeBinomial(Y_train,X_train,exposure=P_train)
res=mod.fit(method="lbfgs")
res.summary()

/opt/anaconda3/lib/python3.7/site-packages/statsmodels/discrete/discrete_model.py:2642: RuntimeWarning: divide by zero encountered in log
  llf = coeff + size*np.log(prob) + endog*np.log(1-prob)
/opt/anaconda3/lib/python3.7/site-packages/statsmodels/discrete/discrete_model.py:2642: RuntimeWarning: invalid value encountered in multiply
  llf = coeff + size*np.log(prob) + endog*np.log(1-prob)
/opt/anaconda3/lib/python3.7/site-packages/statsmodels/base/model.py:568: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  ConvergenceWarning)

The alpha values are close,

1
2
3

p=res.params
nb=sm.families.NegativeBinomial(alpha=p["alpha"])
nb.alpha

0.3024936677385029

Compared to Poisson:
- The regression coefficients are very similar.
- The confidence intervals are much wider.
- The t statistics are smaller because we assume a larger variance of residuals.

In summary, switching from Poisson to Negative Binominal yields stable coefficient estimates, but a higher chance for the null hepothesis to be rejected and more reliable confidence intervals.

The statsmodel.family.NegativeBinomial object knows how to compute deviance, and we use it to calculate the R squared,

1
2
3
4
5
6

mu_train=Y_train.sum()/P_train.sum()
Y_pred=res.predict(X_train,exposure=P_train)
deviance=nb.deviance(Y_train,Y_pred)
null_deviance=nb.deviance(Y_train,mu_train*P_train)
R2=1-deviance/null_deviance
R2

0.5937123540474427

In-sample model Fit

1

covid.plot_time_facets(train_data,res,event,date0)

Then most excitingly, let’s check out the deviance residuals again. The residuals are now well behaved, randomly fall into the range of -1 to 1. Though FL, KY and a few other states seem to be outliers. But no worries, we will address them in our next blog, using a more advanced method - mixed models.

1

plot_deviance_residuals(interested_data,event,glm_res)

Model Prediction

Finally, always take a look at the out-of-sample fit.

Though the R squared is not as good as the in-sample, but there’s still an improvement comparing to Poisson regression and the linear regressions.

1
2
3
4
5
6
7

test_admission_mean=Y_test.sum()/P_test.sum()
Y_pred=res.predict(X_test,exposure=P_test)

null_deviance_test=nb.deviance(Y_test,P_test*test_admission_mean)
deviance_test=nb.deviance(Y_test,Y_pred)
R2 = 1-deviance_test/null_deviance_test
R2

0.4595499732848517

Lately I’ve been working on a project predicting credit losses for mortgages with a bunch of loan-level attributes. To be honest, the data wasn’t rich. In pursuit of good predictive power, I had to try from the simple multi-period regression, to generalized linear model with distributions beyond the Gaussian family, then mixed effect models, and even used some curve fitting techniques like splines at last to fix the oversimple (or wrong) assumptions I made at the beginning. Overall, I think the journey was inspiring as it basically shows that one can almost crack any problem with regression analysis, more importantly, in an elegant way.

When it comes to predictive modeling, people always turn to linear regressions, which isn’t wrong, just so you know that linear regression models and OLS make a number of assumptions about the predictor variables, the response variable, and their relationship. Violating these assumptions can result in biased predictions or imprecise coefficients, in short cause your model to be less predictive. Well, the good news is that there are numerous extensions have been developed that allow each of these assumptions to be relaxed, and in some cases eliminated entirely, but only if you know what they are and when to use them. More commonly, I see people get stuck and choose to live with a just-fine-fit model without knowing that they can easily fix it by employing a more generalized regression models. Now let’s see how powerful regressions can get through a trendy modeling case - project covid admissions based on the vaccinate progress.

Preliminaries

1
2
3
4
5
6
7
8

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.colors
import seaborn as sns
import statsmodels.api as sm
import sys
import covid_analysis as covid # some functions developed for this blog

The data I use includes hospital admissions and vaccinations from the CDC website, and the most up-to-date populations by state from the US Census Bureau.

1

data_dir="../data"

1
2
3
4
5
6
7
8
9
10
11

hospitalizations=pd.read_csv(f"{data_dir}/covid_hospitalizations.csv",parse_dates=["date"])
vaccinations=pd.read_csv(f"{data_dir}/covid_vaccinations.csv",parse_dates=["date"])
population=pd.read_csv(f"{data_dir}/population.csv")

event="admissions"
data=hospitalizations

data=data.merge(vaccinations,on=["date","state"])
data=data.merge(population,on="state")
data["vaccinated"]=data["vaccinated"]/data["population"]
data.head()

Here specifically, I am using weighted least squares rather than ordinary least squares to incorporate the knowledge of state populations. WLS is a generalization of OLS and a specialization of generalized least squares, which can relax the assumptions of constant variance, allowing the variances of the observations to be unequal i.e., heteroscedasticity (btw one of my favorite words in statistics, the pronunciation’s fun).

Single Period Regression

Model Definition

The weighted linear regression model is defined as:
\begin{eqnarray}
\hat{y}_i &= c + \beta x_i \newline
y_i &=\hat{y}_i+ \epsilon_i \newline
\end{eqnarray}
and we minimize the weighted least squares error :
\begin{equation}
E_w = \sum_i P_i(y - \hat{y}_i)^2
\end{equation}
where

• Each observation $i$ represents a single state.
• $\hat{y}_i$ is the predicted probability of hospital admissions per persons on that state.
• $y_i$ is the actual number of events per person observed.
• $x_i$ is the percentage of the state population fully vaccinated.
• $P_i$ is the population of the state.
• $\epsilon_i$ is a Gaussian uncorrelated noise with constant variance $\sigma$
• the parameters $c$ and $\beta$ will be fitted to the observed data to maximize agreement with the model

We first fit the regression model to single day of state hospitalization data.

1

offset=9 # if we do not want to look at the last day, subtract here.

Let’s take Aug.15,

1
2

last_date=data["date"].max()-pd.offsets.Day(offset)
last_data=data[data["date"]==last_date]

1
2
3

X=sm.add_constant(last_data["vaccinated"])
P=last_data["population"]
Y=last_data[event]/last_data["population"]

Use statmodels for Weighted Linear Regression

Let’s first use the Weighted Linear Regression model from statmodels.
The exact same estimation can be done by a different way, which I will show later.

1
2
3

mod=sm.WLS(Y,X,P)
res=mod.fit()
res.params

const         0.000127
vaccinated   -0.000190
dtype: float64

So here’s our estimation of the coefficients and some statistics that are helpful for diagnosis. It’s easy to read that the R-squared is 0.327. But we can also calculate it ourselves.

The R-square is defined as ratio of square loss to target variance,

1
2
3

def square_error(y,y_hat,P):
    err=y-y_hat
    return np.sum(P*err**2)

1
2
3
4

Y_pred=res.predict(X)
Y_bar=(Y*P).sum()/P.sum()
R2=1 - square_error(Y,Y_pred,P)/square_error(Y,Y_bar,P)
R2

0.3266337795171833

Use statmodels GLM for Weighted Linear Regression

Exactly the same model can be written in the more general language of Generalized Linear models as

Model Definition

The linear regression model is defined as:
\begin{eqnarray}
&\eta_i &=c + \beta*x_i \newline
\mathbb{E}(y_i|x_i) =& \hat{y}_i &= \eta_i \newline
&p(y_i|x_i) & = N(\hat{y}_i,\sigma)
\end{eqnarray}

The fact that $\eta_i=\hat{y}_i$ means we are using the identity link between the Gaussian family and the linear model $\eta$. One should be aware here that linear model is just a specical case of GLM.

1
2
3

glm_model=sm.GLM(Y,X,family=sm.families.Gaussian(sm.families.links.identity()),var_weights=P)
glm_res=glm_model.fit()
glm_res.params

const         0.000127
vaccinated   -0.000190
dtype: float64

Results are the same. But with the GLM language we will be able to perform regression many more kinds of variables.

The generalization of square error for a GLM model is called the deviance. This is the square error the Gaussian Family used in linear regression. Let’s calculate the in-sample R-squared in the GLM way and compare it with the previous.

1

glm_res.deviance, square_error(Y,Y_pred,P)

(0.12246143037040369, 0.12246143037040369)

The resuls also include the deviance for the null model fitted to just a constant

1

glm_res.null_deviance,square_error(Y,Y_bar,P)

(0.1818645287591, 0.1818645287591)

The ratio defines an $R^2$ value. Not surprisingly, the results are exactly the same. Overall, I admit a 32.7% R-squared model can barely be said to be a good model. But no worries, all can be fixed.

1

1-glm_res.deviance/glm_res.null_deviance

0.3266337795171833

Now, as we have a better understanding of the data and the modeling methodologies, let’s move on to model the time series.

The data we have here is time series data, indeed panel data, as time series and cross-sectional data are both special cases of panel data that are in one dimension only. But we will deal with the multi-subjects matter in the next blog(probably should have a spoiler tag oops). Anyways, let’s first start with the multi-period regression, here I have trend variables added to the model to represent time.

Multiple Period Regression

Model Definition

\begin{eqnarray}
&\eta_i &=c + \beta_1 x_i + \beta_2 t_i + \beta_3 t_i^2 \newline
\mathbb{E}(y_i|x_i) =& \hat{y}_i &= \eta_i \newline
&p(y_i|x_i) & = N(\hat{y_i},\sigma^2)
\end{eqnarray}
where

• Each observation $i$ represents a single state at one particular date.
• $y_i$ is the number of events on that state.
• $\hat{y}_i$ is the predicted number of events on that state.
• $x_i$ is the percentage of the state population fully vaccinated on that state.
• $t_i$ is the number of days elapsed since the beginning of the training period.
• $\sigma^2$ is the level of noise of the observations (the variance of the residuals).
• the parameters $c$ and $\beta_1,\beta_2,\beta_3$ will be fitted to the observed data to maximize agreement with the model

This can be fitter as a linear model where the inputs are $x_i$, $t_i$ and $t_i^2$

First we select the training period:

1
2
3
4

start_train="2021-07-15"
end_train="2021-08-15"
train_data=data[(data['date'] >= start_train) & (data['date'] <= end_train)].copy()
date0=train_data["date"].min()

Then a 7-day testing period:

1
2
3

test_period=7 # days
test_end=pd.Timestamp(end_train)+pd.DateOffset(days=test_period)
test_data=data[(data["date"]>end_train) & (data["date"]<=test_end)].copy()

I write a function to generate the desired matrix,

1
2
3
4
5
6
7
8
9

def define_variables(data,date0):
        N=len(data)
        T=(data["date"]-date0).dt.days
        X=data["vaccinated"].copy().to_frame()
        X["T"]=T
        X["T2"]=T**2
        P=data["population"]
        X["const"]=np.ones(N)
        return X,P

1
2
3
4

X_train,P_train=covid.define_variables(train_data,date0)
Y_train=train_data[event]/train_data["population"]
X_test,P_test=covid.define_variables(test_data,date0)
Y_test=test_data[event]/test_data["population"]

Weighted Multi-period Linear Regression

1
2
3

mod=sm.WLS(Y_train,X_train,P_train)
res=mod.fit()
res.summary()

1
2
3
4
5
6
7

Y_bar=np.sum(Y_train*P_train)/np.sum(P_train) # weighted average
Y_pred=res.predict(X_train)

deviance=square_error(Y_train,Y_pred,P_train)
null_deviance=square_error(Y_train,Y_bar,P_train)
R2=1-deviance/null_deviance
R2

0.3726019173009557

Extrapolation of Time trends

The linear model projects the admission to be increasing steadily if assuming current vaccination rate.

1

covid.plot_time_extrapolation(train_data,res,event,date0)

In-sample model Fit

To better evaluate the model fit, we visualize the real observations and model prodictions. Here’s some snapshots of the fit for a few selected dates on the training period.

Observations are arranged by state and date.

Don’t have to say much, you can tell the weighted linear model doesn’t fit well as the corelation between vaccinations and admissions looks more like a curve than a straight line.

1

covid.plot_time_facets(train_data,res,event,date0)

We can further plot out the residuals. I choose one date from the dates above, Aug.15 and apply the following function to generate its residual plot.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

def plot_residuals(data,event,res):
    P=data["population"]
    # Create two subplots and unpack the output array immediately
    f, (ax1, ax2) = plt.subplots(1, 2,sharex=True,figsize=(16,4))
    norm=matplotlib.colors.LogNorm(vmin=100_000, vmax=P.max(), clip=False)
    m=ax1.scatter(data["vaccinated"],data[event]/data["population"]*100_000,
                c=P,norm=norm,cmap="Blues",label=event)
    plt.colorbar(m,label="State Population")
    x=np.linspace(0.3,0.75,201)
    x= pd.DataFrame({'vaccinated':x, 'T':31, 'T2':31*31})
    x["const"]=np.ones(201)
    y_pred=res.predict(x)
    ax1.plot(x.vaccinated, y_pred*100_000,"k--",label="predicted") 
    
    idx = data.index
    ax2.scatter(data["vaccinated"], res.resid[idx],
                c=P,norm=norm,cmap="Blues",label=event)
    ax2.set_ybound(lower=-0.00008,upper=0.00008)
    ax2.set_title("Residuals")
    ax2.set_xlabel("vaccination")
    ax2.set_ylabel("residual")

    # select outlier state
    outlier_idx=np.argpartition(np.array(res.resid[idx]), -2)[-2:]
    
    for i in outlier_idx:
        outlier=data.iloc[i]
        ax1.annotate(outlier["state"],(outlier["vaccinated"]+0.01,outlier[event]/outlier["population"]*100_000))
        ax2.annotate(outlier["state"],(outlier["vaccinated"]+0.01,res.resid[idx].iloc[i]))

We can see that

• Predicted events can be negative for large vaccination rates, saying we are over estimating how fast the admission rate is increasing.
• residual variance is negatively correlated with vaccinations rate i.e., smaller vaccination rates have larger magnitude residuals. A Poisson model can fix those issues, which will be the topic of our next blog.
• Some states like FL and KY look like outliers, indicating that they may have different intercepts than the others, which is not surprising as states are different subjects and each may have unique fundamentals. A mixed model would be a good solution for this.

1
2

plot_residuals(last_data,event,res)

Model Prediction

We will now use the model without recalibrating to make predictions about admission rates on new data.

This is out of sample evaluation. It is the only way to make sure the model really works.

Out of Sample Model Fit

Besides greater residuals, the out of sample evaluation is basically telling us a same story.

1

covid.plot_time_facets(test_data,res,event,date0)

Day of the Week Effects

Before we move on to GLM, there’s another trick to play with. The data reporting obviously has day of the week effects,

1
2
3

aggregate=data.groupby(["date"])[event].sum()
aggregate=aggregate.to_frame().reset_index()
aggregate.plot(x="date",y=event)

We use a new function to generate the train data, more dummy variables are added to represent days of a week,

1
2
3
4
5
6
7
8
9
10
11

def define_variables2(data,date0):
        N=len(data)
        T=(data["date"]-date0).dt.days
        X=data["vaccinated"].copy().to_frame()
        X["T"]=T
        X["T2"]=T**2
        X["const"]=np.ones(N)
        Z=pd.get_dummies(data["date"].dt.day_name(),drop_first=True)
        X=pd.concat([X,Z],axis=1)
        P=data["population"]
        return X,P

1
2

X_train,P_train=define_variables2(train_data,date0)
X_train.head()

1
2
3

mod=sm.WLS(Y_train,X_train,P_train)
res=mod.fit()
res.params

vaccinated   -1.510900e-04
T             1.028112e-06
T2           -6.639196e-10
const         8.186223e-05
Monday       -2.235516e-06
Saturday     -8.318830e-07
Sunday       -3.189080e-06
Thursday     -2.233662e-08
Tuesday       3.232847e-07
Wednesday    -2.672040e-07
dtype: float64

1

res.rsquared

0.37572844250241666

$R^2$ is 0.376, higher than the 0.373 we previously have. Day of the week provides a small improvement on $R^2$ compared to original model.