# ------------------------------------------------------------
#
# Analysoidaan viikottaisia kuolleisuuksia
# Created: Thu May 12 13:46:55 2022
# Author: Jari Haukka
#
# ------------------------------------------------------------
library(Epi)
library(pxR)
library(ggplot2)
library(chron)
library(shinyPredict)
# ------------------------------------------------------------
# Read data from Statics Finland
# ------------------------------------------------------------
# apu.dt<-read.px(filename="https://pxnet2.stat.fi:443/PXWeb/sq/786c9a7f-3081-4508-a071-18701dfae142")
apu.dt<-read.px(filename="https://pxdata.stat.fi:443/PxWeb/sq/1ccda3d4-1740-4e1e-9c64-2a478eeb4a55")
apu.dt<-as.data.frame(apu.dt)
# Drop last two weeks
apu.last.2<-as.character(rev(sort(unique(apu.dt$Viikko)))[1:2])
apu.dt<-subset(apu.dt,!(Viikko%in%apu.last.2))
# Add date for each week mid
apu.yr<-as.numeric(substr(apu.dt$Viikko,1,4))
apu.vk<-as.numeric(substr(apu.dt$Viikko,6,7))
apu.dt$pvm<-apu.yr+(apu.vk*7-3.5)/365
apu.dt$vk<-apu.vk
apu.dt$Vuosi<-apu.yr
# Drop week 53 that rarely exists
apu.dt<-subset( apu.dt,vk!=53)
# Väestö 2003 alkaen
# Update in start of year
apu.popul<-read.px(filename="https://pxweb2.stat.fi:443/PxWeb/sq/79514db7-5411-4d58-8b07-3dc290e7cc9a")
apu.popul<-as.data.frame(apu.popul)
# table(Relevel(apu.popul$Ikä,list(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19:20)))
apu.popul.1<-with(apu.popul,aggregate(list(popul=value),
list(Ikä=Relevel(Ikä,list(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19:20)),
Sukupuoli=Sukupuoli,Vuosi=Vuosi),
sum))
levels(apu.popul.1$Ikä)<-levels(apu.dt$Ikä)
apu.popul.1$Vuosi<-as.numeric(as.character(apu.popul.1$Vuosi))
apu.lisa<-subset(apu.popul.1,Vuosi==2021)
apu.lisa$Vuosi<-2022
apu.popul.1<-rbind(apu.popul.1,apu.lisa)
# save(apu.popul.1,file="apupopul1.RData")
# load(file="apupopul1.RData")
Suomi.vk.mort.V3<-merge(x=apu.dt,
y=apu.popul.1,
by=c("Vuosi","Sukupuoli","Ikä"),
all.x=TRUE,all.y=TRUE)
# dim(apu.dt);dim(apu.popul.1);dim(Suomi.vk.mort.V3)
# head(subset(Suomi.vk.mort.V3,Vuosi==2022),4)
Suomi.vk.mort.V3$value.milj<-Suomi.vk.mort.V3$value/(Suomi.vk.mort.V3$popul/1000000)
Suomi.vk.mort.V3<-subset(Suomi.vk.mort.V3,(Vuosi>=2003))#&(Vuosi<2022) )
Suomi.vk.mort.V3$Vuosi.f<-factor(Suomi.vk.mort.V3$Vuosi)
# Mallitus
# Vuosittaista kausivaihtelu kuvataan sin ja cos funktioiden yhdistelmällä
tmp.m1<-glm(I(value/popul)~sin(2*pi*(vk)/52)+cos(2*pi*(vk)/52)+Sukupuoli+Ikä,
data=Suomi.vk.mort.V3,
family=binomial,weights=popul)
# Lisätään vuosi
tmp.m2<-update(tmp.m1,~.+Relevel(Vuosi.f,ref=17))
anova(tmp.m1,tmp.m2,test="Chisq")## Analysis of Deviance Table
##
## Model 1: I(value/popul) ~ sin(2 * pi * (vk)/52) + cos(2 * pi * (vk)/52) +
## Sukupuoli + Ikä
## Model 2: I(value/popul) ~ sin(2 * pi * (vk)/52) + cos(2 * pi * (vk)/52) +
## Sukupuoli + Ikä + Relevel(Vuosi.f, ref = 17)
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 38814 66638
## 2 38795 58193 19 8445.3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Vuoden ja iän yhdysvaikutus
tmp.m3<-update(tmp.m2,~.+Vuosi.f:Ikä)
anova(tmp.m2,tmp.m3,test="Chisq")## Analysis of Deviance Table
##
## Model 1: I(value/popul) ~ sin(2 * pi * (vk)/52) + cos(2 * pi * (vk)/52) +
## Sukupuoli + Ikä + Relevel(Vuosi.f, ref = 17)
## Model 2: I(value/popul) ~ sin(2 * pi * (vk)/52) + cos(2 * pi * (vk)/52) +
## Sukupuoli + Ikä + Relevel(Vuosi.f, ref = 17) + Ikä:Vuosi.f
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 38795 58193
## 2 38453 55722 342 2470.4 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Vuoden ja kausivaihtelun yhdysvaikutus
# tmp.m4<-update(tmp.m3,~.+Vuosi.f:sin(2*pi*(vk)/52)+Vuosi.f:cos(2*pi*(vk)/52))
# anova(tmp.m3,tmp.m4,test="Chisq")
# Mallit on esitetty käyttäen `shinyPredict`-funktiota [`shinyPredict` paketista](https://cran.r-project.org/web/packages/shinyPredict/index.html).
# shinyPredict(list(tmp.m1,tmp.m2,tmp.m3,tmp.m4),data =Suomi.vk.mort.V3[,c("vk","Vuosi.f","Sukupuoli","Ikä")],path = "./" )knitr::kable(ci.exp(tmp.m2,subset=-1), caption="Vaarasuhteet (odds ratio) ja 95% luottamusväli päävaikutusmallista, vertailuvuotena 2019",digits = 3)| exp(Est.) | 2.5% | 97.5% | |
|---|---|---|---|
| sin(2 * pi * (vk)/52) | 1.043 | 1.040 | 1.046 |
| cos(2 * pi * (vk)/52) | 1.062 | 1.059 | 1.065 |
| SukupuoliNaiset | 0.605 | 0.602 | 0.607 |
| Ikä5 - 9 | 0.154 | 0.140 | 0.169 |
| Ikä10 - 14 | 0.180 | 0.165 | 0.196 |
| Ikä15 - 19 | 0.650 | 0.616 | 0.685 |
| Ikä20 - 24 | 1.065 | 1.017 | 1.116 |
| Ikä25 - 29 | 1.125 | 1.075 | 1.178 |
| Ikä30 - 34 | 1.302 | 1.245 | 1.361 |
| Ikä35 - 39 | 1.755 | 1.683 | 1.830 |
| Ikä40 - 44 | 2.629 | 2.527 | 2.735 |
| Ikä45 - 49 | 4.216 | 4.060 | 4.378 |
| Ikä50 - 54 | 6.661 | 6.423 | 6.909 |
| Ikä55 - 59 | 10.125 | 9.769 | 10.494 |
| Ikä60 - 64 | 15.147 | 14.620 | 15.694 |
| Ikä65 - 69 | 22.651 | 21.867 | 23.464 |
| Ikä70 - 74 | 35.357 | 34.137 | 36.620 |
| Ikä75 - 79 | 60.735 | 58.647 | 62.896 |
| Ikä80 - 84 | 112.315 | 108.465 | 116.301 |
| Ikä85 - 89 | 215.481 | 208.098 | 223.127 |
| Ikä90 - | 460.995 | 445.193 | 477.358 |
| Relevel(Vuosi.f, ref = 17)2003 | 1.382 | 1.366 | 1.400 |
| Relevel(Vuosi.f, ref = 17)2004 | 1.312 | 1.296 | 1.328 |
| Relevel(Vuosi.f, ref = 17)2005 | 1.281 | 1.266 | 1.297 |
| Relevel(Vuosi.f, ref = 17)2006 | 1.247 | 1.232 | 1.263 |
| Relevel(Vuosi.f, ref = 17)2007 | 1.240 | 1.225 | 1.255 |
| Relevel(Vuosi.f, ref = 17)2008 | 1.199 | 1.185 | 1.214 |
| Relevel(Vuosi.f, ref = 17)2009 | 1.196 | 1.181 | 1.210 |
| Relevel(Vuosi.f, ref = 17)2010 | 1.183 | 1.169 | 1.198 |
| Relevel(Vuosi.f, ref = 17)2011 | 1.140 | 1.127 | 1.154 |
| Relevel(Vuosi.f, ref = 17)2012 | 1.134 | 1.120 | 1.148 |
| Relevel(Vuosi.f, ref = 17)2013 | 1.103 | 1.090 | 1.116 |
| Relevel(Vuosi.f, ref = 17)2014 | 1.090 | 1.077 | 1.103 |
| Relevel(Vuosi.f, ref = 17)2015 | 1.069 | 1.056 | 1.082 |
| Relevel(Vuosi.f, ref = 17)2016 | 1.067 | 1.054 | 1.079 |
| Relevel(Vuosi.f, ref = 17)2017 | 1.042 | 1.030 | 1.055 |
| Relevel(Vuosi.f, ref = 17)2018 | 1.036 | 1.024 | 1.048 |
| Relevel(Vuosi.f, ref = 17)2020 | 0.999 | 0.987 | 1.011 |
| Relevel(Vuosi.f, ref = 17)2021 | 1.018 | 1.006 | 1.030 |
| Relevel(Vuosi.f, ref = 17)2022 | 1.086 | 1.072 | 1.100 |
apu.RR<-ci.exp(tmp.m2,subset="Vuosi.*2022")
apu.RR.txt<-paste(c(format(apu.RR[1],digits = 3)," (",format(apu.RR[2],digits = 3),
"-",format(apu.RR[3],digits = 3),")"),collapse = "")Alla olevasta taulukosta voidaan todeta, että kuoleman todennäköisyys vuonna 2022 on noin 8.56% korkeampi kuin vuonna 2019.
Suhteellinen riski mallin perusteella on 1.09 (1.07-1.1) (95% luottamusväli).
Tuloksissa on otettu huomioon ikä, sukupuoli ja vuodenaikaisvaihtelu.
Kuolemien lukumäärät viikossa ovat peräisen Tilastokeskuksen tietokannasta. Viimeistä kahta viikko ei ole sisällytetty aineistoon lukujen epävarmuuden takia. Aineisto alkaa vuoden 2003 alusta.
Keskiväestön koko ikäryhmittäin on saatu Tilastokeskuksen tietokannasta. Vuoden 2022 väestönä on käytetty vuoden 2021 väestöä.
Kukin viikon kuolemien lukumäärää mallitettin logistiselle regressiolla, jossa nimittäjänä on keskiväestön koko. Mallissa otettin huomioon vuosittainen kausivaihtelu (sini ja kosini termit), sukupuoli, ikä ja kalenterivuosi. Tulokset esitetään vaarasuhteina (odds ratios), joka on hyvä estimaatti riskille.
sessionInfo()## R version 4.2.1 (2022-06-23)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur ... 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] fi_FI.UTF-8/fi_FI.UTF-8/fi_FI.UTF-8/C/fi_FI.UTF-8/fi_FI.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] shinyPredict_0.1.1 chron_2.3-58 ggplot2_3.3.6 pxR_0.42.4
## [5] plyr_1.8.7 RJSONIO_1.3-1.6 reshape2_1.4.4 stringr_1.4.1
## [9] Epi_2.47
##
## loaded via a namespace (and not attached):
## [1] zoo_1.8-10 tidyselect_1.1.2 xfun_0.32
## [4] bslib_0.4.0 purrr_0.3.4 splines_4.2.1
## [7] etm_1.1.1 lattice_0.20-45 colorspace_2.0-3
## [10] vctrs_0.4.1 generics_0.1.3 htmltools_0.5.3
## [13] yaml_2.3.5 mgcv_1.8-40 utf8_1.2.2
## [16] survival_3.4-0 rlang_1.0.5 jquerylib_0.1.4
## [19] pillar_1.8.1 withr_2.5.0 glue_1.6.2
## [22] lifecycle_1.0.1 munsell_0.5.0 gtable_0.3.1
## [25] evaluate_0.16 knitr_1.40 fastmap_1.1.0
## [28] parallel_4.2.1 fansi_1.0.3 highr_0.9
## [31] Rcpp_1.0.9 scales_1.2.1 cachem_1.0.6
## [34] cmprsk_2.2-11 jsonlite_1.8.0 digest_0.6.29
## [37] stringi_1.7.8 dplyr_1.0.10 numDeriv_2016.8-1.1
## [40] grid_4.2.1 cli_3.3.0 tools_4.2.1
## [43] magrittr_2.0.3 sass_0.4.2 tibble_3.1.8
## [46] pkgconfig_2.0.3 MASS_7.3-57 Matrix_1.4-1
## [49] data.table_1.14.2 rmarkdown_2.16 rstudioapi_0.14
## [52] R6_2.5.1 nlme_3.1-157 compiler_4.2.1