standardsccs: The standard SCCS method

Description

Fits the standard SCCS model where age and exposure effects are represented by piecewise constant functions.

Usage

standardsccs(formula, indiv, astart, aend, aevent, adrug, aedrug,
            expogrp = list(), washout = list(), sameexpopar = list(),
            agegrp = NULL, seasongrp=NULL, dob=NULL, dataformat="stack",
            data)

Value

Relative incidence estimates along with their 95% confidence limits.

Arguments

formula: model formula. The dependent variable should always be "event" e.g. event ~ itp. If age and/or season effects are included, they should always be included as 'age' and 'season', e.g. event ~ itp + age.
indiv: a vector of individual identifiers of cases.
astart: a vector of ages at which the observation periods start.
aend: a vector of ages at end of observation periods.
aevent: a vector of ages at event; an individual can experience multiple events
adrug: a list of vectors of ages at start of exposures or a list of matrices if the exposures have multiple episodes (dataformat multi). Multiple exposures of the same type can be recorded as multiple rows (dataformat stack). One list item per exposure type.
aedrug: a list of vectors of ages at which exposure-related risk ends or a list of matrices if there are multiple episodes (repeat exposures in different columns) of the same exposure type. The dimension of each item of aedrug has to be equal to that of adrug, that is aedrug should be given for each exposure in adrug.
expogrp: list of vectors of days to the start of exposure-related risk, counted from adrug. E.g if the risk period is [adrug+c,aedrug], use expogrp = c or expogrp = list(c). For multiple exposure types expogrp is a list of the same length as list adrug. The DEFAULT is a list of zeros where the exposure-related risk periods are [adrug, aedrug].
washout: list of vectors with days to start of washout periods counted from aedrug, the number of vectors in the list is equal to the number of exposures or the length of list adrug. The default is NULL, no washout periods. The order of the list corresponds to the order of exposures in adrug.
sameexpopar: a vector of logical values. If TRUE (the default) no dose effect is assumed: the same exposure parameters are used for multiple doses/episodes of the same exposure type, presented in dataformat 'multi'. If FALSE different relative incidences are estimated for different doses of the same exposure. The length of the vector is equal to the length of list adrug. The order of the elements of the vector corrosponds to the order of exposures in list adrug.
agegrp: a vector of cut points of the age groups where each value represents the start of an age category. The first element in the vector is the start of the second age group. The first age group starts at the minimum of astart, the start of the observation period. The default is NULL (i.e no age effects included).
seasongrp: a vector of cut points for seasonal effects. The values should be given in ddmm format, representing the first days of each season group. The seasonal effect is a factor, the reference level being the time interval starting at the earliest date in seasongrp. The default is NULL where no seasonal effects are included.
dob: a vector of birth dates of the cases, in ddmmyyyy format. They are used if seasonal effects are included in the model. The default dob is NULL. Required if seasongrp is not NULL.
dataformat: the way the input data are assembled. It accepts "multi" or "stack" (the default), where "multi" refers to a data assembled with one row representing one event and "stack" refers to a data frame where repeated exposures of the same type are stack in one column. In the "multi" dataformat different episodes of the same exposure type are recorded as separate columns in the dataframe.
data: a data frame containing the input data. The data should be in 'stack' or 'multi' (see dataformat).

Author

Yonas Ghebremichael-Weldeselassie, Heather Whitaker, Paddy Farrington.

Details

In the standard SCCS model, originally described in Farrington (1995), age and exposure effects are represented by step functions. Suppose that individual $i$ has $n_i$ events, $n_{ijk}$ occurring in age group $j$ and exposure group $k$ and that the $s^{{th}}$ event falls within age group $j_s$ and exposure group $k_s$. The SCCS likelihood contribution for individual $i$ is

$$l_i = \frac{\prod_{s=1}^{n_i} \exp(\alpha_{j_s} + \beta_{k_s}) }{\Big(\sum_{j=1}^J\sum_{k=1}^K\exp(\alpha_j + \beta_k) e_{ijk}\Big)^{n_i}}.$$

This is a multinomial likelihood with index $n_i$, responses $n_{ijk}$ and probabilities

$$p_{iuv} = \frac{\exp(\alpha_{u} + \beta_{v})e_{iuv} }{\sum_{j=1}^J \sum_{k=1}^K\exp(\alpha_j + \beta_k) e_{ijk}}.$$

The standard SCCS likelihood is equivalent to a product multinomial likelihood.

The SCCS likelihood is equivalent to that of the conditional logistic model for $1:M$ matched case-control studies: see Breslow and Day (1980), Chapter 7. This means that SCCS models can be fit using conditional logistic regression, with a factor for each individual event.

The SCCS R package maximises the likelihood using clogit. Each event is assigned an identifier (indivL), and direct estimation is avoided using the option strata(indivL).

References

Farrington, C.P. (1995). Relative incidence estimation from case series for vaccine safety evaluation. Biometrics 51, 228--235.

Breslow, N. E. and Day, N. E. (1980). Statistical Methods in Cancer Research, volume I: The analysis of case-control studies. IARC Publications No.32.

Farrington, P., Whitaker, H., and Ghebremichael-Weldeselassie, Y. (2018). Self-controlled Case Series studies: A modelling Guide with R. Boca Raton: Chapman & Hall/CRC Press.

Examples

Run this code


# Single exposure-related risk period with no age effect

  itp.mod1 <- standardsccs(event~mmr, indiv=case, astart=sta, aend=end,
                        aevent=itp, adrug=mmr, aedrug=mmr+42, data=itpdat)

  itp.mod1

# Single exposure-related risk period and age effect included

  itp.mod2 <- standardsccs(event~mmr+age, indiv=case, astart=sta, aend=end,
                        aevent=itp, adrug=mmr, aedrug=mmr+42,
                        agegrp=c(427,488,549,610,671), data=itpdat)

  itp.mod2

# Multiple risk periods and age effect included  
  
  itp.mod3 <- standardsccs(event~mmr+age, indiv=case, astart=sta, aend=end,
                        aevent=itp, adrug=mmr, aedrug=mmr+42, expogrp=c(0,15,29),
                        agegrp=c(427,488,549,610,671), data=itpdat)

  itp.mod3

# Multiple risk periods, washout periods and age effects

   ageq <- floor(quantile(hipdat$frac, seq(0.05,0.95,0.05), 
                         names=FALSE)) # Age group 
                                      # cut points
   
   hip.mod1 <- standardsccs(event~ad+age, indiv=case, astart=sta, aend=end,
                           aevent=frac, adrug=ad, aedrug=endad, expogrp=c(0,15,43),
                           washout=c(1,92,182), agegrp=ageq, data=hipdat)


# Multiple/repeat exposures of the same exposure type, dataformat="stack" 

   ageq <- floor(quantile(gidat$bleed[duplicated(gidat$case)==0], 
                seq(0.025,0.975,0.025), names=FALSE))

   gi.mod1 <- standardsccs(event~ns+relevel(age,ref=21), indiv=case, astart=sta,
                         aend=end, aevent=bleed, adrug=ns, aedrug=endns,
                         agegrp=ageq, dataformat="stack", data=gidat)

   gi.mod1
   
#  Multiple doses of a vaccine each with different parameter estimates (sameexpopar=F)
  
   
   ageg <- c(57,85,113,141,169,197,225,253,281,309,337) # age group cut points
   
   dtp.mod2 <- standardsccs(event~dtp+age, indiv=case, astart=sta, aend=end,
                           aevent=conv, adrug=cbind(dtp,dtpd2, dtpd3), 
                           aedrug=cbind(dtp+14,dtpd2+14,dtpd3+14), 
                           expogrp=c(0,4,8),agegrp=ageg, dataformat="multi", 
                           sameexpopar=FALSE, data=dtpdat)

   dtp.mod2
   
# Multiple exposure types 

  ageg <- seq(387,707,20) # Age group cut points
 con.mod <- standardsccs(event~hib+mmr+age, indiv=case, astart=sta, aend=end,
                          aevent=conv, adrug=cbind(hib,mmr), aedrug=cbind(hib+14,mmr+14), 
                          expogrp=list(c(0,8), c(0,8)), agegrp=ageg, data=condat)

  con.mod

   
# Multiple doses/episodes of several exposure types, the doses of each exposure type
# have same paramter


   ageg <- c(57,85,113,141,169,197,225,253,281,309,337)   # age group cut points 

   hib.mod1 <- standardsccs(event~dtp+hib+age, indiv=case, astart=sta, 
                            aend=end, aevent=conv, 
                           adrug=list(cbind(dtp, dtpd2,dtpd3),
                                       cbind(hib,hibd2,hibd3)),
                            aedrug=list(cbind(dtp+14,dtpd2+14,dtpd3+14),
                                        cbind(hib+14,hibd2+14,hibd3+14)),
                            expogrp=list(c(0,4,8),c(0,8)),agegrp=ageg, 
                            dataformat="multi", data=hibdat)
                            
   hib.mod1                   


# Multiple doses/episodes of several exposure types, the doses of "dtp" 
# different parameters and the doses of the second exposure hib have
# same paramters


   ageg <- c(57,85,113,141,169,197,225,253,281,309,337)
                    # age group cut points 

   hib.mod2 <- standardsccs(event~dtp+hib+age, 
                                   indiv=case, astart=sta, aend=end,
                        aevent=conv, adrug=list(cbind(dtp,dtpd2,dtpd3), 
                                      cbind(hib,hibd2,hibd3)), 
                        aedrug=list(cbind(dtp+3,dtpd2+3,dtpd3+3),
                                      cbind(hib+7,hibd2+7,hibd3+7)),
                        sameexpopar=c(FALSE,TRUE), agegrp=ageg, 
                        dataformat="multi", data=hibdat)
   hib.mod2

   
# Season included in a model 

   month <- c(0101,0102,0103,0104,0105,0106,0107,0108,0109,0110,0111,0112) 
   # season cutpoints 
   
   int.mod <- standardsccs(event~opv+age+season, indiv=case, astart=sta,
                          aend=end, aevent=intus, adrug=cbind(opv,opvd2),
                          aedrug=cbind(opv+42,opvd2+42), expogrp=c(0,15,29),
                          agegrp=seq(30,330,30), seasongrp=month,dob=dob, 
                          dataformat="multi", data=intdat)
   int.mod

Run the code above in your browser using DataLab