Exercise 23. Calculating SMRs/SIRs
library(biostat3) # for Surv and survfit
library(dplyr)    # for data manipulation
library(tinyplot) # for some nice plotscalculate_smr = function(data)
    summarise(data,
              Observed = sum(observed),
              Expected = sum(expected)) |>
        mutate(SMR=Observed/Expected,
               poisson.ci(Observed,Expected))(a)
mel <- filter(biostat3::melanoma, stage == "Localised") |> 
        mutate( dead = (status %in% c("Dead: cancer","Dead: other") & surv_mm <= 120)+0, 
                surv_mm = pmin(120, surv_mm))
head(mel) ##      sex age     stage mmdx yydx surv_mm surv_yy       status       subsite
## 1 Female  81 Localised    2 1981    26.5     2.5  Dead: other Head and Neck
## 2 Female  75 Localised    9 1975    55.5     4.5  Dead: other Head and Neck
## 3 Female  78 Localised    2 1978   120.0    14.5  Dead: other         Limbs
## 4 Female  75 Localised    9 1975    19.5     1.5 Dead: cancer         Trunk
## 5 Female  75 Localised   10 1975    65.5     5.5  Dead: other Head and Neck
## 6 Female  80 Localised    6 1980    50.5     4.5  Dead: other Head and Neck
##          year8594         dx       exit agegrp id      ydx    yexit dead
## 1 Diagnosed 75-84 1981-02-02 1983-04-20    75+  1 1981.088 1983.298    1
## 2 Diagnosed 75-84 1975-09-21 1980-05-07    75+  2 1975.720 1980.348    1
## 3 Diagnosed 75-84 1978-02-21 1992-12-07    75+  3 1978.140 1992.934    0
## 4 Diagnosed 75-84 1975-09-03 1977-04-19    75+  6 1975.671 1977.296    1
## 5 Diagnosed 75-84 1975-10-14 1981-03-30    75+  7 1975.783 1981.241    1
## 6 Diagnosed 75-84 1980-06-11 1984-08-27    75+  8 1980.444 1984.654    1## Define the age at start and end of follow-up 
mel <- mutate(mel, adx = age+0.5,   # age at diagnosis  (mid-point approximation) 
              astart = adx, 
              astop  = adx+surv_mm/12)
## Split by age 
mel.split <- survSplit(mel, cut = 1:105, event = "dead", 
                       start = "astart", end = "astop")
## Quick check: the first two ids 
subset(mel.split, id<=2, select = c(id, astart, astop, dead)) ##   id astart    astop dead
## 1  1   81.5 82.00000    0
## 2  1   82.0 83.00000    0
## 3  1   83.0 83.70833    1
## 4  2   75.5 76.00000    0
## 5  2   76.0 77.00000    0
## 6  2   77.0 78.00000    0
## 7  2   78.0 79.00000    0
## 8  2   79.0 80.00000    0
## 9  2   80.0 80.12500    1(b)
# For each age time band from (a), we calculate the start and stop in calendar time 
# We calculate the time since diagnosis as difference between age at start/stop and 
# age at diagnosis, and add that interval to year at diagnosis
mel.split2 <- mutate(mel.split, 
                    ystart = ydx + astart - adx, 
                    ystop  = ydx + astop - adx)
subset(mel.split2, id<=2, select = c(id, adx, astart, astop, dead, ydx, ystart, ystop))##   id  adx astart    astop dead      ydx   ystart    ystop
## 1  1 81.5   81.5 82.00000    0 1981.088 1981.088 1981.588
## 2  1 81.5   82.0 83.00000    0 1981.088 1981.588 1982.588
## 3  1 81.5   83.0 83.70833    1 1981.088 1982.588 1983.296
## 4  2 75.5   75.5 76.00000    0 1975.720 1975.720 1976.220
## 5  2 75.5   76.0 77.00000    0 1975.720 1976.220 1977.220
## 6  2 75.5   77.0 78.00000    0 1975.720 1977.220 1978.220
## 7  2 75.5   78.0 79.00000    0 1975.720 1978.220 1979.220
## 8  2 75.5   79.0 80.00000    0 1975.720 1979.220 1980.220
## 9  2 75.5   80.0 80.12500    1 1975.720 1980.220 1980.345## Now we can split along the calendar time 
## For each of the new age-calendar time bands, we now have to adjust the age at 
## start and end in the same way as above 
mel.split2 <- survSplit( mel.split2, cut = 1970:2000, event = "dead", 
                         start = "ystart", end = "ystop" ) |>
              mutate( astart = adx + ystart - ydx, 
                      astop  = adx + ystop - ydx)
## Quick check: this seems ok 
subset(mel.split2, id<=2, select = c(id, ystart, ystop, astart, astop, dead)) ##    id   ystart    ystop   astart    astop dead
## 1   1 1981.088 1981.588 81.50000 82.00000    0
## 2   1 1981.588 1982.000 82.00000 82.41239    0
## 3   1 1982.000 1982.588 82.41239 83.00000    0
## 4   1 1982.588 1983.000 83.00000 83.41239    0
## 5   1 1983.000 1983.296 83.41239 83.70833    1
## 6   2 1975.720 1976.000 75.50000 75.77993    0
## 7   2 1976.000 1976.220 75.77993 76.00000    0
## 8   2 1976.220 1977.000 76.00000 76.77993    0
## 9   2 1977.000 1977.220 76.77993 77.00000    0
## 10  2 1977.220 1978.000 77.00000 77.77993    0
## 11  2 1978.000 1978.220 77.77993 78.00000    0
## 12  2 1978.220 1979.000 78.00000 78.77993    0
## 13  2 1979.000 1979.220 78.77993 79.00000    0
## 14  2 1979.220 1980.000 79.00000 79.77993    0
## 15  2 1980.000 1980.220 79.77993 80.00000    0
## 16  2 1980.220 1980.345 80.00000 80.12500    1(c)
## We calculate the total person time at risk for each time band
mel.split2 <- mutate( mel.split2, 
                      age  = floor(astart),  # Age at which person time was observed 
                      year = floor(ystart),  # Calendar year during which person time was observed 
                      pt  =  ystop - ystart  # ... or astop - astart, works the same
                    ) 
subset(mel.split2, id<=2, select = c(id, ystart, ystop, astart, astop, dead, age, year, pt))##    id   ystart    ystop   astart    astop dead age year        pt
## 1   1 1981.088 1981.588 81.50000 82.00000    0  81 1981 0.5000000
## 2   1 1981.588 1982.000 82.00000 82.41239    0  82 1981 0.4123864
## 3   1 1982.000 1982.588 82.41239 83.00000    0  82 1982 0.5876136
## 4   1 1982.588 1983.000 83.00000 83.41239    0  83 1982 0.4123864
## 5   1 1983.000 1983.296 83.41239 83.70833    1  83 1983 0.2959470
## 6   2 1975.720 1976.000 75.50000 75.77993    0  75 1975 0.2799255
## 7   2 1976.000 1976.220 75.77993 76.00000    0  75 1976 0.2200745
## 8   2 1976.220 1977.000 76.00000 76.77993    0  75 1976 0.7799255
## 9   2 1977.000 1977.220 76.77993 77.00000    0  76 1977 0.2200745
## 10  2 1977.220 1978.000 77.00000 77.77993    0  76 1977 0.7799255
## 11  2 1978.000 1978.220 77.77993 78.00000    0  77 1978 0.2200745
## 12  2 1978.220 1979.000 78.00000 78.77993    0  77 1978 0.7799255
## 13  2 1979.000 1979.220 78.77993 79.00000    0  78 1979 0.2200745
## 14  2 1979.220 1980.000 79.00000 79.77993    0  78 1979 0.7799255
## 15  2 1980.000 1980.220 79.77993 80.00000    0  79 1980 0.2200745
## 16  2 1980.220 1980.345 80.00000 80.12500    1  79 1980 0.1250000## Now tabulate: sum of person time across all combinations of age & year 
## (for some years, ages) 
xtabs(pt ~ age + year, data=mel.split2, subset = age>=50 & age<60 & year>=1980 & year<1990)##     year
## age      1980     1981     1982     1983     1984     1985     1986     1987
##   50 17.75739 16.10950 22.14116 21.53778 28.80944 35.87093 32.43368 36.35429
##   51 17.84611 21.10601 19.73652 25.74929 25.34206 33.81678 36.15502 36.99009
##   52 22.66005 23.20135 26.14500 23.00520 27.86951 26.86491 34.68995 38.43219
##   53 17.68110 27.56296 28.56208 30.58475 27.26166 31.79285 26.89462 36.09376
##   54 16.85111 19.77104 32.88876 33.04496 35.43436 30.98549 33.83041 28.67146
##   55 15.34311 23.06560 24.02609 35.78792 38.21690 37.88843 33.36533 37.75103
##   56 19.62087 17.07697 29.49875 29.69993 40.84652 40.83995 40.13706 38.25580
##   57 20.87159 23.91900 21.89273 33.33116 35.31722 44.38920 43.96677 41.94074
##   58 22.71497 23.32395 30.11496 28.08525 36.47228 38.15391 45.32161 47.98965
##   59 22.12377 29.11800 29.50430 35.09614 31.73203 40.31984 37.83583 46.55349
##     year
## age      1988     1989
##   50 43.11127 46.04857
##   51 39.73722 47.18098
##   52 37.43768 40.44656
##   53 42.05711 40.60234
##   54 42.20236 49.13991
##   55 33.63208 48.03838
##   56 45.57863 38.02925
##   57 42.72208 54.65228
##   58 45.11685 49.34837
##   59 54.18740 52.15563## Same: sum up 0/1 alive/dead for total count of deaths 
xtabs(dead ~ age + year, data=mel.split2, subset = age>=50 & age<60 & year>=1980 & year<1990) ##     year
## age  1980 1981 1982 1983 1984 1985 1986 1987 1988 1989
##   50    1    1    3    0    0    2    1    1    1    1
##   51    1    0    1    1    1    1    0    0    1    0
##   52    1    1    0    0    0    0    1    1    0    0
##   53    2    1    2    1    1    2    1    2    0    0
##   54    0    1    2    0    3    0    0    1    1    0
##   55    1    0    1    2    0    1    2    1    0    2
##   56    1    2    1    1    3    2    0    1    0    1
##   57    0    3    0    1    0    3    1    2    2    0
##   58    1    1    2    1    2    3    1    0    0    0
##   59    1    0    3    2    2    2    0    3    2    3(d)
mel.split2 <- mutate(mel.split2, 
                     age10  = cut(age, seq(0, 110 ,by=10), right=FALSE), 
                     year10 = cut(year, seq(1970, 2000, by=5), right=FALSE) 
                    ) 
sr <- survRate(Surv(pt, dead) ~ sex + age10 + year10, data=mel.split2) 
rownames(sr) <-1:nrow(sr) ## Simple rownames for display 
head(sr, n = 20) ##     sex   age10      year10      tstop event        rate        lower
## 1  Male  [0,10) [1980,1985)   6.509911     0 0.000000000 0.0000000000
## 2  Male  [0,10) [1985,1990)   4.490089     0 0.000000000 0.0000000000
## 3  Male [10,20) [1975,1980)   4.458333     1 0.224299065 0.0056787607
## 4  Male [10,20) [1980,1985)   6.842898     0 0.000000000 0.0000000000
## 5  Male [10,20) [1985,1990)  19.640948     0 0.000000000 0.0000000000
## 6  Male [10,20) [1990,1995)  24.219691     0 0.000000000 0.0000000000
## 7  Male [10,20) [1995,2000)   4.171463     0 0.000000000 0.0000000000
## 8  Male [20,30) [1975,1980)  26.398108     1 0.037881503 0.0009590766
## 9  Male [20,30) [1980,1985)  57.999389     4 0.068966244 0.0187909804
## 10 Male [20,30) [1985,1990) 143.213782     1 0.006982568 0.0001767833
## 11 Male [20,30) [1990,1995) 180.307318     3 0.016638260 0.0034312092
## 12 Male [20,30) [1995,2000)  27.414736     1 0.036476732 0.0009235109
## 13 Male [30,40) [1975,1980)  84.052956     6 0.071383569 0.0261965118
## 14 Male [30,40) [1980,1985) 292.510103     8 0.027349483 0.0118075654
## 15 Male [30,40) [1985,1990) 412.675104    15 0.036348207 0.0203438154
## 16 Male [30,40) [1990,1995) 434.040225     9 0.020735405 0.0094815477
## 17 Male [30,40) [1995,2000)  70.554946     3 0.042520053 0.0087686571
## 18 Male [40,50) [1975,1980) 160.375465    11 0.068589045 0.0342394041
## 19 Male [40,50) [1980,1985) 470.379381    15 0.031889153 0.0178481168
## 20 Male [40,50) [1985,1990) 846.907111    26 0.030699943 0.0200542220
##         upper
## 1  0.56665587
## 2  0.82156048
## 3  1.24971441
## 4  0.53908148
## 5  0.18781575
## 6  0.15230910
## 7  0.88431320
## 8  0.21106222
## 9  0.17658098
## 10 0.03890438
## 11 0.04862406
## 12 0.20323534
## 13 0.15537198
## 14 0.05388938
## 15 0.05995084
## 16 0.03936226
## 17 0.12426164
## 18 0.12272475
## 19 0.05259631
## 20 0.04498253(e)
pt <- mutate(mel.split2, sex = unclass(sex)) |>    # make sex integer to be in line with popmort 
    group_by(sex, age, year)                 |>    # aggregate by sex, age, year 
    summarise(pt = sum(pt), observed = sum(dead), .groups="keep") |> # sum the person time, deaths 
    ungroup()  # For convenience
head(pt) ## # A tibble: 6 × 5
##     sex   age  year    pt observed
##   <int> <dbl> <dbl> <dbl>    <dbl>
## 1     1     4  1980 0.5          0
## 2     1     4  1983 0.5          0
## 3     1     5  1980 0.213        0
## 4     1     5  1981 0.787        0
## 5     1     5  1983 0.297        0
## 6     1     5  1984 0.703        0##   sex    prob        rate age year
## 1   1 0.96429 0.036363177   0 1951
## 2   1 0.99639 0.003616547   1 1951
## 3   1 0.99783 0.002172384   2 1951
## 4   1 0.99842 0.001581249   3 1951
## 5   1 0.99882 0.001180690   4 1951
## 6   1 0.99893 0.001070595   5 1951##       sex           prob             rate               age       
##  Min.   :1.0   Min.   :0.5238   Min.   :0.000090   Min.   :  0.0  
##  1st Qu.:1.0   1st Qu.:0.9055   1st Qu.:0.001181   1st Qu.: 26.0  
##  Median :1.5   Median :0.9926   Median :0.007468   Median : 52.5  
##  Mean   :1.5   Mean   :0.9278   Mean   :0.084049   Mean   : 52.5  
##  3rd Qu.:2.0   3rd Qu.:0.9988   3rd Qu.:0.099218   3rd Qu.: 79.0  
##  Max.   :2.0   Max.   :0.9999   Max.   :0.646626   Max.   :105.0  
##       year     
##  Min.   :1951  
##  1st Qu.:1963  
##  Median :1976  
##  Mean   :1976  
##  3rd Qu.:1988  
##  Max.   :2000## Joining with `by = join_by(sex, age, year)`## # A tibble: 6 × 7
##     sex   age  year    pt observed  prob     rate
##   <int> <dbl> <dbl> <dbl>    <dbl> <dbl>    <dbl>
## 1     1     4  1980 0.5          0  1.00 0.000410
## 2     1     4  1983 0.5          0  1.00 0.000280
## 3     1     5  1980 0.213        0  1.00 0.000350
## 4     1     5  1981 0.787        0  1.00 0.000250
## 5     1     5  1983 0.297        0  1.00 0.000250
## 6     1     5  1984 0.703        0  1.00 0.000250## # A tibble: 6 × 8
##     sex   age  year    pt observed  prob     rate  expected
##   <int> <dbl> <dbl> <dbl>    <dbl> <dbl>    <dbl>     <dbl>
## 1     1     4  1980 0.5          0  1.00 0.000410 0.000205 
## 2     1     4  1983 0.5          0  1.00 0.000280 0.000140 
## 3     1     5  1980 0.213        0  1.00 0.000350 0.0000744
## 4     1     5  1981 0.787        0  1.00 0.000250 0.000197 
## 5     1     5  1983 0.297        0  1.00 0.000250 0.0000744
## 6     1     5  1984 0.703        0  1.00 0.000250 0.000176(f)
| Observed | Expected | SMR | 2.5 % | 97.5 % | 
|---|---|---|---|---|
| 1580 | 775.6643 | 2.036964 | 1.937751 | 2.139939 | 
| sex | Observed | Expected | SMR | 2.5 % | 97.5 % | 
|---|---|---|---|---|---|
| 1 | 818 | 394.1975 | 2.075102 | 1.935316 | 2.222316 | 
| 2 | 762 | 381.4668 | 1.997553 | 1.858221 | 2.144564 | 
| year | Observed | Expected | SMR | 2.5 % | 97.5 % | 
|---|---|---|---|---|---|
| 1975 | 2 | 1.861354 | 1.074486 | 0.1301253 | 3.881415 | 
| 1976 | 14 | 4.818816 | 2.905278 | 1.5883425 | 4.874563 | 
| 1977 | 31 | 7.004271 | 4.425871 | 3.0071647 | 6.282171 | 
| 1978 | 34 | 10.789242 | 3.151287 | 2.1823579 | 4.403608 | 
| 1979 | 44 | 15.461017 | 2.845867 | 2.0678114 | 3.820444 | 
| 1980 | 52 | 19.133171 | 2.717793 | 2.0297778 | 3.564024 | 
| 1981 | 53 | 21.050479 | 2.517757 | 1.8859731 | 3.293289 | 
| 1982 | 79 | 26.296027 | 3.004256 | 2.3784978 | 3.744200 | 
| 1983 | 75 | 31.357280 | 2.391789 | 1.8812939 | 2.998128 | 
| 1984 | 83 | 36.451780 | 2.276981 | 1.8136025 | 2.822660 | 
| 1985 | 79 | 39.894258 | 1.980235 | 1.5677706 | 2.467963 | 
| 1986 | 79 | 42.358990 | 1.865011 | 1.4765471 | 2.324361 | 
| 1987 | 86 | 47.820352 | 1.798398 | 1.4384858 | 2.221006 | 
| 1988 | 108 | 52.241548 | 2.067320 | 1.6958627 | 2.495958 | 
| 1989 | 99 | 55.155101 | 1.794938 | 1.4588372 | 2.185273 | 
| 1990 | 104 | 57.597529 | 1.805633 | 1.4753330 | 2.187827 | 
| 1991 | 100 | 56.277542 | 1.776908 | 1.4457631 | 2.161196 | 
| 1992 | 105 | 59.106054 | 1.776468 | 1.4529746 | 2.150523 | 
| 1993 | 126 | 62.302058 | 2.022405 | 1.6847175 | 2.407931 | 
| 1994 | 93 | 64.749787 | 1.436298 | 1.1592780 | 1.759562 | 
| 1995 | 134 | 63.937681 | 2.095791 | 1.7559794 | 2.482176 | 
joint <- mutate(joint, age_group = cut(age, seq(0, 110, by=10), right = FALSE))
SMR_byAge <- group_by(joint, age_group) |> calculate_smr()
SMR_byAge |> kable("html")| age_group | Observed | Expected | SMR | 2.5 % | 97.5 % | 
|---|---|---|---|---|---|
| [0,10) | 0 | 0.0046936 | 0.000000 | 0.0000000 | 785.943356 | 
| [10,20) | 2 | 0.0688809 | 29.035636 | 3.5163502 | 104.886700 | 
| [20,30) | 16 | 0.9887579 | 16.181919 | 9.2493647 | 26.278422 | 
| [30,40) | 73 | 4.9153573 | 14.851413 | 11.6411368 | 18.673422 | 
| [40,50) | 144 | 19.1854079 | 7.505704 | 6.3298848 | 8.836586 | 
| [50,60) | 215 | 51.7963269 | 4.150874 | 3.6145389 | 4.744364 | 
| [60,70) | 323 | 129.0811235 | 2.502302 | 2.2368231 | 2.790623 | 
| [70,80) | 390 | 257.7721031 | 1.512964 | 1.3665152 | 1.670833 | 
| [80,90) | 342 | 250.9696149 | 1.362715 | 1.2221008 | 1.515071 | 
| [90,100) | 74 | 60.0807906 | 1.231675 | 0.9671297 | 1.546255 | 
| [100,110) | 1 | 0.8012821 | 1.248000 | 0.0315966 | 6.953411 | 
plot(SMR~age_group, SMR_byAge, xlab="Age group (years)")
abline( h = 1:2, lty = 2)  # two reference lines at 1 & 2## joint$sex: 1
## 
##  Exact Poisson test
## 
## data:  sum(data$observed) time base: sum(data$expected)
## number of events = 818, time base = 394.2, p-value < 2.2e-16
## alternative hypothesis: true event rate is not equal to 1
## 95 percent confidence interval:
##  1.935316 2.222317
## sample estimates:
## event rate 
##   2.075102 
## 
## ------------------------------------------------------------ 
## joint$sex: 2
## 
##  Exact Poisson test
## 
## data:  sum(data$observed) time base: sum(data$expected)
## number of events = 762, time base = 381.47, p-value < 2.2e-16
## alternative hypothesis: true event rate is not equal to 1
## 95 percent confidence interval:
##  1.858221 2.144564
## sample estimates:
## event rate 
##   1.997553(g)
joint2 <- transform(joint,
                    sex  = factor(sex, levels = 1:2, labels = c("m", "f")),
                    year =  factor(year) |> relevel(ref = "1985"),  # mid-study
                    age_group = relevel(age_group, ref = "[70,80)")  # Close to one already
                   )
## Model & parameters
summary(fit <- glm(observed ~ sex + year + age_group + offset(log(expected)), data=joint2, family=poisson)) ## 
## Call:
## glm(formula = observed ~ sex + year + age_group + offset(log(expected)), 
##     family = poisson, data = joint2)
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          0.24931    0.12399   2.011 0.044354 *  
## sexf                 0.19732    0.05232   3.771 0.000162 ***
## year1975            -0.73554    0.71610  -1.027 0.304348    
## year1976             0.29944    0.29010   1.032 0.301975    
## year1977             0.69624    0.21205   3.283 0.001026 ** 
## year1978             0.39484    0.20517   1.925 0.054291 .  
## year1979             0.32616    0.18814   1.734 0.082996 .  
## year1980             0.28623    0.17860   1.603 0.109023    
## year1981             0.22373    0.17758   1.260 0.207711    
## year1982             0.40014    0.15913   2.515 0.011917 *  
## year1983             0.17359    0.16124   1.077 0.281660    
## year1984             0.13835    0.15719   0.880 0.378769    
## year1986            -0.04542    0.15913  -0.285 0.775291    
## year1987            -0.07071    0.15586  -0.454 0.650058    
## year1988             0.07989    0.14809   0.539 0.589558    
## year1989            -0.05832    0.15090  -0.386 0.699130    
## year1990            -0.05248    0.14929  -0.352 0.725201    
## year1991            -0.04901    0.15068  -0.325 0.744989    
## year1992            -0.03529    0.14910  -0.237 0.812920    
## year1993             0.10860    0.14366   0.756 0.449680    
## year1994            -0.22877    0.15316  -1.494 0.135246    
## year1995             0.16654    0.14207   1.172 0.241115    
## age_group[0,10)     -7.66192  317.84697  -0.024 0.980768    
## age_group[10,20)     2.95588    0.70911   4.168 3.07e-05 ***
## age_group[20,30)     2.38059    0.25515   9.330  < 2e-16 ***
## age_group[30,40)     2.29433    0.12768  17.969  < 2e-16 ***
## age_group[40,50)     1.63236    0.09788  16.677  < 2e-16 ***
## age_group[50,60)     1.03105    0.08554  12.053  < 2e-16 ***
## age_group[60,70)     0.52924    0.07571   6.990 2.74e-12 ***
## age_group[80,90)    -0.11656    0.07472  -1.560 0.118759    
## age_group[90,100)   -0.21350    0.12769  -1.672 0.094519 .  
## age_group[100,110)  -0.18214    1.00359  -0.181 0.855982    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 3205.1  on 3268  degrees of freedom
## Residual deviance: 2540.6  on 3237  degrees of freedom
## AIC: 5037.3
## 
## Number of Fisher Scoring iterations: 14##                       exp(beta)         2.5 %        97.5 %
## (Intercept)        1.283137e+00  1.006314e+00  1.636109e+00
## sexf               1.218128e+00  1.099401e+00  1.349677e+00
## year1975           4.792457e-01  1.177646e-01  1.950300e+00
## year1976           1.349103e+00  7.640366e-01  2.382189e+00
## year1977           2.006194e+00  1.323953e+00  3.039998e+00
## year1978           1.484151e+00  9.927525e-01  2.218786e+00
## year1979           1.385637e+00  9.582975e-01  2.003542e+00
## year1980           1.331396e+00  9.381684e-01  1.889442e+00
## year1981           1.250729e+00  8.830989e-01  1.771403e+00
## year1982           1.492038e+00  1.092268e+00  2.038124e+00
## year1983           1.189572e+00  8.672427e-01  1.631702e+00
## year1984           1.148381e+00  8.438903e-01  1.562737e+00
## year1986           9.555917e-01  6.995579e-01  1.305332e+00
## year1987           9.317305e-01  6.864679e-01  1.264621e+00
## year1988           1.083166e+00  8.102976e-01  1.447924e+00
## year1989           9.433463e-01  7.018197e-01  1.267993e+00
## year1990           9.488769e-01  7.081704e-01  1.271399e+00
## year1991           9.521713e-01  7.086851e-01  1.279313e+00
## year1992           9.653299e-01  7.207160e-01  1.292967e+00
## year1993           1.114719e+00  8.411598e-01  1.477245e+00
## year1994           7.955094e-01  5.892228e-01  1.074017e+00
## year1995           1.181209e+00  8.941165e-01  1.560483e+00
## age_group[0,10)    4.704029e-04 1.320197e-274 1.676105e+267
## age_group[10,20)   1.921857e+01  4.787698e+00  7.714632e+01
## age_group[20,30)   1.081125e+01  6.556804e+00  1.782625e+01
## age_group[30,40)   9.917758e+00  7.722012e+00  1.273786e+01
## age_group[40,50)   5.115943e+00  4.222904e+00  6.197837e+00
## age_group[50,60)   2.804008e+00  2.371198e+00  3.315818e+00
## age_group[60,70)   1.697641e+00  1.463531e+00  1.969199e+00
## age_group[80,90)   8.899779e-01  7.687412e-01  1.030335e+00
## age_group[90,100)  8.077550e-01  6.289151e-01  1.037450e+00
## age_group[100,110) 8.334822e-01  1.165831e-01  5.958774e+00## Single term deletions
## 
## Model:
## observed ~ sex + year + age_group + offset(log(expected))
##           Df Deviance    AIC    LRT  Pr(>Chi)    
## <none>         2540.6 5037.3                     
## sex        1   2554.8 5049.5  14.17 0.0001673 ***
## year      20   2590.8 5047.5  50.19 0.0002079 ***
## age_group 10   3135.6 5612.3 594.95 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1