Statin use is associated with lower risk of dementia in stroke patients: a community-based cohort study with inverse probability weighted marginal structural model analysis

Current evidence is inconclusive on cognitive benefits or harms of statins among stroke patients, who have high risk of dementia. This observational cohort study investigated the association between statin use and post-stroke dementia using data from the Clinical Practice Research Datalink. Patients without prior dementia who had an incident stroke but received no statins in the preceding year were followed for up to 10 years. We used inverse probability weighted marginal structural models to estimate observational analogues of intention-to-treat (ITT, statin initiation vs. no initiation) and per-protocol (PP, sustained statin use vs. no use) effects on the risk of dementia. To explore potential impact of unmeasured confounding, we examined the risks of coronary heart disease (CHD, positive control outcome), fracture and peptic ulcer (negative control outcomes). In 18,577 statin initiators and 14,613 non-initiators (mean follow-up of 4.2 years), the adjusted hazard ratio (aHR) for dementia was 0.70 (95% confidence interval [CI] 0.64–0.75) in ITT analysis and 0.55 (95% CI 0.50–0.62) in PP analysis. The corresponding aHRITT and aHRPP were 0.87 (95% CI 0.79–0.95) and 0.70 (95% CI 0.62–0.80) for CHD, 1.03 (95% CI 0.82–1.29) and 1.09 (95% CI 0.77–1.54) for peptic ulcer, and 0.88 (95% CI 0.80–0.96) and 0.86 (95% CI 0.75–0.98) for fracture. Statin initiation after stroke was associated with lower risk of dementia, with a potentially greater benefit in patients who persisted with statins over time. The observed association of statin use with post-stroke dementia may in part be overestimated due to unmeasured confounding shared with the association between statin use and fracture. Supplementary Information The online version contains supplementary material available at 10.1007/s10654-022-00856-7.


ICD-10 Codes for stroke used in HES and ONS ICD Code Description
Fracture of malar and maxillary bones S02.6 Fracture of mandible S02.7 Multiple fractures involving skull and facial bones S02.8 Fractures of other skull and facial bones S02.9 Fracture of skull and facial bones, part unspecified S12 Fracture of neck S12.0 Fracture of first cervical vertebra S12. 1 Fracture of second cervical vertebra S12. 2 Fracture of other specified cervical vertebra S12.7 Multiple fractures of cervical spine S12.8 Fracture of other parts of neck S12.9 Fracture of neck, part unspecified S22 Fracture of rib(s), sternum and thoracic spine S22.0 Fracture of thoracic vertebra S22.

Appendix 2. Potential confounders
Depending on the time to measure the confounder covariates, three sets of covariates measured using the CPRD and HES data were considered in our analysis for adjustment: baseline characteristics, characteristics by the 90 th day of stroke, and time-varying characteristics.

Baseline characteristics
Baseline covariates included demographics, comorbidities, lifestyle, healthcare utilisation and co-medications prior to or on the date of index stroke (Figure 1), unless otherwise specified below. These covariates were used to adjust for between-group imbalance in baseline characteristics and 90-day post-stroke loss to follow-up (i.e. baseline selection).
Demographic covariates included age at index stroke, gender and socioeconomic status.
Index of Multiple Deprivation (IMD) categorised into quintiles was used as a marker for socioeconomic status. It includes seven domains: income; employment; health and disability; education, skills and training; barriers to housing and services; crime; and living environment. When the patient-level IMD was missing, we replaced it with the practice-level IMD. Stroke subtype, year of stroke diagnosis, and healthcare utilisation (defined as the number of primary care consultation within one year before index stroke) were also included as baseline covariates.
We used previously developed algorithms to identify the following baseline comorbidities, which are considered to be associated with dementia: atrial fibrillation, alcohol use disorders, coronary heart disease, diabetes, heart failure, hyperlipidaemia, hypertension, peripheral artery disease, transient ischemic attack, anxiety, asthma, chronic obstructive pulmonary disease, depression, epilepsy, Parkinson's disease and rheumatoid arthritis. Lifestyle variables, including smoking (current, former, and never smoking) and body mass index (a continuous variable), were defined using the most recent available values before the index stroke.
Co-medications included anticoagulant, antiplatelet, antidiabetic, antihypertensive and nonstatin lipid-lowering treatments, which may be associated with the risk of dementia. We defined these medications via product codes recorded within the 365 days prior to the index stroke.

A48
These post-stroke covariates included comorbidities, lifestyle and co-medications. We identified these covariates using the same approaches as mentioned above but extended the time window of ascertainment up to the 90 th day after stroke. These covariates, as well as baseline demographics, stroke subtype, year of stroke diagnosis and healthcare utilisation, were used to adjust for censoring associated to characteristics at the start of the time-at-risk period.

Time-varying characteristics
Time-varying covariates included comorbidities, lifestyle and co-medications. We identified these using the same approaches as mentioned above in 30-day intervals of follow-up since the index stroke ( Figure 1). For covariates not measured in a particular 30-day interval, we carried the last observation forward to impute the value. These covariates, as well as baseline demographics, stroke subtype, year of stroke diagnosis, and healthcare utilisation, were used to adjust for censoring associated with time-varying confounding during the first 90 days (baseline selection) and the subsequent time-at-risk period.

Appendix 4. Details on effect estimation using inverse probability weighted marginal structural models
We estimated the treatment effects via a two-stage process. Firstly, we estimated the four stabilised weights using logistic regression. In the second stage, we used weighted pooled logistic regression models to estimate the observational analogues of intention-to-treat (ITT) and per-protocol (PP) effects.

First stage: weight estimation
Four types of weights were used: baseline treatment weights, baseline selection weights, follow-up weights and treatment persistence weightsto adjust for baseline confounding, baseline selection bias, loss to follow-up, and non-persistence with the initial treatment, respectively.
To reduce the variability of the original weights, we estimated the stabilised version for each weight. We estimated the stabilised weights by fitting two logistic regression models for denominator and numerator, respectively. Covariates included in denominator and numerator models for each type of weight are summarised in  The models of weights were fitted separately by treatment group, except for baseline treatment weight.

Second stage: ITT and PP effect estimation
In this stage, we used weighted pooled logistic regression models to estimate two types of

Mathematical details
Details about the models used in the two stages are listed as below.

Baseline treatment weight
Inverse probability weighting derived from propensity scores was used to adjust for the imbalance in the baseline characteristics between the two treatment groups (i.e., propensity score weighting).
Let In the denominator model, the propensity score, denoted as PSd, was estimated using a logistic regression model in which the treatment indicator (A = 1 or 0) was regressed on the all the observed baseline covariates listed above. In the numerator model, we calculated the probability of a subject initiating statin treatment conditional only on baseline age (as restricted cubic smoothing splines with four knots) and gender, denoted as PSn, using the same modelling approach to stabilise the weights.
The stabilised inverse probability of treatment weights was calculated as: The baseline treatment weight was the same for each 30-day interval of a subject through the time-at-risk period.

Baseline selection weight
Potential selection bias due to death or transfer-out within the first 90 days following index stroke before the time-at-risk period was accounted for via the use of inverse probability of censoring weights.
Let C(t) be a dichotomous variable taking the value 1 if a subject is censored in time interval t and 0 otherwise (for this type of weight calculation, t=1, 2, 3). ̅ ( ) denotes censoring history (that is, the vector of values of C(k) from k = 1 to k = t − 1). ̅ ( − 1) denotes time-varying covariate history (that is, the matrix of values of L(k-1) from k = 1 to k = t).To adjust for baseline selection bias in the marginal structural model, we derived the probability of remaining uncensored up to time t, and calculated the stabilised version of the baseline selection weight as: In this equation (2), A denotes initial treatment, which is a dichotomous variable, taking the value 1 if a subject received a statin within the first 90 days of stroke and 0 otherwise (i.e., the weight was estimated separately by treatment group). We treated it as a baseline variable, instead of a time-varying variable over this period, because the research question of interest is "starting a stain within the first 3 months of stroke" vs "never starting a statin in the first 3 months". V denotes a subset of the baseline covariates and time, which included age (as A53 restricted cubic smoothing splines with four knots), gender, and 30-day interval (as indicator variables) in this model. L(k-1) denotes the values of all time-varying covariates at k=t. As an exception, L(0) was defined to be the complement baseline covariate set of V.

Follow-up weight
Potential bias due to death or transfer-out during the time-at-risk period was accounted for via the use of inverse probability of follow-up weights. Using similar approach mentioned above, we derived the probability of remaining uncensored up to time interval t, and calculated the stabilised version of the follow-up weight as: (3) The differences between (2) and (3) were as below. The weight was calculated from the fourth 30-day interval, i.e., t=4, 5,…, T). V included age (as restricted cubic smoothing splines with four knots), gender, and time (as restricted cubic smoothing splines with four knots). Specifically, L(3) denotes the values of all time-varying covariates measured up to the third 30-day interval following stroke. The weight was estimated separately by treatment group.

Treatment persistence weight
In the observational analogue of PP analysis, we artificially censored patients discontinuing their initial treatments. To account for the treatment non-persistence during the time-at-risk period, we derived the probability of remaining on initial treatment up to time interval t, and calculated the stabilised version of the treatment persistence weight as: (4) In the equation (4), CA(t) was defined as a dichotomous variable taking the value 1 if a subject changed the initial treatment in time interval t and 0 otherwise (t=4, 5,…, t). ̅ ( ) denotes persistence history (that is, the vector of values of CA(k) from k = 4 to k = t − 1). A denotes initial treatment, which is a dichotomous variable, taking the value 1 if a subject received a statin within the first 90 days of stroke and 0 otherwise (i.e., a=0, 1). All the patients were considered to persist with the initial treatment over the time-at-risk period until they discontinued statins (statin initiation group) or start on statins (statin non-initiation group) (as defined in the Follow-up section of the main paper) and thus were artificially censored. ̅ ( − 1) denotes history of time-varying confounders for a subject (that is, the matrix of values of L(k-1) from k = 4 to k = t). In this model, V included age (as restricted A54 cubic smoothing splines with four knots), gender, and time (as restricted cubic smoothing splines with four knots). L(k-1) denotes the values of all time-varying covariates at k=t. As an exception, L(3) was defined to be the complement covariate set of V measured up to the third 30-day interval following stroke. The weight was estimated separately by treatment group.

Estimating the outcome model
Using the weights calculated as above, we estimated the hazard ratio for the effects of interest using a weighted pooled logistic regression model, which approximates a time-dependent Cox model [1].
We defined the outcome event at each time point over the time-at-risk period, = 4, 5, … , , to be ( ) . ( ) = 0 persisted for all 30-day intervals until the interval in which the outcome event occurs, at which point ( ) = 1. The risk of the outcome was evaluated with respect to the initial treatment, i.e., ( ) depends on , which was assumed unchanged over the whole time-at-risk period in the ITT analysis. For this analysis, we fitted three different outcome models using the pooled logistic regression as: In model 1, we adjusted for the between-group imbalance in baseline characteristics. The inverse probability of weights used in the model was .
In model 2, we additionally adjusted for potential selection bias within the first 90 days of stroke. The inverse probability of weights used in the model was × ( ).
In model 3, we additionally adjusted for loss to follow-up during the time-at-risk period after the first 90 days of stroke. The inverse probability of weights used in the model was For the PP analysis in which patients deviated from their initial treatment were artificially censored, we fitted model 1 and 2 using the same approach above. In model 3, we additionally adjusted for treatment non-persistence during the time-at-risk period after the first 90 days of stroke as:
In (5) to (8), was the parameter of interest. Because the stabilisation of the weight means that treatment initiation is balanced conditional on all variables (age, gender and time), these variables included in the numerator model of the stabilised weight were included in the outcome model. The marginal structural models also accounted for the dependence between observations from the same subject which were introduced by the weighting process, and therefore the variance was estimated using a robust variance estimator [2].    The incidence rates of dementia in stroke patients observed in our study were 1.2 to 2.8 times as high as that in the general population from the Cognitive Function and Ageing Studies II [1], approximately in line with the risk ratios summarised in the meta-analysis, which ranged from 1.42 (95% CI 1.20-1.67) to 3.28 (1.92-5.62) [2].  For PP analysis, artificial censoring due to deviation from initial treatment during time-at-risk period was additionally accounted for. *Extreme values were truncated to 1 st and 99 th percentiles. Abbreviations: CHD, coronary heart disease; ITT, intention-to-treat; max, maximum weight; min, minimum weight; PP, per-protocol; PSD, post-stroke dementia; SD, standard deviation.  weight calculation), 29,601 patients who survived at month 3 were included in the outcome models. CHD: Of 28,946 eligible patients (contributing to weight calculation), 26,098 patients who survived at month 3 were included in the outcome models.

Appendix 8. Sensitivity analysis
Fracture: Of 22,850 eligible patients (contributing to weight calculation), 20,545 patients who survived at month 3 were included in the outcome models. Peptic ulcer: Of 31,042 eligible patients (contributing to weight calculation), 27,647 patients who survived at month 3 were included in the outcome models. All the outcome models adjusted for baseline characteristics, baseline selection, loss to follow-up and months of follow-up. For PP analysis, artificial censoring due to deviation from initial treatment during time-at-risk period was additionally accounted for. Abbreviations: aHR, adjusted hazard ratio; CHD, coronary heart disease; CI, confidence interval; IMD, Index of Multiple Deprivation; ITT, intention-to-treat; PP, per-protocol.  ,700 patients who survived at month 3 were included in the outcome models. **38,034 eligible patients contributed to weight calculation, of whom 33,700 patients who survived at month 3 were included in the outcome models. #Deviation from initial treatment during follow-up was artificially censored in all the models. a. Model 1: adjusted for baseline characteristics and months of follow-up. b. Model 2: adjusted for baseline selection plus Model 1. c. Model 3: adjusted for loss to follow-up plus Model 2. For PP analysis, artificial censoring due to deviation from initial treatment during time-at-risk period was additionally accounted for. Abbreviations: BMI, body mass index; CI, confidence interval; HR, hazard ratio; ITT, intention-to-treat; PP, per-protocol. ,601 patients who survived at month 3 were included in the outcome models. **33,190 eligible patients contributed to weight calculation, of whom 29,601 patients who survived at month 3 were included in the outcome models. #Deviation from initial treatment during follow-up was artificially censored in all the models. a. Model 1: adjusted for baseline characteristics and months of follow-up. b. Model 2: adjusted for baseline selection plus Model 1. c. Model 3: adjusted for loss to follow-up plus Model 2. For PP analysis, artificial censoring due to deviation from initial treatment during time-at-risk period was additionally accounted for. Abbreviations: CI, confidence interval; HR, hazard ratio; ITT, intention-to-treat; PP, per-protocol.  ,317 patients who survived at month 3 were included in the outcome models. **32,906 eligible patients contributed to weight calculation. Of whom, 29,317 patients who survived at month 3 were included in the outcome model. #Deviation from initial treatment during follow-up was artificially censored in all the models. a. Model 1: adjusted for baseline characteristics and months of follow-up. b. Model 2: adjusted for baseline selection plus Model 1. c. Model 3: adjusted for loss to follow-up plus Model 2. For PP analysis, artificial censoring due to deviation from initial treatment during time-at-risk period was additionally accounted for. Abbreviations: CI, confidence interval; HR, hazard ratio; ITT, intention-to-treat; PP, per-protocol. .59 (0.52-0.68) *Of 32,520 eligible patients, 28,931 patients who survived at month 3 were included in the outcome models. **32,520 eligible patients contributed to weight calculation. Of whom, 28,931 patients who survived at month 3 were included in the outcome model. #Deviation from initial treatment during follow-up was artificially censored in all the models. a. Model 1: adjusted for baseline characteristics and months of follow-up. b. Model 2: adjusted for baseline selection plus Model 1. c. Model 3: adjusted for loss to follow-up plus Model 2. For PP analysis, artificial censoring due to deviation from initial treatment during time-at-risk period was additionally accounted for. Abbreviations: CI, confidence interval; HR, hazard ratio; ITT, intention-to-treat; PP, per-protocol.  were included in the outcome models. **32,906 eligible patients contributed to weight calculation. Of them, 26,205 patients who survived at month 6 were included in the outcome models. #Deviation from initial treatment during follow-up was artificially censored in all the models. a. Model 1: adjusted for baseline characteristics and months of follow-up. b. Model 2: adjusted for baseline selection plus Model 1. c. Model 3: adjusted for loss to follow-up plus Model 2. For PP analysis, artificial censoring due to deviation from initial treatment during time-at-risk period was additionally accounted for. Abbreviations: aHR, adjusted hazard ratio; cHR, crude hazard ratio; CI, confidence interval; ITT, intention-totreat; PP, per-protocol.  ,601 patients who survived at month 3 were included in the outcome models. **33,190 eligible patients contributed to weight calculation, 29,601 patients who survived at month 3 were included in the outcome models. #Deviation from initial treatment during follow-up was artificially censored in all the models. a. Model 1: adjusted for baseline characteristics and months of follow-up. b. Model 2: adjusted for baseline selection plus Model 1. c. Model 3: adjusted for loss to follow-up plus Model 2. For PP analysis, artificial censoring due to deviation from initial treatment during time-at-risk period was additionally accounted for. Abbreviations: CI, confidence interval; HR, hazard ratio.