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