New SSC package: -sccsdta- for generating the dataset for the self-controlled case series method (SCCS)

[Niels Henrik Bruun]

Thanks to Kit Baum, there is a new package -sccsdta- at the SSC.

The SCCS method examines the association between a time-varying exposure and an event outcome.
The study samples only cases, and it requires that an event has occurred during the observation period.
The method doesn't compare incidences for cases with incidences for references.
Instead, it contrasts incidences in periods of risk with incidences in periods where the case is not at risk.
In this approach, cases serve as their control for fixed confounders.
It's feasible to adjust for time effects such as age.
The intervals are marked by risk, time, and individual.
For each interval, the number of incidences and the width of the interval are determined.

An introduction is: https://www.medicine.mcgill.ca/epide...l/Tutorial.pdf or https://www.bmj.com/content/354/bmj.i4515

Below is a reproduction of example 3.1 in the tutorial from the first reference:

Hospital records indicate an association between MMR vaccination and viral meningitis.
Specifically, using a certain live mumps vaccine, known as the Urabe strain, has been linked to an increased risk of viral meningitis.
Instances of viral meningitis were identified in 10 children during their second year of life.

Event day (day of meningitis) and a day for exposure (day of vaccination) for 10 children are listed in the dataset

* Example generated by -dataex-. For more info, type help dataex
clear
input float(eventday exday)
398 458
413 392
449 429
455 433
472 432
474 395
485 470
524 496
700 428
399 716
end

The observation period was from the 366th to the 730th day of age.
Evidence led to the definition of the risk period as the 15th to the 35th day following the administration of the MMR vaccine.
Age groups were 366 to 547 days and 548 to 730 days.

Using the sccsdta command:

. sccsdta eventday exday, enter(365) exit(730) riskpoints(14 35) timepoints(547)

transforms the dataset into (first 3 children):

. list if _id < 4, sepby(_id) noobs

  +---------------------------------------------------------------------------------------+
  | eventday   exday   _id   _start   _stop   _nevents      _exgr        _tmgr   _inter~l |
  |---------------------------------------------------------------------------------------|
  |      398     458     1      365     472          1         no   ]365; 547]        107 |
  |      398     458     1      472     493          0   ]14; 35]   ]365; 547]         21 |
  |      398     458     1      493     547          0         no   ]365; 547]         54 |
  |      398     458     1      547     730          0         no   ]547; 730]        183 |
  |---------------------------------------------------------------------------------------|
  |      413     392     2      365     406          0         no   ]365; 547]         41 |
  |      413     392     2      406     427          1   ]14; 35]   ]365; 547]         21 |
  |      413     392     2      427     547          0         no   ]365; 547]        120 |
  |      413     392     2      547     730          0         no   ]547; 730]        183 |
  |---------------------------------------------------------------------------------------|
  |      449     429     3      365     443          0         no   ]365; 547]         78 |
  |      449     429     3      443     464          1   ]14; 35]   ]365; 547]         21 |
  |      449     429     3      464     547          0         no   ]365; 547]         83 |
  |      449     429     3      547     730          0         no   ]547; 730]        183 |
  +---------------------------------------------------------------------------------------+

We estimate the incidence rate ratio of the risk period versus the no-risk period using a Poisson regression and looking at the i._exgr estimate:

. poisson _nevents i._exgr i._tmgr i._id, exposure(_interval) irr

We get a better estimate report using the command regmat (ssc install matrixtools - works from version 13.1)

. regmat, o(_nevents) e(i._exgr) a("i._tmgr i._id") eform d(3) label btext(irr) names("") base: poisson, exposure(_interval)
  
-----------------------------------------------------------------------------------
                                  irr  se(irr)  Lower 95% CI  Upper 95% CI  P value
-----------------------------------------------------------------------------------
Events(#)  At risk (no)         1.000                                              
           At risk (]14; 35])  12.037    2.031         3.002        48.259    0.000
-----------------------------------------------------------------------------------

Enjoy

[Niels Henrik Bruun]

Thanks to Kit Baum, the package package -sccsdta- was updated at the SSC.
The help Is extended and rewritten.
I've also revised the options for setting the dates.
Hopefully, it is more clear now, what you get.
Enjoy

[gunnu singh]

Hello Neils,
I used the following code to calculate IRR. I tried to generate similar results using sccsdta but to no avail. Can you please help me implement it using sccsdta?

clear
input indiv str11 vaccination_dt_str str11 hospitalization_dt_str str11 birthdt_str
1 "11feb2004" "18mar2004" "01jan2004"
1 "22nov2004" "28nov2004" "01jan2004"
2 "11feb2004" "22feb2004" "01jan2004"
2 "11feb2004" "12feb2004" "01jan2004"
2 "22nov2004" "28nov2004" "01jan2004"
3 "11feb2004" "25mar2004" "01jan2004"
3 "22nov2004" "26nov2004" "01jan2004"
3 "22nov2004" "30nov2004" "01jan2004"
3 "22nov2004" "21feb2005" "01jan2004"
4 "11feb2004" "25may2004" "01jan2004"
4 "22nov2004" "15dec2004" "01jan2004"
5 "11feb2004" "26jun2004" "01jan2004"
end

* Convert date variables
gen vaccination_dt = date(vaccination_dt, "DMY")
gen hospitalization_dt = date(hospitalization_dt, "DMY")
gen birthdt = date(birthdt, "DMY")

gen admtday = hospitalization_dt - birthdt
//Generate exposure group cut points
//Risk period : 0-28 days, control period: 35-83 days
gen cutp40 = vaccination_dt - birthdt -1 //Include vaccination day
gen cutp41 = cutp40 + 29
gen cutp42 = cutp40 + 35
gen cutp43 = cutp40 + 84

drop birthdt vaccination_dt hospitalization_dt

reshape long cutp,i(indiv  admtday) j(type)

sort indiv admtday cutp type

//Number of hospitalization events within each interval
by indiv : gen int nevents = 1 if admtday > cutp[_n-1]+0.5 & admtday <= cutp[_n]+0.5

// A second ED visit\hospitalization during same risk or control window will be counted as a separate event if it occurred >7 days after the first.;
sort indiv type nevents admtday cutp 

collapse (sum) nevents, by(indiv  cutp type)
//Generate intervals
by indiv : gen int interval = cutp[_n] - cutp[_n-1]

//Exposure groups
gen exgr = type - 40 if type > 40
count if exgr >= .
local nmiss = r(N)
local nchange = 1
while `nchange' > 0 {
    by indiv: replace exgr = exgr[_n+1] if exgr >= . & !missing(exgr[_n-1])
    count if exgr >= .
    local nchange = `nmiss' - r(N)
    local nmiss = r(N)
}
replace exgr = . if type == 40

replace exgr = 0 if exgr == 3 //Control period
replace exgr = . if inlist(exgr ,2) //Washout period

drop if inlist(interval ,.,0) | exgr == .

sort indiv cutp exgr
drop cutp type

generate loginterval = log(interval)

xi: xtpoisson nevents i.exgr, fe i(indiv ) offset(loginterval ) irr

i.exgr            _Iexgr_0-1          (naturally coded; _Iexgr_0 omitted)
note: 1 group (2 obs) dropped because of all zero outcomes

Iteration 0:   log likelihood = -11.203621  
Iteration 1:   log likelihood = -8.1165048  
Iteration 2:   log likelihood = -8.1153088  
Iteration 3:   log likelihood = -8.1153087  

Conditional fixed-effects Poisson regression    Number of obs     =         16
Group variable: indiv                           Number of groups  =          4

                                                Obs per group:
                                                              min =          4
                                                              avg =        4.0
                                                              max =          4

                                                Wald chi2(1)      =       4.91
Log likelihood  = -8.1153087                    Prob > chi2       =     0.0266

------------------------------------------------------------------------------
     nevents |        IRR   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
    _Iexgr_1 |   5.913793   4.741583     2.22   0.027     1.228532    28.46726
 loginterval |          1  (offset)
------------------------------------------------------------------------------

I tried using sccsdta like below:

gen eventday = hospitalization_dt - birthdt
gen expday = vaccination_dt - birthdt -1
gen enday = vaccination_dt - birthdt -1
sccsdta eventday expday, en(enday)  r(0 28 34) t(84) 
. xtpoisson _nevents i._exgr , fe i(_rowid) exposure(_interval) irr
note: 3 groups (9 obs) dropped because of all zero outcomes

Iteration 0:   log likelihood = -8.7278736  
Iteration 1:   log likelihood = -4.9278887  
Iteration 2:   log likelihood = -4.8018012  
Iteration 3:   log likelihood = -4.7742483  
Iteration 4:   log likelihood = -4.7689667  
Iteration 5:   log likelihood = -4.7677265  
Iteration 6:   log likelihood = -4.7674282  
Iteration 7:   log likelihood = -4.7673685  
Iteration 8:   log likelihood = -4.7673588  
Iteration 9:   log likelihood = -4.7673564  
Iteration 10:  log likelihood = -4.7673559  

Conditional fixed-effects Poisson regression    Number of obs     =         27
Group variable: _rowid                          Number of groups  =          9

                                                Obs per group:
                                                              min =          3
                                                              avg =        3.0
                                                              max =          3

                                                Wald chi2(2)      =       5.22
Log likelihood  = -4.7673559                    Prob > chi2       =     0.0734

------------------------------------------------------------------------------
    _nevents |        IRR   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _exgr |
    ]0; 28]  |   6.248765   5.009883     2.29   0.022     1.298231    30.07712
   ]28; 34]  |   4.37e-07   .0013496    -0.00   0.996            0           .
ln(_inter~l) |          1  (exposure)
------------------------------------------------------------------------------

[Niels Henrik Bruun]

Hello Gunnu
I think your problem is that your start and end days are the same.
I did set the end day to 365 in the option timepoints.

cls
. clear

. input indiv str11 vac_dt_str str11 hosp_dt_str str11 birth_dt_str

          indiv   vac_dt_str  hosp_dt_str  birth_dt_~r
  1. 1 "11feb2004" "18mar2004" "01jan2004"
  2. 1 "22nov2004" "28nov2004" "01jan2004"
  3. 2 "11feb2004" "22feb2004" "01jan2004"
  4. 2 "11feb2004" "12feb2004" "01jan2004"
  5. 2 "22nov2004" "28nov2004" "01jan2004"
  6. 3 "11feb2004" "25mar2004" "01jan2004"
  7. 3 "22nov2004" "26nov2004" "01jan2004"
  8. 3 "22nov2004" "30nov2004" "01jan2004"
  9. 3 "22nov2004" "21feb2005" "01jan2004"
 10. 4 "11feb2004" "25may2004" "01jan2004"
 11. 4 "22nov2004" "15dec2004" "01jan2004"
 12. 5 "11feb2004" "26jun2004" "01jan2004"
 13. end

. 
. * Convert date variables
. gen vac_days = date(vac_dt_str, "DMY") - date(birth_dt_str, "DMY") - 1

. gen hosp_days = date(hosp_dt_str, "DMY") - date(birth_dt_str, "DMY")

. drop *_str

. sccsdta hosp_days vac_days, enter(0) riskpoints(0 28 34) t(84 365)

The SCCS regression summary table

                                 |       IRR       [95%        CI]   P(IRR=1) 
---------------------------------+-------------------------------------------
At risk                          |                                           
                            ctrl |     1.000          .          .          . 
                         ]0; 28] |    25.099      5.657    111.363      0.000 
                        ]28; 34] |     0.000      0.000          .      0.996 
---------------------------------+-------------------------------------------
Time group                       |                                           
               ]enter; enter+84] |     1.000          .          .          . 
           ]enter+84; enter+365] |     1.617      0.258     10.123      0.607