[med-svn] [r-cran-epi] 01/07: New upstream version 2.19
Andreas Tille
tille at debian.org
Mon Oct 9 15:16:49 UTC 2017
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository r-cran-epi.
commit 467a8b3e9767e5d0f0c049f663c4fb13230b25d1
Author: Andreas Tille <tille at debian.org>
Date: Mon Oct 9 16:42:58 2017 +0200
New upstream version 2.19
---
CHANGES | 89 ++++++-
DESCRIPTION | 18 +-
MD5 | 101 ++++----
NAMESPACE | 13 +-
R/LCa.fit.R | 35 ++-
R/Ns.r | 2 +-
R/addCov.Lexis.R | 113 +++++++++
R/ci.lin.R | 53 +++-
R/clogistic.R | 5 +-
R/crr.Lexis.r | 18 +-
R/erl.R | 21 +-
R/lexis.R | 3 +
R/lls.R | 19 +-
R/mcutLexis.R | 151 ++++++++++++
R/plotCIF.R | 54 +++++
R/stackedCIF.R | 72 ++++++
R/xgrep.R | 3 +
R/zzz.R | 4 -
build/vignette.rds | Bin 276 -> 309 bytes
data/DMepi.rda | Bin 0 -> 76086 bytes
inst/doc/Follow-up.pdf | Bin 206532 -> 206908 bytes
inst/doc/index.html | 3 +
inst/doc/simLexis.R | 40 +--
inst/doc/simLexis.pdf | Bin 352900 -> 358803 bytes
inst/doc/simLexis.rnw | 26 +-
inst/doc/yll.R | 206 ++++++++++++++++
inst/doc/yll.pdf | Bin 0 -> 291815 bytes
inst/doc/yll.rnw | 641 +++++++++++++++++++++++++++++++++++++++++++++++++
man/DMepi.Rd | 41 ++++
man/Icens.Rd | 2 +-
man/LCa.fit.Rd | 38 +--
man/Lexis.Rd | 2 +
man/Ns.Rd | 2 +-
man/addCov.Lexis.Rd | 133 ++++++++++
man/cal.yr.Rd | 2 +-
man/ci.cum.Rd | 21 +-
man/ci.lin.Rd | 42 +++-
man/cutLexis.Rd | 11 +-
man/erl.Rd | 30 +--
man/fit.mult.Rd | 2 +-
man/foreign.Lexis.Rd | 2 +-
man/gen.exp.Rd | 22 +-
man/lgrep.Rd | 39 +++
man/mcutLexis.Rd | 96 ++++++++
man/ncut.Rd | 2 +-
man/nice.Rd | 2 +-
man/plotCIF.Rd | 122 ++++++++++
man/plotEst.Rd | 2 +-
man/splitLexis.Rd | 6 +-
man/stack.Lexis.Rd | 2 +-
man/summary.Lexis.rd | 4 +-
man/transform.Lexis.Rd | 12 +-
src/chinv2.c | 4 +-
src/cholesky2.c | 4 +-
src/chsolve2.c | 4 +-
src/clogit.c | 6 +-
src/init.c | 18 ++
vignettes/index.html | 3 +
vignettes/simLexis.rnw | 26 +-
vignettes/yll.rnw | 641 +++++++++++++++++++++++++++++++++++++++++++++++++
60 files changed, 2807 insertions(+), 226 deletions(-)
diff --git a/CHANGES b/CHANGES
index c678b6d..d206de6 100644
--- a/CHANGES
+++ b/CHANGES
@@ -1,3 +1,88 @@
+Changes in 2.19
+
+o Typos in documentation fixed, LCa.fit
+
+o Bug in knot calculation in LCa.fit fixed. Meaningless models emerged
+ if explicit knots were supplied for cohort effects. Prior to 2.19
+ only supplying *number* of knots for effects would give meaningful
+ models if a cohort effect were included (with or without
+ age-interaction).
+
+WISH: gen.exp has now got a wrapper, genExp.Lexis, explicitly using the
+ Lexis structure.
+
+Changes in 2.18
+
+o addCov.Lexis was fundamentally flawed, re-written, argument names and
+ order changed too.
+
+o Documentation links between addCov.Lexis and gen.exp are introduced.
+
+Changes in 2.16
+
+o Function addCov.Lexis added. Allows addition of covariates (clinical
+ mesurements) taken at a particular time to be added to a Lexis
+ object.
+
+Changes in 2.15
+
+o typo in stackedCIF code corrected (caused a crash with ony one group)
+
+Changes in 2.14
+
+o plotCIF, stackedCIF plotting Nelson-Aalen-Johansen estimators of
+ cumulative risks added, courtesy Esa L��r�
+
+o Convenice wrappers for grep to select elements: fgrep, ngrep, lgrep1.
+
+o A bug in Ns has been fixed, thanks to Lars J Diaz (DK) and Stephen
+ Wade (AUS).
+
+o surv1, surv2, yll, erl: NAs in input rates are now changed to 0
+ (with a warning) instead of crashing the function.
+
+o documentation of ci.cum groomed.
+
+Changes in 2.12
+
+o New function ci.ratio to compute RR with CIs from independent
+ estimates of rates
+
+Changes in 2.11
+
+o Small errors in calculation of knots in the simLexis macro corrected
+
+o A severe bug in mcutLexis which caused omission of certain cuts has
+ been fixed.
+
+Changes in 2.10
+
+o Bug in lls() caused a crash when objects had funny names (such as
+ '[.Lexis'). Fixed.
+
+Changes in 2.9
+
+o Grooming of code and documentation for mcutLexis.
+
+Changes in 2.8
+
+o A function, mcutLexis, to cut at several different event times,
+ preserving intermediat event histories has been added.
+
+o Errors in the erl.Rd corrected: Description of the argument "immune"
+ was wrong, as were the description of the timepoints where rates
+ were supposed given.
+
+o Inaccuracies in the vignette for simLexis patched.
+
+o lls() now also lists the size of objects
+
+o For illustrative purposes the DMepi dataset has been included
+
+Changes in 2.7
+
+o '[.Lexis' redefined to comply with data.table as used from popEpi
+
Changes in 2.6
o Added function erl computing Expected Residual Lifetime in an
@@ -7,8 +92,8 @@ o Added function erl computing Expected Residual Lifetime in an
Changes in 2.5
o Argument "timeScales" added to summary.Lexis, printing names of
- timescales and which (if any) are defined as time since enty into a
- state.
+ timescales and which of them (if any) are defined as time since enty
+ into a state.
o rm.tr added; removes transitions from a Lexis object
diff --git a/DESCRIPTION b/DESCRIPTION
index e510afb..29518c9 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,16 +1,16 @@
Package: Epi
-Version: 2.7
-Date: 2016-10-04
+Version: 2.19
+Date: 2017-08-01
Title: A Package for Statistical Analysis in Epidemiology
Authors at R: c(person("Bendix", "Carstensen", role = c("aut", "cre"),
- email = "bxc at steno.dk"),
+ email = "b at bxc.dk"),
person("Martyn", "Plummer", role = "aut",
email = "plummerm at iarc.fr"),
person("Esa", "Laara", role = "ctb"),
person("Michael", "Hills", role = "ctb"))
Depends: R (>= 3.0.0), utils
Imports: cmprsk, etm, splines, MASS, survival, plyr, Matrix, numDeriv,
- data.table
+ data.table, zoo
Suggests: mstate, nlme, lme4
Description: Functions for demographic and epidemiological analysis in
the Lexis diagram, i.e. register and cohort follow-up data, in
@@ -19,16 +19,16 @@ Description: Functions for demographic and epidemiological analysis in
'mstate', 'etm' and 'cmprsk' packages.
Also contains functions for Age-Period-Cohort and Lee-Carter
modeling and a function for interval censored data and some useful
- functions for tabulation and plotting, as well some epidemiological
- datasets.
+ functions for tabulation and plotting, as well as a number of
+ epidemiological data sets.
License: GPL-2
URL: http://BendixCarstensen.com/Epi/
NeedsCompilation: yes
-Packaged: 2016-10-03 19:59:16 UTC; bendix
+Packaged: 2017-08-08 20:03:05 UTC; bendix
Author: Bendix Carstensen [aut, cre],
Martyn Plummer [aut],
Esa Laara [ctb],
Michael Hills [ctb]
-Maintainer: Bendix Carstensen <bxc at steno.dk>
+Maintainer: Bendix Carstensen <b at bxc.dk>
Repository: CRAN
-Date/Publication: 2016-10-04 14:35:49
+Date/Publication: 2017-08-09 04:39:25 UTC
diff --git a/MD5 b/MD5
index f116753..bd0069b 100644
--- a/MD5
+++ b/MD5
@@ -1,21 +1,22 @@
-203c582ea761fd50ca3bb7549f87f920 *CHANGES
-b812b1d33d655d933f3c57e067a687f9 *DESCRIPTION
-98ff4454ab2bc144d016e0cf6a00b8a0 *NAMESPACE
+9a824270cc5a762fea93b89f7be6af50 *CHANGES
+7ec191f7f9198e00e243799ac541a98a *DESCRIPTION
+3889ddbc9a5ac1ba4f04ec84d8dfc26e *NAMESPACE
fac3c7f01ab0bd6930cf3cbea20e6bcc *R/Aplot.R
a031997be7388601dd1df6a28fc55328 *R/Cplot.R
6a1c4f4cbbf509d1f3c40b9f0adb6399 *R/Icens.R
-ecdc4193570616dd8b42ca451941fb29 *R/LCa.fit.R
+6dcc1bf2ed458a8898d366791de02f7d *R/LCa.fit.R
d802fc697f55ab0ef5ccd45e51776e46 *R/Lexis.diagram.R
c6b992622881e60fb82ee4b8c3357c9d *R/Lexis.lines.R
05c38381cc48ab25b681e8615ac31ea8 *R/Life.lines.R
a25e4901b0899a6b00285713117e3f21 *R/N2Y.r
8743609003b32036a7a1d02fb4b8e17b *R/NArray.R
-9a46235d94dabfbf52b324083f573492 *R/Ns.r
+876bd6183d0fdb8a180e462715675e0c *R/Ns.r
e9f582b5dc17d07bf2df96de74441963 *R/Pplot.R
b1c0f04443499205645ad578a3dd8573 *R/ROC.R
0a7bbd37e75920b54fd0fbe00a86511d *R/Relevel.R
0e3e3d23b8ab79dcf3046828190462b2 *R/Termplot.R
28a0ec0c8584e22759b516ff69e5bd78 *R/Wald.R
+6f35204baee4552a741c7262ac33a05d *R/addCov.Lexis.R
b1f0000430d6e15f59fe949ed8f1beb1 *R/apc.LCa.R
847f61901b2c9123367be3ebef30abe6 *R/apc.fit.R
3c53b3d667196c204fab811d2231ccd5 *R/apc.frame.R
@@ -25,21 +26,21 @@ f4e2a039dda786e4cad0b2961e3dd599 *R/boxes.MS.R
44e5f888817c94b8a190eae4aa964a2c *R/cal.yr.R
6988e76eb64f004a8c5029c23f8603af *R/ccwc.R
d5171eb23f0b431a7cc89b69d03308a6 *R/ci.cum.R
-1e69aa59a2367037c3a054bfa1a2759e *R/ci.lin.R
+ab5fd1bbdc59ef536f103bf533a5f143 *R/ci.lin.R
6f80c7b2b2120068cd16c043fbc29db3 *R/ci.mat.R
886f87cbd08266fa281e2ed7b3b3c253 *R/ci.pd.R
7637da11fdfb3612df5d73a64ee044e4 *R/clear.R
-89f2fe2a79bba0565efc02fd4aaeca9a *R/clogistic.R
+e7e271d9c47316877e9726fdaaa9699f *R/clogistic.R
d365bc4739ebb104b0318c1c9ee226fd *R/contr.2nd.R
bcd4142d3bf3e4517397387fe43fdf95 *R/contr.cum.R
8064804ee2b1cfa0cea2b44a58106b0a *R/contr.diff.R
8ffc3eabb90b6b9b710109015779d853 *R/contr.orth.R
-84ac435b69e71e6b8844d267570c6d28 *R/crr.Lexis.r
+efb8f1a3b97b41244275534855707b8b *R/crr.Lexis.r
281380e913564ab522664f2255b2e8a8 *R/cutLexis.R
14270f90ed00f4f043286d50e2de9f26 *R/detrend.R
a86250cce524f9fa2215c0c0cac4c6e7 *R/effx.match.r
9b46438b1c6fdc926ac1a253a85096aa *R/effx.r
-fac08dcafcfddfe8e9857f2546b5f0ea *R/erl.R
+abd3fdeb059a214f253e2f3e4658a302 *R/erl.R
822e02862de39625e5ed9ef6f84bdf39 *R/expand.data.r
0466dba42d10766e0b9cd220ad39dcf4 *R/factorize.R
67ce5a28408f42c7b65084195897b995 *R/fit.add.r
@@ -49,11 +50,13 @@ fac08dcafcfddfe8e9857f2546b5f0ea *R/erl.R
23a68a19f8f8e70c5664eb5110b3e807 *R/foreign.R
da7b5d41fe7876ef6bab01839b613939 *R/ftrend.R
f0f84fc26294b1cb3016a9e23974c545 *R/gen.exp.R
-52ccaf61645561d9c050e5d5665fbdcd *R/lexis.R
-9c575049aa6c3841e7e92f563a6fc450 *R/lls.R
+06f8a1ac3afc69610b9edfebc638f17b *R/lexis.R
+f8eba91de17a8207c6303ffb7b34f0a7 *R/lls.R
+d593af0ab833c8655469a24891e9116f *R/mcutLexis.R
79d6e0b02a2c4e98fbbfcd936b21d377 *R/mh.R
9e48e47769874a34f63a2e2d4e6ffffc *R/ncut.r
2a44a66ac4fe7eccab191f94957265e2 *R/pctab.R
+ff45a8401b6862061389134357b610cc *R/plotCIF.R
bb1ce8ea65f95b7e1720c2686e4de799 *R/plotEst.R
57eecdd415b7d0f7dfad39825644add5 *R/plotevent.r
c0229c2819d977a8c557a22dea1dd5fa *R/print.Icens.r
@@ -64,14 +67,16 @@ cca51283f637927702fbe217367011eb *R/print.floated.R
804e11d9d564a860fce275258f40db3d *R/simLexis.R
358b6fa60093d8bbee969a76a88797c0 *R/splitLexis.R
e08261ab194f4aca8502ee249afd309e *R/stack.Lexis.R
+8796d2af48a79385ad22b03ed0f4059c *R/stackedCIF.R
fb143d532001643d3f74fbd2e114e976 *R/stattable.R
889282ea37a5a9d9bf7155d790cd9036 *R/summary.Icens.r
8ceccd0ffa9694670f38c34b65d5be10 *R/summary.Lexis.r
88983ba8314a37910993d469a982d800 *R/twoby2.R
-497daef89f1ea705cd2c755132bec69d *R/zzz.R
-19da589577365d1c7f87fa6460dfefdc *build/vignette.rds
+53904292629435d35ec42f11329b2c60 *R/xgrep.R
+1ffde7eb6c78d26bcb1e51a138e63485 *build/vignette.rds
b9d3112271b9545896c98ec6406f5f02 *data/B.dk.rda
666e336a1aefb9c4298abb83a81a3300 *data/DMconv.rda
+ea83ed46a5581ea8489646abd3d1c0ae *data/DMepi.rda
e200a06eee9305b2fa2c59af5982d62b *data/DMlate.rda
e8f6f529b77cb3a860b9d53f5f65c3bb *data/DMrand.rda
ba4f4ad111ab2f488d1373c1cfcb7398 *data/M.dk.rda
@@ -96,18 +101,22 @@ deb94384c7380312f9b33f75cf6d8559 *data/nickel.rda
1953ddd7d750051c5682fa08b11cb777 *data/thoro.rda
9bd39178a387f935acc54519109c9f4a *inst/CITATION
e61ded63f726562d6f3172838e4f7337 *inst/doc/Follow-up.R
-f6e45b6606e4c70ae545b651bc3ea59b *inst/doc/Follow-up.pdf
+d8bfa83259c458592bdf1e8cabca4247 *inst/doc/Follow-up.pdf
a6978ca7d45a578a4fb047fb92eeed73 *inst/doc/Follow-up.rnw
-fff6643f55321f70c23cb36b9ee587e0 *inst/doc/index.html
-a9ec600c3c54c2b2e0a252f599bda85a *inst/doc/simLexis.R
-6ccba2ccef7004dd6f2c129f0e80afce *inst/doc/simLexis.pdf
-c38c581cd2ac1acdcce9e60e56dc37ca *inst/doc/simLexis.rnw
+364611344a5e791ae93f506e961b5fcc *inst/doc/index.html
+8d64b0389082fb54478f8e9667934265 *inst/doc/simLexis.R
+441f8e7242cf6b1ef990ffd642ff56fa *inst/doc/simLexis.pdf
+37f8cd2365c025fbb1d348052c67f712 *inst/doc/simLexis.rnw
+3cf7875e75ee2ef80cf2eda579764d6f *inst/doc/yll.R
+b4a034b2c6be5a4e5b63d7f38007a142 *inst/doc/yll.pdf
+5b597c102643a5597a70ec91abc1899c *inst/doc/yll.rnw
317bb6ebeaec5e7376b609d7e382c166 *man/B.dk.Rd
3abc5c20e62c874d70dcb2688ad2cce5 *man/DMconv.Rd
+e9a988b029749169403ab0bf749d53da *man/DMepi.Rd
362200b221355a7bcbcd1c7c8df5243f *man/DMlate.Rd
-38f9ae52c47f70d01f46eb437fca0f16 *man/Icens.Rd
-39d7b1fc6a514f8a6fd7d5013893a4c5 *man/LCa.fit.Rd
-1ab9326a322e35e9855be902cc2a5da8 *man/Lexis.Rd
+c87e6a0d0a5c461ecf4df07c2c96f1c1 *man/Icens.Rd
+c20a70a6f6c4f815e9a7c90d26529cab *man/LCa.fit.Rd
+e90e2ce2e2bc9daf1283492261c8ebbd *man/Lexis.Rd
c97dabafce248e4d47c5b18b6f04efe0 *man/Lexis.diagram.Rd
2e71ca58213788c8a358708e776e2112 *man/Lexis.lines.Rd
b4b2ca144feacece204f86bc6640d51a *man/Life.lines.Rd
@@ -115,12 +124,13 @@ daf66f950320e836b8ce1d2bd37cdc41 *man/M.dk.Rd
6ae7494a426801b9944c37190c09b3ee *man/N.dk.Rd
0983e4031974eea1a30e67db90310c25 *man/N2Y.Rd
424aecfdd664b3123ad4815880aa3761 *man/NArray.Rd
-5d1d5dcea53632bf7ea6f18adb174014 *man/Ns.Rd
+03f9885a6c523d7a221cf0a2272628fc *man/Ns.Rd
335b75739b8a25293690da804cd5e522 *man/ROC.Rd
ca827755eae8ee8286860a715cbd1744 *man/Relevel.Rd
bb2162d557dfb7685edac64ac3213b6c *man/S.typh.Rd
9351888a37017068f2c9ac8b843c846d *man/Termplot.Rd
b4879df831de32e9e298b7b3d90389a9 *man/Y.dk.Rd
+b813883bd7caca8d97ee45f456932eb4 *man/addCov.Lexis.Rd
62c36534ab2b6043f68306dd33407891 *man/apc.LCa.Rd
d5976837ba976d5969aaaedb98b215eb *man/apc.fit.Rd
8e42054ad19c34f2bf77210cfafc560e *man/apc.frame.Rd
@@ -130,71 +140,76 @@ d74a4f6e66ce91c73cd2817bb690446a *man/births.Rd
8dc11840302a388e4f26892335d5d888 *man/blcaIT.Rd
6b9f6bc7ca862808b89f6439914b8a0a *man/boxes.MS.Rd
18d942d375c5d5f6899291bf9cd1fe11 *man/brv.Rd
-a1a09c9591d54235cc026e44ca4ddad8 *man/cal.yr.Rd
+a2bedd44df3606b6585fd1da897de7ec *man/cal.yr.Rd
effba2c17ade00f0afd10184a073d34c *man/ccwc.Rd
-c47967c97d376b54706caacce489a9a2 *man/ci.cum.Rd
-fd00c162c3b545d86546615dfc24750f *man/ci.lin.Rd
+90337e7772c67e4e140e21165c30f079 *man/ci.cum.Rd
+5fca264afe8e1bfd25f93d365d9156fa *man/ci.lin.Rd
e7477402aa4746d653a0c72b2239d185 *man/ci.pd.Rd
128c0b1d00a69e41d131a6f74a1bae23 *man/clogistic.Rd
cf5e5a6e18e07b7bb8d70dc0417e6166 *man/contr.cum.Rd
24343adc886e90c75f7af178201bd6cd *man/crr.Lexis.rd
-b6ece2eb3b57ea902b737e87db9d9fb3 *man/cutLexis.Rd
+7d42e6a788c81361f03f03af15a2187c *man/cutLexis.Rd
a1ab8c14b721360ba0cb17da598bfa84 *man/detrend.Rd
386d0a4aa3515cf754ef458e1f32089c *man/diet.Rd
64dd6599bb830e98f4717b1a34ee12ab *man/effx.Rd
3c3cc3b1ecabece6b7709c1b147bab4e *man/effx.match.Rd
-139ed9b0278b94d417e9464472c3ef5c *man/erl.Rd
+a4922d13042b2d343868518a072e64d3 *man/erl.Rd
0dc62e726cdff6694774fab66595ee81 *man/ewrates.Rd
eb0c50c4a2572fd514c475358165fa57 *man/expand.data.rd
344adb137962cfa8fc7c2cf2858de3a3 *man/fit.add.Rd
3b44fed4ae70541c965c4cfbed966a68 *man/fit.baseline.rd
-a318bd705eddc9d721e59c95a6932ef0 *man/fit.mult.Rd
+b92589537d0e1c4abebdf1aad6abf88b *man/fit.mult.Rd
bca8b8c3452180e63d4d22885b04e3ce *man/float.Rd
-c96cdd9ce1f87335dd8eb454341f44ba *man/foreign.Lexis.Rd
+3afd2a7cb9c812dfa783cb916828b7cf *man/foreign.Lexis.Rd
e5a57a203cb91e9680dedc9023779e59 *man/ftrend.Rd
-9a1a04510b887f39bab1c2db4b6d2cfd *man/gen.exp.Rd
+e7e33d1ae13dec7a5cf3c3f96e1c085e *man/gen.exp.Rd
cbbc2a23f902d83ba2527b27b5ed3adb *man/gmortDK.Rd
d9f1d31e109d6058a2160ca0aa49bf81 *man/hivDK.Rd
9c3b059921d728298ccff73298b4dd31 *man/lep.Rd
+470bf279897b9c9ee2600126127d8ada *man/lgrep.Rd
ec0e27968d4ee35f390161e49abd903c *man/lines.apc.Rd
78125ba49d54905fb727765dc8997796 *man/lls.Rd
661ec298ffcd0f39cf55c5d033e04fd1 *man/lungDK.Rd
+cb1f50d16fe21259603273d2004d14b3 *man/mcutLexis.Rd
3c462c4ce7ca1ce4b8c6fef50b3c7623 *man/merge.Lexis.Rd
4c26c9edf5d64b0af3e63ad791b2a280 *man/merge.data.frame.Rd
34c530c36a3a1a7447ff8eea58f42084 *man/mh.Rd
2ecfec2d86e80b9bbcc2c3a253348edd *man/mortDK.Rd
-0dbb7c541b3ebd70a1cd104dbf0df25d *man/ncut.Rd
-dec4748ef3ff7c520df1075ec8296b25 *man/nice.Rd
+9a531412399076898e3e39c1c20eb69a *man/ncut.Rd
+880466737da2c1037259c3f712a7e9fe *man/nice.Rd
05bd25340716d0805cbc7eb433896206 *man/nickel.Rd
0e0ff2e0894523b3878213d03afcc8cf *man/occup.Rd
fb7ba7ea4cc2a9a67529348a4d971a4f *man/pc.lines.Rd
86cac5433077a16960bca5d5865586a4 *man/pctab.Rd
50c0516a536004d7f4ce91ef1bbfd8c9 *man/plot.Lexis.Rd
d30b27f4c3ef397de514614ba3e6ba3e *man/plot.apc.Rd
-43909fc8f0b297bbc879237717e18d4b *man/plotEst.Rd
+41c46ab16d54e331e163a455e67f9d20 *man/plotCIF.Rd
+5e3bda25091001df87ded934e0fd929f *man/plotEst.Rd
b0047a7bab6a20b153485db9bd338394 *man/plotevent.rd
feb2228ef38635fcbf9083253c2f6d21 *man/projection.ip.rd
a04f5841de6a7a9cf9a21e5caa8e54a5 *man/rateplot.Rd
78a9dbec1798f1b4e5735c71756f0ec6 *man/rbind.Lexis.Rd
c5df7f725aed8ec3520ebd9735deeccb *man/rm.tr.Rd
9fee82808de1916808de45242ebafd78 *man/simLexis.Rd
-a76a901889d88cb7e0301de3c8d07680 *man/splitLexis.Rd
-c146b743be57cba543822ae2f613b08b *man/stack.Lexis.Rd
+a22aa8901f1f717b6d3a7ed086bc8e02 *man/splitLexis.Rd
+00fd0080aeb0551b2af69b6d918265e2 *man/stack.Lexis.Rd
eaad8ae22bd68a18750a94352f65f9c3 *man/start.Lexis.Rd
4f43436fc863747dad4fcdd566948039 *man/stattable.Rd
73cf4597fcb4a9d019bda47431b095be *man/stattable.funs.Rd
f10fc89dadf571e07f3b13923660cff6 *man/subset.Lexis.Rd
-2639c132d73f56d51bca04d8d5e22b66 *man/summary.Lexis.rd
+75b79b094aa4e94218df58cd5bec86d9 *man/summary.Lexis.rd
14711d23ae0b89d9d22b9f3b07211093 *man/testisDK.Rd
e3446824b09ea4ad8e8a0ddfc1987eef *man/thoro.Rd
6ea9cf64475187fa2f402275c40573a7 *man/time.band.Rd
2edb93539fa7b9f47a3083b624443f9d *man/time.scales.Rd
-b827f35f4206cd67af318cc526d98f2f *man/transform.Lexis.Rd
+f80e2b3fe7e1d9df9c40dcd577a3f499 *man/transform.Lexis.Rd
02fc311541fae16a5d3a580a04925050 *man/twoby2.Rd
-3b649b94df7f0c0863bb4ef0a7129855 *src/chinv2.c
-c21af4e77e577fae093bb3dfc6cbe9a9 *src/cholesky2.c
-33be867717f019d03821e58b662a94c4 *src/chsolve2.c
-ce9e73f4028463cf8b84a7783f282ac3 *src/clogit.c
+bcc0556f8f0a158637592731928cf835 *src/chinv2.c
+3e3a51847b5767274f0945350e947726 *src/cholesky2.c
+09bf10fed80e805e7e58db9794e5edcc *src/chsolve2.c
+872ef5e7ef0c1941c17d7dd612f6b1c3 *src/clogit.c
+d77bd820f6ce7dc7883c9e2813a44c89 *src/init.c
a6978ca7d45a578a4fb047fb92eeed73 *vignettes/Follow-up.rnw
-fff6643f55321f70c23cb36b9ee587e0 *vignettes/index.html
-c38c581cd2ac1acdcce9e60e56dc37ca *vignettes/simLexis.rnw
+364611344a5e791ae93f506e961b5fcc *vignettes/index.html
+37f8cd2365c025fbb1d348052c67f712 *vignettes/simLexis.rnw
+5b597c102643a5597a70ec91abc1899c *vignettes/yll.rnw
diff --git a/NAMESPACE b/NAMESPACE
index 6c9a086..6a8f70b 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,3 +1,4 @@
+useDynLib(Epi, .registration=TRUE)
export(
as.Date.cal.yr,
apc.frame,
@@ -18,6 +19,7 @@ lines.apc,
ci.exp,
ci.cum,
ci.pred,
+ ci.ratio,
ci.mat,
lls,
clear,
@@ -39,6 +41,11 @@ lines.apc,
fit.add,
fit.mult,
ftrend,
+ fgrep,
+ ngrep,
+ lgrep,
+ plotCIF,
+ stackedCIF,
show.apc.LCa,
apc.LCa,
LCa.fit,
@@ -67,6 +74,7 @@ lines.apc,
Relevel.Lexis,
factorize.Lexis,
cutLexis,
+ mcutLexis,
countLexis,
stack.Lexis,
tmat.Lexis,
@@ -77,6 +85,7 @@ lines.apc,
etm.Lexis,
crr.Lexis,
simLexis,
+ addCov.Lexis,
subset.stacked.Lexis,
transform.stacked.Lexis,
plot.pState,
@@ -136,6 +145,7 @@ importFrom( MASS, mvrnorm )
importFrom( survival, clogit )
importFrom( numDeriv, hessian )
importFrom( Matrix, nearPD )
+importFrom( zoo, na.locf )
importFrom("grDevices", "gray", "rainbow")
importFrom("graphics", "abline", "arrows", "axis", "box", "layout",
"lines", "locator", "matlines", "matplot", "matpoints",
@@ -148,7 +158,7 @@ importFrom("stats", ".getXlevels", "addmargins", "anova", "approx",
"is.empty.model", "median", "model.extract", "model.matrix",
"model.offset", "model.response", "nlm", "pchisq", "pnorm",
"poisson", "predict", "qnorm", "qt", "quantile", "runif",
- "termplot", "update", "vcov", "weighted.mean", "xtabs")
+ "termplot", "update", "vcov", "weighted.mean", "xtabs" )
importFrom("utils", "flush.console", "str")
# register S3 methods
@@ -177,6 +187,7 @@ S3method( Relevel, Lexis)
S3method( Relevel, factor)
S3method( Relevel, default)
S3method(factorize, Lexis)
+S3method( addCov, Lexis)
S3method( stack, Lexis)
S3method( tmat, Lexis)
S3method( boxes, Lexis)
diff --git a/R/LCa.fit.R b/R/LCa.fit.R
index b573043..b842b99 100644
--- a/R/LCa.fit.R
+++ b/R/LCa.fit.R
@@ -23,10 +23,10 @@ function( data, A, P, D, Y,
if( !(model %in% c("APa","ACa","APaC","APCa","APaCa")) )
stop( '"model" must be one of "APa","ACa","APaC","APCa","APaCa", but is', model,'\n' )
-# Which main effects and interactiosn are in the model
+# Which main effects and interactions are in the model
intP <- as.logical(length(grep("Pa",model)))
-mainP <- as.logical(length(grep("P" ,model))) # Also includes the age-period interaction
intC <- as.logical(length(grep("Ca",model)))
+mainP <- as.logical(length(grep("P" ,model))) # Also includes the age-period interaction
mainC <- as.logical(length(grep("C" ,model))) # Also includes the age-cohort product
# if a dataframe is supplied, fish out data and put in the function's environment
@@ -43,11 +43,10 @@ if( !missing(data) )
D <- data$D
Y <- data$Y
} else { # if single vectors supplied, check they are all there
- nm <- logical(4)
- nm[1] <- missing(A)
- nm[2] <- missing(P)
- nm[3] <- missing(D)
- nm[4] <- missing(Y)
+ nm <- c(missing(A),
+ missing(P),
+ missing(D),
+ missing(Y))
if (any(nm))
stop("Variable", if (sum(nm) > 1)
"s", paste(c(" A", " P", " D", " Y")[nm], collapse = ","),
@@ -68,26 +67,26 @@ if( is.list(npar) ) {
# Check if names is a named list
if( is.null(names(npar)) ) stop( "If npar= is a list, it must be a *named* list.\n" )
a.kn <- if( length(npar$a )>1 ) npar$a else quantile( rep( A,D), probs=eqqnt(npar$a ) )
- pi.kn <- if( length(npar$pi)>1 ) npar$pi else quantile( rep( A,D), probs=eqqnt(npar$pi) )
p.kn <- if( length(npar$p )>1 ) npar$p else quantile( rep(P ,D), probs=eqqnt(npar$p ) )
+ c.kn <- if( length(npar$c )>1 ) npar$c else quantile( rep(P-A,D), probs=eqqnt(npar$c ) )
+ pi.kn <- if( length(npar$pi)>1 ) npar$pi else quantile( rep( A,D), probs=eqqnt(npar$pi) )
ci.kn <- if( length(npar$ci)>1 ) npar$ci else quantile( rep( A,D), probs=eqqnt(npar$ci) )
- c.kn <- if( length(npar$b )>1 ) npar$c else quantile( rep(P-A,D), probs=eqqnt(npar$c ) )
}
else { # if npar is too short fill it up
npar <- rep( npar, 5 )[1:5]
# if not named, name it and notify
if( is.null(names(npar)) ) names(npar) <- c("a","p","c","pi","ci")
a.kn <- quantile( rep( A,D), probs=eqqnt(npar["a"] ) )
- pi.kn <- quantile( rep( A,D), probs=eqqnt(npar["pi"]) )
p.kn <- quantile( rep(P ,D), probs=eqqnt(npar["p"] ) )
- ci.kn <- quantile( rep( A,D), probs=eqqnt(npar["ci"]) )
c.kn <- quantile( rep(P-A,D), probs=eqqnt(npar["c"] ) )
+ pi.kn <- quantile( rep( A,D), probs=eqqnt(npar["pi"]) )
+ ci.kn <- quantile( rep( A,D), probs=eqqnt(npar["ci"]) )
}
# Reference points
if( missing( p.ref) ) p.ref <- median( rep(P ,D) )
-if( missing(pi.ref) ) pi.ref <- median( rep( A,D) )
if( missing( c.ref) ) c.ref <- median( rep(P-A,D) )
+if( missing(pi.ref) ) pi.ref <- median( rep( A,D) )
if( missing(ci.ref) ) ci.ref <- median( rep( A,D) )
############################################################################
@@ -482,8 +481,8 @@ if( mt$intC ) {
# First fitted values from mod.at
# Note that the model object mod.at always has the same number of
-# parameters, for some of the model eitehr period or cohort parameters
-# are 0, so not used.
+# parameters, for some of the models either period or cohort parameters
+# are 0, hence not used.
pr0 <- ci.exp( object$mod.at, alpha=alpha, ctr.mat=cbind(Ma,Mp*pi,Mc*ci) )
# Then fitted values from mod.b
@@ -495,10 +494,10 @@ pr0 <- cbind( pr0, exp( lp.b %*% ci.mat(alpha=alpha) ) )
# label the estimates
colnames( pr0 )[c(1,4)] <- c("at|b Est.","b|at Est.")
-# This gives confidence intervals based on the conditional models,
-# so if we want proper intervals we should simulate instead, using the
-# posterior distribuion of all parameters, albeit under the slightly
-# fishy assumption that the joint posterior is normal...
+# The doings above gives confidence intervals based on the conditional
+# models, so if we want proper intervals we should simulate instead,
+# using the posterior distribuion of all parameters, albeit under the
+# slightly fishy assumption that the joint posterior is normal...
if( sim ) # also renders TRUE if sim is numerical (and not 0)
{
if( is.logical(sim) & sim ) sim <- 1000
diff --git a/R/Ns.r b/R/Ns.r
index 89cce6a..7ba7a35 100644
--- a/R/Ns.r
+++ b/R/Ns.r
@@ -129,7 +129,7 @@ Ns <- function( x, ref = NULL,
}
if( detrend & any(!is.na(fixsl)) ) {
warning( "detrend= specified, hence fixsl argument is ignored")
- fixsl=c(NA,NA)
+ fixsl=c(FALSE,FALSE)
}
## Here is the specification of the spline basis
## df= specified
diff --git a/R/addCov.Lexis.R b/R/addCov.Lexis.R
new file mode 100644
index 0000000..88eb34e
--- /dev/null
+++ b/R/addCov.Lexis.R
@@ -0,0 +1,113 @@
+# The addCov method
+addCov <- function (x, ...) UseMethod("addCov")
+
+addCov.default <-
+addCov.Lexis <-
+function( Lx,
+ clin,
+ timescale = 1,
+ exnam,
+ tfc = "tfc",
+ addScales = FALSE )
+{
+# Function to add clinically measured covariates to a Lexis object
+
+# The point is to cut the Lexis diagram at the examination dates
+# and subsequently add the clinical records
+# ...but first the usual cheking of paraphernalia
+
+if( !inherits(Lx ,"Lexis") ) stop( "Lx must be a Lexis object.\n" )
+if( inherits(clin,"Lexis") ) stop( "clin cannot be a Lexis object.\n" )
+
+# Is the timescale argument a timescale in Lx and is it a variable in clin?
+ts <- if( is.numeric(timescale) ) timeScales( Lx )[timescale] else timescale
+if( !( ts %in% timeScales(Lx) ) )
+ stop( "timescale argument (", ts, ") must be one of timescales in in the Lexis object ",
+ deparse(substitute(Lx)),":", timeScales(Lx), ".\n" )
+
+clin.nam <- deparse(substitute(clin))
+if( !( ts %in% names(clin) & "lex.id" %in% names(clin) ) )
+ stop( "'lex.id' and timescale '", ts, "' must be a variables in the clin object ",
+ clin.nam, "\n" )
+
+# variables to merge by
+mvar <- c( "lex.id", ts )
+
+# order clin to get the possible construction of examination names ok
+clin <- clin[order(clin$lex.id,clin[,ts]),]
+
+# check that examination dates are unique within persons
+if( any( dd <- duplicated(clin[,c("lex.id",ts)]) ) )
+ {
+ warning( "Examination dates must be unique within persons\n",
+ sum(dd), " records with duplicate times from clin object ", clin.nam,
+ " excluded.")
+ clin <- clin[!dd,]
+ }
+
+# the variable holding the name of the examination
+if( missing(exnam) ) exnam <- "exnam"
+# and if it is not there, construct it
+if( !(exnam %in% names(clin)) )
+ clin[,exnam] <- paste( "ex",
+ ave( clin$lex.id,
+ clin$lex.id,
+ FUN = function(x) cumsum(x/x) ),
+ sep="." )
+
+# Add copy of the time of examination to be carried forward
+clin[,tfc] <- clin[,ts]
+
+# Clinical variables to be merged in --- note we take examination date
+# and name as a cinical variable too
+cvar <- setdiff( names(clin), mvar )
+
+# A data frame of cutting times
+cfr <- data.frame( lex.id = clin$lex.id,
+ cut = clin[,ts],
+ new.state = clin[,exnam] )
+
+# Now cut Lx --- this is really inefficient
+mc <- Lx
+for( st in levels(cfr$new.state) )
+mc <- cutLexis( mc,
+ cut = cfr[cfr$new.state==st,],
+ timescale = ts,
+ precursor.states = NULL )
+
+# Merge in the old states
+mx <- Lx[,mvar]
+mx$org.Cst <- Lx$lex.Cst
+mx$org.Xst <- Lx$lex.Xst
+mx <- merge( mx, mc, by = mvar, all.y = TRUE, sort = TRUE )
+
+# Complete the state variables
+( wh <- which(is.na(mx$org.Cst)) )
+mx$org.Cst[wh] <- na.locf( mx$org.Cst, nx.rm=FALSE )[wh]
+mx$org.Xst[wh] <- na.locf( mx$org.Xst, nx.rm=FALSE )[wh]
+# but - oops - the earlier lex.Xst should be as lex.Cst
+mx$org.Xst[wh-1] <- mx$org.Cst[wh-1]
+# overwrite the useless ones and get rid of the intermediate
+mx$lex.Cst <- mx$org.Cst
+mx$lex.Xst <- mx$org.Xst
+wh.rm <- match( c("org.Cst","org.Xst"), names(mx) )
+mx <- mx[,-wh.rm]
+
+# Merge in the clinical variables
+mx <- merge( mx, clin, by=mvar, all.x=TRUE, sort=TRUE )
+
+# And carry them forward within each lex.id
+# locf within each person (should be easier in data.table)
+locf.df <- function( df ) as.data.frame( lapply( df, na.locf, na.rm=FALSE ) )
+# ave does not like character variables so we convert to factors
+wh <- which( sapply( mx[,cvar], is.character ) )
+for( j in wh ) mx[,cvar[j]] <- factor( mx[,cvar[j]] )
+# then we can carry forward
+mx[,cvar] <- ave( mx[,cvar], mx$lex.id, FUN=locf.df )
+
+# Finally update the time from clinical measurement
+mx[,tfc] <- mx[,ts] - mx[,tfc]
+
+# Done!
+mx
+}
diff --git a/R/ci.lin.R b/R/ci.lin.R
index b9664a2..d3a03e6 100644
--- a/R/ci.lin.R
+++ b/R/ci.lin.R
@@ -211,21 +211,64 @@ else
ci.lin( ..., Exp=FALSE )[,if(pval) c(1,5,6,4) else c(1,5,6),drop=FALSE]
}
-# Wrapper for predict.glm
+# Wrapper for predict.glm to give estimates and confidnece intervals
ci.pred <-
function( obj, newdata,
Exp = NULL,
- alpha = 0.05,
- df = Inf )
+ alpha = 0.05 )
{
if( !inherits( obj, "glm" ) ) stop("Not usable for non-glm objects")
-# get prediction and se on the link scale
+# get the prediction and se on the link scale
zz <- predict( obj, newdata=newdata, se.fit=TRUE, type="link" )
# compute ci on link scale
-zz <- cbind( zz$fit, zz$se.fit ) %*% ci.mat( alpha=alpha, df=df )
+zz <- cbind( zz$fit, zz$se.fit ) %*% ci.mat( alpha=alpha )
# transform as requested
if( missing(Exp) ) { return( obj$family$linkinv(zz) )
} else { if( Exp ) { return( exp(zz) )
} else if( !Exp ) return( zz )
}
}
+
+ci.ratio <-
+function( r1, r2,
+ se1 = NULL, # standard error of rt1
+ se2 = NULL, # standard error of rt2
+ log.tr = !is.null(se1) & !is.null(se2), # is this log-rates?
+ alpha = 0.05,
+ pval = FALSE )
+{
+# Function to calculate RR with CIs from independent rates with CIs;
+# r1 and r2 are assumed to be vectors or 2 or 3-column matrices with
+# rate, lower and upper confidence limits repectively.
+
+if( is.matrix(r1) & !is.null(se1) ) warning("r1 is matrix, se1 is ignored")
+if( is.matrix(r2) & !is.null(se2) ) warning("r2 is matrix, se2 is ignored")
+
+# if supplied as 1-column matrix chnage to vector
+if( is.matrix(r1) ) if( ncol(r1)==1 ) r1 <- as.vector( r1 )
+if( is.matrix(r2) ) if( ncol(r2)==1 ) r2 <- as.vector( r2 )
+
+# move to log scale
+if( !log.tr ) {
+ r1 <- log( r1 )
+ r2 <- log( r2 )
+ }
+
+# how wide are the condidence intervals
+if( is.matrix(r1) ) if( ncol(r1)>1 ) rg1 <- t( apply(r1,1,range) )
+if( is.matrix(r2) ) if( ncol(r2)>1 ) rg2 <- t( apply(r2,1,range) )
+
+# get the estimates on the log-scale
+R1 <- if( is.matrix(r1) ) apply( rg1, 1, mean ) else r1
+R2 <- if( is.matrix(r2) ) apply( rg2, 1, mean ) else r2
+if( is.null(se1) ) se1 <- apply( rg1, 1, diff ) / (2*qnorm(1-alpha/2))
+if( is.null(se2) ) se2 <- apply( rg2, 1, diff ) / (2*qnorm(1-alpha/2))
+
+# compute the RR and the c.i. and optionally the p-value
+ lrr <- R1 - R2
+slrr <- sqrt( se1^2 + se2^2 )
+ rr <- cbind(lrr,slrr) %*% ci.mat(alpha=alpha)
+if( !log.tr ) rr <- exp( rr )
+if( pval ) return( cbind( rr, 1-pchisq( (lrr/slrr)^2, 1 ) ) )
+ else return( rr )
+}
diff --git a/R/clogistic.R b/R/clogistic.R
index b4eeb97..72b5ce7 100644
--- a/R/clogistic.R
+++ b/R/clogistic.R
@@ -72,9 +72,8 @@ fitClogit <- function(X, y, offset, strata, init, iter.max, eps, toler.chol)
stop("There are no informative observations")
}
- ans <- .Call("clogit", Xsplit, ysplit, osplit, as.double(init),
- as.integer(iter.max), as.double(eps), as.double(toler.chol),
- PACKAGE="Epi")
+ ans <- .Call(C_clogit, Xsplit, ysplit, osplit, as.double(init),
+ as.integer(iter.max), as.double(eps), as.double(toler.chol))
ans$informative <- info
return(ans)
}
diff --git a/R/crr.Lexis.r b/R/crr.Lexis.r
index 5778c6c..8e39b79 100644
--- a/R/crr.Lexis.r
+++ b/R/crr.Lexis.r
@@ -15,17 +15,17 @@ mod[2] <- NULL
cv <- model.matrix( mod, data=obj )[,-1]
# Then do it
M <- crr( ftime = obj$lex.dur,
- fstatus = obj$lex.Xst,
- failcode = fc,
- cencode = cn,
- cov1 = cv,
- ... )
-# A table of the no of transitions
+ fstatus = obj$lex.Xst,
+ failcode = fc,
+ cencode = cn,
+ cov1 = cv,
+ ... )
+# A table of the no. of transitions
N <- with( obj, table( Relevel( lex.Xst, c(fc,cc,cn) ) ) )
names( N ) <- paste( rep( c("Event:"," comp:"," cens:"),
c(length(fc),length(cc),length(cn)) ),
names(N) )
-# add model an table to the resulting object
+# add model and table to the resulting object
M <- c( M, list( model.Lexis = md,
transitions = cbind(N) ) )
# remember the class attribute (lost by doing "c")
@@ -36,5 +36,5 @@ cat( "crr analysis of event", paste('"',fc,'"',sep=''),
"\n versus", paste('"',cc,'"',sep=''),
"\n with", paste('"',cn,'"',sep=''),
"as censoring.\n" )
- M
- }
+M
+}
diff --git a/R/erl.R b/R/erl.R
index 05d7914..a6cd1d9 100644
--- a/R/erl.R
+++ b/R/erl.R
@@ -25,6 +25,12 @@ sum( ( x[-length(x)] + x[-1] ) / 2, na.rm=TRUE )
# Check sensibility
if( !immune & is.null(lam) ) stop( "'lam' is required when immune=FALSE\n" )
+# Replace NAs with 0s
+if( !is.null(lam) ) {
+if( any(is.na(lam)) ){ lam[is.na(lam)] <- 0 ; warning("NAs in agument 'lam' set to 0") } }
+if( any(is.na(muD)) ){ muD[is.na(muD)] <- 0 ; warning("NAs in agument 'muD' set to 0") }
+if( any(is.na(muW)) ){ muW[is.na(muW)] <- 0 ; warning("NAs in agument 'muW' set to 0") }
+
# Survival functions
sD <- surv1( int=int, muD, age.in = age.in, A = A )
if( immune ) sW <- surv1( int=int, muW, age.in = age.in, A = A )
@@ -75,6 +81,9 @@ function( int, mu, age.in=0, A=NULL )
# int and mu should be in compatible units that is T and T^-1 for
# some unit T (months, years, ...)
+# impute 0s for NAs
+if( any(is.na(mu)) ){ mu[is.na(mu)] <- 0 ; warning("NAs in agument 'mu' set to 0") }
+
# age-class boundaries
age <- 0:length(mu)*int + age.in
@@ -131,6 +140,12 @@ if( length(muW) != length(muD) |
", length(muD)=", length(muD),
", length(lam)=", length(lam) )
+# Replace NAs with 0s
+if( !is.null(lam) ) {
+if( any(is.na(lam)) ){ lam[is.na(lam)] <- 0 ; warning("NAs in agument 'lam' set to 0") } }
+if( any(is.na(muD)) ){ muD[is.na(muD)] <- 0 ; warning("NAs in agument 'muD' set to 0") }
+if( any(is.na(muW)) ){ muW[is.na(muW)] <- 0 ; warning("NAs in agument 'muW' set to 0") }
+
# First the workhorse that computes the survival function for a
# person in Well assuming that the mortality rate from this state is
# muW, disease incidence is in lam, and mortality in the diseased
@@ -139,8 +154,8 @@ if( length(muW) != length(muD) |
wsurv2 <-
function( int, muW, muD, lam, age.in=0, A=0 )
{
-# age-class boundaries - note one longer that rate vectors refers to
-# boundaries of intervals not midpoints
+# age-class boundaries - note one longer than rate vectors as it
+# refers to boundaries of intervals not midpoints
age <- 0:length(muW)*int + age.in
# cumulative rates at the boundaries, given survival to A
@@ -181,7 +196,7 @@ if( !is.null(A) )
for( j in 1:length(A) )
{
surv <- cbind( surv,
- wsurv2( int, muW, muD, lam, age.in=age.in, A=A[j] )[,2] )
+ wsurv2( int, muW, muD, lam, age.in=age.in, A=A[j] )[,2] )
}
}
Al <- A
diff --git a/R/lexis.R b/R/lexis.R
index d7c48ee..12d5b47 100644
--- a/R/lexis.R
+++ b/R/lexis.R
@@ -770,3 +770,6 @@ transform.Lexis <- function(`_data`, ... )
attributes(y) <- save.at
y
}
+
+# order.Lexis <- function( x ) order( x$lex.id, lex[,timeScales(x)[1]] )
+# sort.Lexis <- function( x, ... ) x[order.Lexis(x,...),]
diff --git a/R/lls.R b/R/lls.R
index 68acdd3..ac93f49 100644
--- a/R/lls.R
+++ b/R/lls.R
@@ -11,20 +11,25 @@ if( length(lll) > 0 )
{
obj.mode <-
obj.clas <-
+obj.dimx <-
obj.size <- character(0)
# Then find mode, class, name and dimension of them and return it
for(i in 1:length(lll))
{
-obj.mode[i] <- eval( parse(text = paste( "mode(", lll[i], ")")))
-obj.clas[i] <- paste( eval( parse(text = paste("class(", lll[i], ")"))), collapse=" " )
-obj.size[i] <- paste( eval( parse(text = paste( "dimx(", lll[i], ")"))), collapse=" " )
+obj.mode[i] <- eval( parse(text = paste( "mode(`", lll[i], "`)",sep="")))
+obj.clas[i] <- paste( eval( parse(text = paste( "class(`", lll[i], "`)",sep=""))), collapse=" " )
+obj.dimx[i] <- paste( eval( parse(text = paste( "dimx(`", lll[i], "`)",sep=""))), collapse=" " )
+obj.size[i] <- formatC( eval( parse(text = paste("unclass(object.size(`", lll[i], "`))",sep="")))/2^10,
+ format="f", digits=1, big.mark=",", width=14, flag=" " )
}
dfr <-
-data.frame( name=lll,
- mode=obj.mode,
- class=obj.clas,
- size=obj.size,
+data.frame( name = lll,
+ mode = obj.mode,
+ class = obj.clas,
+ dim = obj.dimx,
+ sizeKbytes = obj.size,
stringsAsFactors=FALSE )
+names( dfr )[5] <- " size(Kb)"
print( invisible( dfr ), right=FALSE )
}
}
diff --git a/R/mcutLexis.R b/R/mcutLexis.R
new file mode 100644
index 0000000..9b4f111
--- /dev/null
+++ b/R/mcutLexis.R
@@ -0,0 +1,151 @@
+mcutLexis <-
+function( L0, # A Lexis object
+ timescale = 1, # the time scale referred to by L0[,wh]
+ wh, # indices/names of columns holding dates of state entries (events)
+ new.states = NULL, # Names of the event types (states)
+precursor.states = NULL,
+ seq.states = TRUE, # Should state names reflect ordering of events
+ new.scales = NULL, # Time-scales referring to time since
+ ties.resolve = FALSE # Are tied event times accepted?
+ )
+{
+### we rely on referring to the timescale and event time variables by name
+if( is.numeric(timescale) ) timescale <- timeScales(L0)[timescale]
+if( is.numeric(wh) ) wh <- names(L0)[wh]
+
+### don't be silly
+if( length(wh)==1 )
+# return( docut( L0, osv ) ) # old cutLexis should be absorbed here
+ stop( "mcutLexis not needed for one type of event - use cutLexis\n" )
+
+### states
+if( is.null(new.states) )
+ {
+ new.states <- wh
+ cat( "NOTE: Name of new states set to\n", new.states )
+ }
+if( length(wh) != length(new.states) )
+ stop( "wh and new.states must have same length, but lengths are",
+ "wh:", length(wh), "and new.states:", length(new.states), "\n" )
+
+### timescales
+# either all or none
+if( is.logical(new.scales) )
+ if( any( new.scales ) )
+ {
+ new.scales <- paste( "tf", new.states, sep="" )
+ cat( "NOTE: new.scales set to: ", new.scales, "\n" )
+ }
+if( is.character(new.scales) & length(new.scales) != length(wh) )
+ {
+ new.scales <- paste( "tf", new.states, sep="" )
+ warning( "new.scales not of same length as wh. Set to: ",
+ new.scales, "\n" )
+ }
+if( is.character(new.scales) & length(intersect(new.states,timeScales(L0))) )
+ stop( "Names of new time scales must be different from names of timescales:\n",
+ timeScales(L0) )
+
+### Tied transition times untied
+has.ties <- any( wh.tied <- apply( L0[,wh], 1,
+ function(x) any(diff(sort(x[!is.na(x)]))==0) ) )
+if( has.ties & is.logical(ties.resolve) & !ties.resolve )
+ stop( "Tied event times not allowed with ties.resolve=FALSE:\n",
+ "there were", length(wh.tied), "records with tied event times.")
+if( has.ties & is.logical(ties.resolve) & ties.resolve ) ties.resolve = 1/100
+if( has.ties & is.numeric(ties.resolve) )
+ {
+ rnd <- L0[wh.tied,wh]*0
+ rnd[,] <- runif(rnd,-1,1) * ties.resolve
+ L0[wh.tied,wh] <- L0[wh.tied,wh] + rnd
+ cat( "NOTE: ", length(wh.tied),
+ "records with tied events times resolved.\n",
+ "Results only reproducible if the seed for the random number generator is set.")
+ }
+# End of checks
+
+# The object to return initiated as NULL
+Lcut <- NULL
+
+# Utility function returning sequences of ocurrences as paste of numbers
+NAorder <-
+function (x)
+ {
+ oo <- order(x, na.last = T)
+ on <- (1:length(oo))[oo]
+ on[is.na(x[oo])] <- NA
+ paste(on[!is.na(on)], collapse = "-")
+ }
+
+# where do the different sequences of events actually occur in data
+L0$whseq <- apply( L0[,wh], 1, NAorder )
+
+# Loop through the actually occurring orders of event occurrences
+for( sq in unique(L0$whseq) )
+{
+# Persons with none of the events occurring transferred to result
+if( sq=="" ) Lcut <- rbind( Lcut, L0[L0$whseq=="",] )
+else {
+
+# Extract the subset of persons with a given sequence of events
+Ltmp <- L0[L0$whseq==sq,]
+
+# The numerical sequence of states (refer to the elements of wh)
+ost <- as.numeric( strsplit( sq, "-" )[[1]] )
+nxst <- ""
+prst <- precursor.states
+for( cs in ost )
+ {
+ nxst <- ifelse( cs==ost[1], new.states[cs],
+ paste( nxst, new.states[cs], sep="-" ) )
+ Ltmp <- cutLexis( Ltmp, cut = Ltmp[,wh[cs]],
+ timescale = timescale,
+ new.state = nxst,
+ precursor.states = prst )
+ # include the created state among the precursor states for next cut
+ prst <- c(prst,nxst)
+ } # end of for loop through events in this sequence (cs)
+
+# Attach it to the end of the Lexis object
+Lcut <- rbind( Lcut, Ltmp )
+ } # end of the else clause
+
+} # end of for loop through sequences (sq)
+
+# Do we want the sequences or just the unordered set of previous events
+if( !seq.states )
+ {
+ lvl <- levels( Lcut )
+ ulvl <- sapply( lapply( strsplit(lvl,"-"),
+ sort ),
+ paste,
+ collapse="+" )
+ levels( Lcut$lex.Cst ) <-
+ levels( Lcut$lex.Xst ) <- ulvl
+ }
+
+# Did we ask for timescales as time since events?
+if( !is.null(new.scales) )
+ {
+ # insert columns for the new time scales
+ Lcut <- Lcut[,c(rep("whseq",length(new.scales)),names(Lcut))]
+ names( Lcut )[1:length(new.scales)] <- new.scales
+ for( i in 1:length(wh) )
+ Lcut[,i] <- ifelse( Lcut[,timescale] - Lcut[,wh[i]] < 0,
+ NA,
+ Lcut[,timescale] - Lcut[,wh[i]] )
+ # set attributes
+ attr( Lcut, "time.scales" ) <- c( attr( Lcut, "time.scales" ), new.scales )
+ attr( Lcut, "time.since" ) <- c( attr( Lcut, "time.since" ), new.states )
+ }
+
+# return the cut object without the auxilary variable
+rmcol <- grep( "whseq", names(Lcut) )
+Lcut[,-rmcol]
+ }
+
+
+
+
+
+
diff --git a/R/plotCIF.R b/R/plotCIF.R
new file mode 100644
index 0000000..a8e0cf4
--- /dev/null
+++ b/R/plotCIF.R
@@ -0,0 +1,54 @@
+plotCIF <-
+########
+ function( x, event=1,
+ xlab = "Time", ylab = "Cumulative incidence",
+ ylim = c(0,1), lty = NULL, col = NULL,
+ ... )
+#-----------------------------------------------------------
+# Function that plots cumulative incidences for one of
+# two competing outcome events for each category of a
+# grouping factor.
+
+# Main arguments:
+# x = a survfit object with Surv( time, status, type = 'mstate')
+# event = factor for the event; first level = censoring
+# lty = vector, its length = no. of groups (default: 1 for all)
+# col = vector, its length = no. of groups (default: all "black")
+# Other arguments, like in any plot() function
+#------------------------------------------------------------
+{
+
+ng <- length(x$n)
+ne <- dim(x$pstate)[2]-1
+if (event %in% 1:ne)
+{
+time <- x$time
+gr <- rep(1, length(time) )
+if (ng > 1)
+ for (g in 2:ng)
+ for (t in (cumsum(x$strata)[g-1]+1): cumsum(x$strata)[g] )
+ gr[t] <- g
+
+lt <- rep(1, ng)
+co <- rep(1, ng)
+if (!is.null(lty)) lt = lty
+if (!is.null(col)) co = col
+
+CI <- x$pstate[, 1]
+for (e in 1:ne)
+ if (event==e) CI <- x$pstate[, e]
+
+plot( c(0, time[gr==1], max(time[gr==1]) ),
+ c(0, CI[gr==1], max(CI[gr==1]) ),
+ type = 's', ylim = ylim,
+ xlab = xlab, ylab = ylab ,
+ col = co[1], lty = lt[1], ... )
+if (ng > 1) for (g in 2:ng )
+ lines( c(0, time[gr==g], max(time[gr==g]) ),
+ c(0, CI[gr==g], max(CI[gr==g]) ),
+ type = 's', lty=lt[g], col = co[g], ... )
+}
+else print(paste("Error: event must be an integer from 1 to", ne))
+}
+
+
diff --git a/R/stackedCIF.R b/R/stackedCIF.R
new file mode 100644
index 0000000..0cd2ff4
--- /dev/null
+++ b/R/stackedCIF.R
@@ -0,0 +1,72 @@
+###################################################################
+##
+## stacked:
+## Function that draws stacked CIF curves from a survfit object
+## which has been created with Surv( ..., type = 'mstate').
+##
+## main arguments:
+## - x: survfit object the number of columns of its subobject 'prev'
+## depending on how many competing events were included
+## - group: integer showing the level of the grouping factor
+## possibly appearing after ~ in the survfit model formula
+## - colour: a character vector indicating the colours to be
+## used for shading the areas pertinent to the separate outcome
+## - ylim: as usual
+##------------------------------------------------------------------
+
+stackedCIF <- function(x, group = 1, colour = NULL, ylim = c(0,1),
+ xlab = "Time", ylab = "Cumulative incidence", ... )
+{########
+
+ne <- ncol(x$pstate)-1 # number of competing events
+ng <- length(x$n) # number of groups
+
+if (group %in% 1:ng)
+{
+r1 <- ifelse(group==1 | ng==1, 1, cumsum(x$strata)[group-1]+1 )
+r2 <- ifelse(ng==1, length(x$time), cumsum(x$strata)[group] )
+
+# Creating the matrix containing estimated probabilities of
+# different outcome states at each time point
+pSt0 <- matrix(0, nrow = nrow(x$pstate[r1:r2, ]), ncol = ne+2)
+pSt0[ , 1] <- x$pstate[r1:r2 , 1]
+for (c in 2:ne)
+ pSt0[ , c] <- pSt0[ , c-1] + x$pstate[r1:r2 , c]
+pSt0[ , ne + 1] <- 1
+pSt0[ , ne + 2] <- x$time[r1:r2]
+
+cols <- rep("white", ne)
+if (!is.null(colour)) cols = colour
+
+plot(as.numeric(pSt0[, ne + 2]), pSt0[, 2], type = "s",
+ ylim = ylim, yaxs = "i",
+ xlab = xlab, ylab = ylab, ...)
+
+for (c in 2:ne)
+lines(as.numeric(pSt0[, ne + 2]), pSt0[, c], type = "s")
+
+col <- colour
+bord = "black"
+border <- rep(bord, ne)[1:ne]
+ pSt <- cbind(0, pSt0)
+ pSt2 <- matrix( 0, nrow = 2*(dim(pSt)[1]-1), ncol = dim(pSt)[2])
+ pSt2[1, ne+3] <- pSt[1, ne+3]
+ pSt2[ nrow(pSt2), ne+3] <- pSt[ nrow(pSt), ne+3]
+ for ( j in 1:( dim(pSt)[1]-1 ) )
+ { pSt2[ 2*j-1 , ne+3] <- pSt[j, ne+3]
+ pSt2[ 2*j, ne+3] <- pSt[j+1, ne+3]
+ pSt2[ 2*j-1 , 1:(ne+2)] <- pSt[j, 1:(ne+2)]
+ pSt2[ 2*j, 1:(ne+2)] <- pSt[j, 1:(ne+2)] }
+ for (i in 2:(ne+1)) {
+ polygon(c(pSt2[, ne+3], rev(pSt2[, ne+3])),
+ c(pSt2[, i], rev(pSt2[, i - 1])), col = cols[i - 1],
+ border = NULL, ... )
+ }
+abline(v = max(x$time[r1:r2]), lwd = 2, col = "white")
+box()
+} ## end of if
+else print(paste("Error: group indicator must be an integer from 1 to", ng))
+} ## end of function
+
+#------------------------------------------------------------------
+
diff --git a/R/xgrep.R b/R/xgrep.R
new file mode 100644
index 0000000..90047a7
--- /dev/null
+++ b/R/xgrep.R
@@ -0,0 +1,3 @@
+fgrep <- function( pattern, x, ... ) x [grep( pattern, x , ... )]
+ngrep <- function( pattern, x, ... ) names(x)[grep( pattern, names(x), ... )]
+lgrep <- function( pattern, x, ... ) levels(x)[grep( pattern, levels(x), ... )]
diff --git a/R/zzz.R b/R/zzz.R
deleted file mode 100644
index 72df4cd..0000000
--- a/R/zzz.R
+++ /dev/null
@@ -1,4 +0,0 @@
-.onLoad <- function(lib, pkg)
-{
- library.dynam("Epi", pkg, lib, local=FALSE)
-}
diff --git a/build/vignette.rds b/build/vignette.rds
index 0d21d86..998f1a2 100644
Binary files a/build/vignette.rds and b/build/vignette.rds differ
diff --git a/data/DMepi.rda b/data/DMepi.rda
new file mode 100644
index 0000000..aa6fe91
Binary files /dev/null and b/data/DMepi.rda differ
diff --git a/inst/doc/Follow-up.pdf b/inst/doc/Follow-up.pdf
index 720fc81..a9c459e 100644
Binary files a/inst/doc/Follow-up.pdf and b/inst/doc/Follow-up.pdf differ
diff --git a/inst/doc/index.html b/inst/doc/index.html
index 5d67758..9da1be2 100644
--- a/inst/doc/index.html
+++ b/inst/doc/index.html
@@ -10,6 +10,9 @@
<li> <a href="simLexis.pdf">
Simulation in multistate models with multiple timescales
</a>, and the corresponding <a href="simLexis.R"> R-code </a>.
+ <li> <a href="yll.pdf">
+ Calculation of years of life lost,
+ </a>, and the corresponding <a href="yll.R"> R-code </a>.
</ul>
Here is the website for the <a href="http://BendixCarstensen.com/Epi/">Epi package</a>.
</html>
diff --git a/inst/doc/simLexis.R b/inst/doc/simLexis.R
index fe8db38..10eaad6 100644
--- a/inst/doc/simLexis.R
+++ b/inst/doc/simLexis.R
@@ -51,11 +51,11 @@ print( subset( Si, lex.id==97 )[,1:10], digits=6 )
### code chunk number 6: knots
###################################################
nk <- 5
-( ai.kn <- with( subset(Si,lex.Xst=="Ins"),
+( ai.kn <- with( subset(Si,lex.Xst=="Ins" & lex.Cst!=lex.Xst ),
quantile( Age+lex.dur , probs=(1:nk-0.5)/nk ) ) )
( ad.kn <- with( subset(Si,lex.Xst=="Dead"),
quantile( Age+lex.dur , probs=(1:nk-0.5)/nk ) ) )
-( di.kn <- with( subset(Si,lex.Xst=="Ins"),
+( di.kn <- with( subset(Si,lex.Xst=="Ins" & lex.Cst!=lex.Xst ),
c(0,quantile( DMdur+lex.dur, probs=(1:(nk-1))/nk ) )) )
( dd.kn <- with( subset(Si,lex.Xst=="Dead"),
c(0,quantile( DMdur+lex.dur, probs=(1:(nk-1))/nk ) )) )
@@ -121,7 +121,7 @@ str( pr.rates )
###################################################
-### code chunk number 11: simLexis.rnw:509-510
+### code chunk number 11: simLexis.rnw:515-516
###################################################
ci.pred
@@ -389,13 +389,13 @@ box( lwd=3 )
###################################################
-### code chunk number 27: simLexis.rnw:1027-1028
+### code chunk number 27: simLexis.rnw:1033-1034
###################################################
-options( keep.source=TRUE )
+source( "../R/simLexis.R", keep.source=TRUE )
###################################################
-### code chunk number 28: simLexis.rnw:1044-1047
+### code chunk number 28: simLexis.rnw:1050-1053
###################################################
cbind(
attr( ini, "time.scale" ),
@@ -403,55 +403,55 @@ attr( ini, "time.since" ) )
###################################################
-### code chunk number 29: simLexis.rnw:1072-1073
+### code chunk number 29: simLexis.rnw:1078-1079
###################################################
simLexis
###################################################
-### code chunk number 30: simLexis.rnw:1090-1091
+### code chunk number 30: simLexis.rnw:1096-1097
###################################################
-Epi:::simX
+simX
###################################################
-### code chunk number 31: simLexis.rnw:1103-1104
+### code chunk number 31: simLexis.rnw:1109-1110
###################################################
-Epi:::sim1
+sim1
###################################################
-### code chunk number 32: simLexis.rnw:1116-1117
+### code chunk number 32: simLexis.rnw:1122-1123
###################################################
-Epi:::lint
+lint
###################################################
-### code chunk number 33: simLexis.rnw:1127-1128
+### code chunk number 33: simLexis.rnw:1133-1134
###################################################
-Epi:::get.next
+get.next
###################################################
-### code chunk number 34: simLexis.rnw:1137-1138
+### code chunk number 34: simLexis.rnw:1143-1144
###################################################
-Epi:::chop.lex
+chop.lex
###################################################
-### code chunk number 35: simLexis.rnw:1155-1156
+### code chunk number 35: simLexis.rnw:1161-1162
###################################################
nState
###################################################
-### code chunk number 36: simLexis.rnw:1165-1166
+### code chunk number 36: simLexis.rnw:1171-1172
###################################################
pState
###################################################
-### code chunk number 37: simLexis.rnw:1170-1172
+### code chunk number 37: simLexis.rnw:1176-1178
###################################################
plot.pState
lines.pState
diff --git a/inst/doc/simLexis.pdf b/inst/doc/simLexis.pdf
index 5028e54..47a6a91 100644
Binary files a/inst/doc/simLexis.pdf and b/inst/doc/simLexis.pdf differ
diff --git a/inst/doc/simLexis.rnw b/inst/doc/simLexis.rnw
index aa12586..e5c9abf 100644
--- a/inst/doc/simLexis.rnw
+++ b/inst/doc/simLexis.rnw
@@ -366,21 +366,27 @@ knots so that the number of events is the same between each pair of
knots, but only half of this beyond each of the boundary knots,
whereas for the timescales \texttt{DMdur} and \texttt{tIns} where we
have observation from a well-defined 0, we put knots at 0 and place the
-remaining knots so that the number of events is the same between them
-and beyond the last:
+remaining knots so that the number of events is the same between each
+pair of knots as well as outside the boundary knots.
<<knots>>=
nk <- 5
-( ai.kn <- with( subset(Si,lex.Xst=="Ins"),
+( ai.kn <- with( subset(Si,lex.Xst=="Ins" & lex.Cst!=lex.Xst ),
quantile( Age+lex.dur , probs=(1:nk-0.5)/nk ) ) )
( ad.kn <- with( subset(Si,lex.Xst=="Dead"),
quantile( Age+lex.dur , probs=(1:nk-0.5)/nk ) ) )
-( di.kn <- with( subset(Si,lex.Xst=="Ins"),
+( di.kn <- with( subset(Si,lex.Xst=="Ins" & lex.Cst!=lex.Xst ),
c(0,quantile( DMdur+lex.dur, probs=(1:(nk-1))/nk ) )) )
( dd.kn <- with( subset(Si,lex.Xst=="Dead"),
c(0,quantile( DMdur+lex.dur, probs=(1:(nk-1))/nk ) )) )
( ti.kn <- with( subset(Si,lex.Xst=="Dead(Ins)"),
c(0,quantile( t.Ins+lex.dur, probs=(1:(nk-1))/nk ) )) )
@ %
+Note that when we tease out the event records for transition to
+\emph{transient} states (in this case ``Ins'', that is
+verb|lex.Xst=="Ins"|), we should add \verb|lex.Cst!=lex.Xst|, to
+include only transition records and avoiding including records of
+sojourn time in the transient state.
+
We then fit Poisson models to transition rates, using the wrapper
\texttt{Ns} from the \texttt{Epi} package to simplify the
specification of the rates:
@@ -1025,7 +1031,7 @@ follows:
This section explains the actually existing code for
\texttt{simLexis}, as it is in the current version of \texttt{Epi}.
<<echo=FALSE,results=hide>>=
-options( keep.source=TRUE )
+source( "../R/simLexis.R", keep.source=TRUE )
@ %
The function \texttt{simLexis} takes a \texttt{Lexis} object as
input. This defines the initial state(s) and times of the start for a
@@ -1088,7 +1094,7 @@ and objects of class \texttt{function} is made.
The dataset on which this prediction is done has
\texttt{length(time.pts)} rows per person.
<<>>=
-Epi:::simX
+simX
@ %
As we see, \texttt{simX} calls \texttt{sim1} which simulates the
transition for one person.
@@ -1101,7 +1107,7 @@ The predicted cumulative intensities are fed, person by person, to
appended to the \texttt{nd} data frame. This way we have simulated
\emph{one} transition (time and state) for each person:
<<>>=
-Epi:::sim1
+sim1
@ %
The \texttt{sim1} function uses \texttt{lint} to do linear interpolation.
@@ -1114,7 +1120,7 @@ a custom made linear interpolator that does the interpolation
exploiting the fact the the \texttt{ci} is non-decreasing and
\texttt{tt} is strictly monotonously increasing:
<<>>=
-Epi:::lint
+lint
@
\subsection{\texttt{get.next}}
@@ -1125,7 +1131,7 @@ defined as time since entry to one of these states be initialized to 0
before a call to \texttt{simX} is made for these persons. This
accomplished by the function \texttt{get.next}:
<<>>=
-Epi:::get.next
+get.next
@
\subsection{\texttt{chop.lex}}
@@ -1135,7 +1141,7 @@ after \emph{each} new entry to a transient state. In order to groom
the output data we use \texttt{chop.lex} to censor all persons at the
same designated time after \emph{initial} entry.
<<>>=
-Epi:::chop.lex
+chop.lex
@
\section{Probabilities from simulated \texttt{Lexis} objects}
diff --git a/inst/doc/yll.R b/inst/doc/yll.R
new file mode 100644
index 0000000..8d5310d
--- /dev/null
+++ b/inst/doc/yll.R
@@ -0,0 +1,206 @@
+### R code from vignette source 'yll.rnw'
+### Encoding: UTF-8
+
+###################################################
+### code chunk number 1: yll.rnw:146-149
+###################################################
+options( width=90,
+ SweaveHooks=list( fig=function()
+ par(mar=c(3,3,1,1),mgp=c(3,1,0)/1.6,las=1,bty="n") ) )
+
+
+###################################################
+### code chunk number 2: states
+###################################################
+getOption("SweaveHooks")[["fig"]]()
+library( Epi )
+TM <- matrix(NA,4,4)
+rownames(TM) <-
+colnames(TM) <- c("Well","DM","Dead","Dead(DM)")
+TM[1,2:3] <- TM[2,4] <- 1
+zz <- boxes( TM, boxpos=list(x=c(20,80,20,80),y=c(80,80,20,20)), wm=1.5, hm=4 )
+
+
+###################################################
+### code chunk number 3: states
+###################################################
+getOption("SweaveHooks")[["fig"]]()
+zz$Arrowtext <- c( expression(lambda),
+ expression(mu[W]),
+ expression(mu[D][M]) )
+boxes( zz )
+
+
+###################################################
+### code chunk number 4: yll.rnw:392-393
+###################################################
+data( DMepi )
+
+
+###################################################
+### code chunk number 5: yll.rnw:398-400
+###################################################
+str( DMepi )
+head( DMepi )
+
+
+###################################################
+### code chunk number 6: yll.rnw:420-424
+###################################################
+DMepi <- transform( subset( DMepi, A>30 ),
+ D.T = D.nD + D.DM,
+ Y.T = Y.nD + Y.DM )
+head(DMepi)
+
+
+###################################################
+### code chunk number 7: yll.rnw:430-456
+###################################################
+# Knots used in all models
+( a.kn <- seq(40,95,,6) )
+( p.kn <- seq(1997,2015,,4) )
+( c.kn <- seq(1910,1976,,6) )
+# Check the number of events between knots
+ae <- xtabs( cbind(D.nD,D.DM,X) ~ cut(A,c(30,a.kn,Inf)) + sex, data=DMepi )
+ftable( addmargins(ae,1), col.vars=3:2 )
+pe <- xtabs( cbind(D.nD,D.DM,X) ~ cut(P,c(1990,p.kn,Inf)) + sex, data=DMepi )
+ftable( addmargins(pe,1), col.vars=3:2 )
+ce <- xtabs( cbind(D.nD,D.DM,X) ~ cut(P-A,c(-Inf,c.kn,Inf)) + sex, data=DMepi )
+ftable( addmargins(ce,1), col.vars=3:2 )
+# Fit an APC-model for all transitions, seperately for men and women
+mW.m <- glm( D.nD ~ -1 + Ns(A ,knots=a.kn,int=TRUE) +
+ Ns( P,knots=p.kn,ref=2005) +
+ Ns(P-A,knots=c.kn,ref=1950),
+ offset = log(Y.nD),
+ family = poisson,
+ data = subset( DMepi, sex=="M" ) )
+mD.m <- update( mW.m, D.DM ~ . , offset=log(Y.DM) )
+mT.m <- update( mW.m, D.T ~ . , offset=log(Y.T ) )
+lW.m <- update( mW.m, X ~ . )
+# Model for women
+mW.f <- update( mW.m, data = subset( DMepi, sex=="F" ) )
+mD.f <- update( mD.m, data = subset( DMepi, sex=="F" ) )
+mT.f <- update( mT.m, data = subset( DMepi, sex=="F" ) )
+lW.f <- update( lW.m, data = subset( DMepi, sex=="F" ) )
+
+
+###################################################
+### code chunk number 8: yll.rnw:463-500
+###################################################
+a.ref <- 30:90
+p.ref <- 1996:2016
+aYLL <- NArray( list( type = c("Imm","Tot","Sus"),
+ sex = levels( DMepi$sex ),
+ age = a.ref,
+ date = p.ref ) )
+str( aYLL )
+system.time(
+for( ip in p.ref )
+ {
+ nd <- data.frame( A = seq(30,90,0.2)+0.1,
+ P = ip,
+ Y.nD = 1,
+ Y.DM = 1,
+ Y.T = 1 )
+ muW.m <- ci.pred( mW.m, nd )[,1]
+ muD.m <- ci.pred( mD.m, nd )[,1]
+ muT.m <- ci.pred( mT.m, nd )[,1]
+ lam.m <- ci.pred( lW.m, nd )[,1]
+ muW.f <- ci.pred( mW.f, nd )[,1]
+ muD.f <- ci.pred( mD.f, nd )[,1]
+ muT.f <- ci.pred( mT.f, nd )[,1]
+ lam.f <- ci.pred( lW.f, nd )[,1]
+ aYLL["Imm","M",,paste(ip)] <- yll( int=0.2, muW.m, muD.m, lam=NULL,
+ A=a.ref, age.in=30, note=FALSE )[-1]
+ aYLL["Imm","F",,paste(ip)] <- yll( int=0.2, muW.f, muD.f, lam=NULL,
+ A=a.ref, age.in=30, note=FALSE )[-1]
+ aYLL["Tot","M",,paste(ip)] <- yll( int=0.2, muT.m, muD.m, lam=NULL,
+ A=a.ref, age.in=30, note=FALSE )[-1]
+ aYLL["Tot","F",,paste(ip)] <- yll( int=0.2, muT.f, muD.f, lam=NULL,
+ A=a.ref, age.in=30, note=FALSE )[-1]
+ aYLL["Sus","M",,paste(ip)] <- yll( int=0.2, muW.m, muD.m, lam=lam.m,
+ A=a.ref, age.in=30, note=FALSE )[-1]
+ aYLL["Sus","F",,paste(ip)] <- yll( int=0.2, muW.f, muD.f, lam=lam.f,
+ A=a.ref, age.in=30, note=FALSE )[-1]
+ } )
+round( ftable( aYLL[,,seq(1,61,10),], col.vars=c(3,2) ), 1 )
+
+
+###################################################
+### code chunk number 9: imm
+###################################################
+getOption("SweaveHooks")[["fig"]]()
+plyll <- function(wh){
+par( mfrow=c(1,2), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, bty="n", las=1 )
+
+matplot( a.ref, aYLL[wh,"M",,],
+ type="l", lty=1, col="blue", lwd=1:2,
+ ylim=c(0,12), xlab="Age",
+ ylab="Years lost to DM", yaxs="i" )
+abline(v=50,h=1:10,col=gray(0.7))
+text( 90, 11, "Men", col="blue", adj=1 )
+text( 40, aYLL[wh,"M","40","1996"], "1996", adj=c(0,0), col="blue" )
+text( 43, aYLL[wh,"M","44","2016"], "2016", adj=c(1,1), col="blue" )
+
+matplot( a.ref, aYLL[wh,"F",,],
+ type="l", lty=1, col="red", lwd=1:2,
+ ylim=c(0,12), xlab="Age",
+ ylab="Years lost to DM", yaxs="i" )
+abline(v=50,h=1:10,col=gray(0.7))
+text( 90, 11, "Women", col="red", adj=1 )
+text( 40, aYLL[wh,"F","40","1996"], "1996", adj=c(0,0), col="red" )
+text( 43, aYLL[wh,"F","44","2016"], "2016", adj=c(1,1), col="red" )
+}
+plyll("Imm")
+
+
+###################################################
+### code chunk number 10: tot
+###################################################
+getOption("SweaveHooks")[["fig"]]()
+plyll("Tot")
+
+
+###################################################
+### code chunk number 11: sus
+###################################################
+getOption("SweaveHooks")[["fig"]]()
+plyll("Sus")
+
+
+###################################################
+### code chunk number 12: yll.rnw:584-585
+###################################################
+source( "../R/erl.R", keep.source=TRUE )
+
+
+###################################################
+### code chunk number 13: yll.rnw:595-596
+###################################################
+surv1
+
+
+###################################################
+### code chunk number 14: yll.rnw:600-601
+###################################################
+erl1
+
+
+###################################################
+### code chunk number 15: yll.rnw:608-609
+###################################################
+surv2
+
+
+###################################################
+### code chunk number 16: yll.rnw:613-614
+###################################################
+erl
+
+
+###################################################
+### code chunk number 17: yll.rnw:618-619
+###################################################
+yll
+
+
diff --git a/inst/doc/yll.pdf b/inst/doc/yll.pdf
new file mode 100644
index 0000000..3ea3355
Binary files /dev/null and b/inst/doc/yll.pdf differ
diff --git a/inst/doc/yll.rnw b/inst/doc/yll.rnw
new file mode 100644
index 0000000..5da0094
--- /dev/null
+++ b/inst/doc/yll.rnw
@@ -0,0 +1,641 @@
+\SweaveOpts{results=verbatim,keep.source=TRUE,include=FALSE,eps=FALSE}
+%\VignetteIndexEntry{Years of life lost: simLexis}
+\documentclass[a4paper,twoside,12pt]{report}
+
+% ----------------------------------------------------------------------
+% General information for the title page and the page headings
+\newcommand{\Title}{Years of Life Lost (YLL)
+ to disease\\Diabetes in DK as example}
+\newcommand{\Tit}{YLL}
+\newcommand{\Version}{Version 1.2}
+\newcommand{\Dates}{February 2017}
+\newcommand{\Where}{SDC}
+\newcommand{\Homepage}{\url{http://bendixcarstensen.com/Epi}}
+\newcommand{\Faculty}{\begin{tabular}{rl}
+Bendix Carstensen
+ & Steno Diabetes Center, Gentofte, Denmark\\
+ & {\small \& Department of Biostatistics,
+ University of Copenhagen} \\
+ & \texttt{b at bxc.dk}\\
+ & \url{http://BendixCarstensen.com} \\[1em]
+ \end{tabular}}
+% Packages
+\usepackage[utf8]{inputenc}
+\usepackage[T1]{fontenc}
+\usepackage[english]{babel}
+\usepackage[font=it,labelfont=normalfont]{caption}
+\usepackage[colorlinks,urlcolor=blue,linkcolor=red,,citecolor=Maroon]{hyperref}
+\usepackage[ae,hyper]{Rd}
+\usepackage[dvipsnames]{xcolor}
+\usepackage[super]{nth}
+\usepackage[noae]{Sweave}
+\usepackage{makeidx,floatflt,amsmath,amsfonts,amsbsy,enumitem,dcolumn,needspace}
+\usepackage{ifthen,calc,eso-pic,everyshi}
+\usepackage{booktabs,longtable,rotating,graphicx,subfig}
+\usepackage{pdfpages,verbatim,fancyhdr,datetime,%
+afterpage}
+\usepackage[abspath]{currfile}
+% \usepackage{times}
+\renewcommand{\textfraction}{0.0}
+\renewcommand{\topfraction}{1.0}
+\renewcommand{\bottomfraction}{1.0}
+\renewcommand{\floatpagefraction}{0.9}
+% \usepackage{mslapa}
+\definecolor{blaa}{RGB}{99,99,255}
+\DeclareGraphicsExtensions{.png,.pdf,.jpg}
+% Make the Sweave output nicer
+\DefineVerbatimEnvironment{Sinput}{Verbatim}{fontsize=\small,fontshape=sl,formatcom=\color{Blue}}
+\DefineVerbatimEnvironment{Soutput}{Verbatim}{fontsize=\small,formatcom=\color{Maroon},xleftmargin=0em}
+\DefineVerbatimEnvironment{Scode}{Verbatim}{fontsize=\small}
+\fvset{listparameters={\setlength{\topsep}{-0.1ex}}}
+\renewenvironment{Schunk}%
+{\renewcommand{\baselinestretch}{0.85} \vspace{\topsep}}%
+{\renewcommand{\baselinestretch}{1.00} \vspace{\topsep}}
+\providecommand{\ptxt}[1]{\Pp\left\{\text{#1}\right\}}
+\providecommand{\dif}{{\,\mathrm d}}
+\DeclareMathOperator{\YLL}{YLL}
+\DeclareMathOperator{\Pp}{P}
+
+%----------------------------------------------------------------------
+% Set up layout of pages
+\oddsidemargin 1mm
+\evensidemargin 1mm
+\topmargin -10mm
+\headheight 8mm
+\headsep 5mm
+\textheight 240mm
+\textwidth 165mm
+%\footheight 5mm
+\footskip 15mm
+\renewcommand{\topfraction}{0.9}
+\renewcommand{\bottomfraction}{0.9}
+\renewcommand{\textfraction}{0.1}
+\renewcommand{\floatpagefraction}{0.9}
+\renewcommand{\headrulewidth}{0.1pt}
+\setcounter{secnumdepth}{4}
+% \setcounter{tocdepth}{2}
+
+%----------------------------------------------------------------------
+% How to insert a figure in a .rnw file
+\newcommand{\rwpre}{./graph/gr}
+\newcommand{\insfig}[3]{
+\begin{figure}[h]
+ \centering
+ \includegraphics[width=#2\textwidth]{\rwpre-#1}
+ \caption{#3}
+ \label{fig:#1}
+% \afterpage{\clearpage}
+\end{figure}}
+
+%----------------------------------------------------------------------
+% Here is the document starting with the titlepage
+\begin{document}
+
+%----------------------------------------------------------------------
+% The title page
+\setcounter{page}{1}
+\pagenumbering{roman}
+\pagestyle{plain}
+\thispagestyle{empty}
+% \vspace*{0.05\textheight}
+\flushright
+% The blank below here is necessary in order not to muck up the
+% linespacing in title if it has more than 2 lines
+{\Huge \bfseries \Title
+
+}\ \\[-1.5ex]
+\noindent\textcolor{blaa}{\rule[-1ex]{\textwidth}{5pt}}\\[2.5ex]
+\large
+\Where \\
+\Dates \\
+\Homepage \\
+\Version \\[1em]
+\normalsize
+Compiled \today,\ \currenttime\\
+from: \texttt{\currfileabspath}\\[1em]
+% \input{wordcount}
+\normalsize
+\vfill
+\Faculty
+% End of titlepage
+% \newpage
+
+%----------------------------------------------------------------------
+% Table of contents
+\tableofcontents
+
+%----------------------------------------------------------------------
+% General text layout
+\raggedright
+\parindent 1em
+\parskip 0ex
+\cleardoublepage
+
+%----------------------------------------------------------------------
+% General page style
+\pagenumbering{arabic}
+\setcounter{page}{1}
+\pagestyle{fancy}
+\renewcommand{\chaptermark}[1]{\markboth{\textsl{#1}}{}}
+\renewcommand{\sectionmark}[1]{\markright{\thesection\ \textsl{#1}}{}}
+\fancyhead[EL]{\bf \thepage \quad \rm \leftmark}
+\fancyhead[ER,OL]{\sl \Tit}
+\fancyhead[OR]{\rm \rightmark \quad \bf \thepage}
+\fancyfoot{}
+
+<<echo=FALSE>>=
+options( width=90,
+ SweaveHooks=list( fig=function()
+ par(mar=c(3,3,1,1),mgp=c(3,1,0)/1.6,las=1,bty="n") ) )
+@ %
+\renewcommand{\rwpre}{./yll}
+
+%----------------------------------------------------------------------
+% Here comes the substance part
+\chapter{Theory and technicalities}
+
+This vignette for the \texttt{Epi} package describes the
+probabilistic/demographic background for and technical implementation
+of the \texttt{erl} and \texttt{yll} functions that computes the
+expected residual life time and years of life lost in an illness-death
+model.
+
+\section{Years of life lost (YLL)}
+
+\ldots to diabetes or any other disease for that matter.
+
+The general concept in calculation of ``years lost to\ldots'' is the
+comparison of the expected lifetime between two groups of persons; one
+with and one without disease (in this example DM). The expected
+lifetime is the area under the survival curve, so basically the
+exercise requires that two survival curves that are deemed relevant be
+available.
+
+The years of life lost is therefore just the area between the survival
+curves for those ``Well'', $S_W(t)$, and for those ``Diseased'',
+$S_D(t)$:
+\[
+ \YLL = \int_0^\infty S_W(t) - S_D(t) \dif t
+\]
+The time $t$ could of course be age, but it could also be ``time after
+age 50'' and the survival curves compared would then be survival
+curves \emph{conditional} on survival till age 50, and the YLL would
+be the years of life lost for a 50-year old person with diabetes.
+
+If we are referring to the expected lifetime we will more precisely use
+the label expected residual lifetime, ERL.
+
+\section{Constructing the survival curves}
+
+YLL can be computed in two different ways, depending on the way the
+survival curve and hence the expected lifetime of a person
+\emph{without} diabetes is computed:
+\begin{itemize}
+\item Assume that the ``Well'' persons are \emph{immune} to disease
+ --- using only the non-DM mortality rates throughout for calculation
+ of expected life time.
+\item Assume that the ``Well'' persons \emph{can} acquire the disease and
+ thereby see an increased mortality, thus involving all three rates
+ shown in figure \ref{fig:states}.
+\end{itemize}
+The former gives a higher YLL because the comparison is to persons
+assumed immune to DM (and yet with the same mortality as non-immune
+prior to diagnosis), the latter gives a more realistic picture of the
+comparison of group of persons with and without diabetes at a given
+age that can be interpreted in the real world.
+
+The differences can be illustrated by figure \ref{fig:states}; the
+immune approach corresponds to an assumption of $\lambda(t)=0$ in the
+calculation of the survival curve for a person in the ``Well'' state.
+
+Calculation of the survival of a diseased person already in the ``DM''
+state is unaffected by assumptions about $\lambda$.
+
+<<states,fig=TRUE,height=5,width=5,echo=FALSE>>=
+library( Epi )
+TM <- matrix(NA,4,4)
+rownames(TM) <-
+colnames(TM) <- c("Well","DM","Dead","Dead(DM)")
+TM[1,2:3] <- TM[2,4] <- 1
+zz <- boxes( TM, boxpos=list(x=c(20,80,20,80),y=c(80,80,20,20)), wm=1.5, hm=4 )
+@
+<<states,fig=TRUE,height=4,width=5,echo=FALSE>>=
+zz$Arrowtext <- c( expression(lambda),
+ expression(mu[W]),
+ expression(mu[D][M]) )
+boxes( zz )
+@ %
+\insfig{states}{0.7}{Illness-death model describing diabetes incidence
+ and -mortality.}
+
+\subsection{Total mortality --- a shortcut?}
+
+A practical crude shortcut could be to compare the ERL in the diabetic
+population to the ERL for the \emph{entire} population (that is use
+the total mortality ignoring diabetes status).
+
+Note however that this approach also counts the mortality of persons
+that acquired the disease earlier, thus making the comparison
+population on average more ill than the population we aim at, namely
+those well at a given time, which only then become more gradually ill.
+
+How large these effects are can however be empirically explored, as we
+shall do later.
+
+\subsection{Disease duration}
+
+In the exposition above there is no explicit provision for the effect of
+disease duration, but if we were able to devise mortality rates for
+any combination of age and duration, this could be taken into account.
+
+There are however severe limitations in this as we in principle would
+want to have duration effects as long as the age-effects --- in
+principle for all $(a,d)$ where $d\leq A$, where $A$ is the age at
+which we condition. So even if we were only to compute ERL from
+age, say, 40 we would still need duration effects up to 60 years
+(namely to age 100).
+
+The incorporation of duration effects is in principle trivial from a
+computational point of view, but we would be forced to entertain
+models predicting duration effects way beyond what is actually
+observed disease duration in any practical case.
+
+\subsection{Computing integrals}
+
+The practical calculations of survival curves, ERL and YLL involves
+calculation of (cumulative) integrals of rates and functions of these
+as we shall see below. This is easy if we have a closed form
+expression of the function, so its value may be computed at any time
+point --- this will be the case if we model rates by smooth parametric
+functions.
+
+Computing the (cumulative) integral of a function is done as follows:
+\begin{itemize}
+\item Compute the value of the function (mortality rate for example)
+ at the midpoints of a sequence of narrow equidistant intervals ---
+ for example one- or three month intervals of age, say.
+\item Take the cumulative sum of these values multiplied by the
+ interval length --- this will be a very close approximation to the
+ cumulative integral evaluated at the end of each interval.
+\item If the intervals are really small (like 1/100 year), the
+ distinction between the value at the middle and at the end of each
+ interval becomes irrelevant.
+\end{itemize}
+Note that in the above it is assumed that the rates are given in units
+corresponding to the interval length --- or more precisely, as the
+cumulative rates over the interval.
+
+\section{Survival functions in the illness-death model}
+
+The survival functions for persons in the ``Well'' state can be
+computed under two fundamentally different scenarios, depending on
+whether persons in the ``Well'' state are assumed to be immune to the
+disease ($\lambda(a)=0$) or not.
+
+\subsection{Immune approach}
+
+In this case both survival functions for person in the two states are
+the usual simple transformation of the cumulative mortality rates:
+\[
+ S_W(a) = \exp\left(-\int_0^a\!\!\mu_W(u) \dif u \right), \qquad
+ S_D(a) = \exp\left(-\int_0^a\!\!\mu_D(u) \dif u \right)
+\]
+
+\subsubsection{Conditional survival functions}
+
+If we want the \emph{conditional} survival functions given survival to
+age $A$, say, they are just:
+\[
+ S_W(a|A) = S_W(a)/S_W(A), \qquad S_D(a|A) = S_D(a)/S_D(A)
+\]
+
+\subsection{Non-immune approach}
+
+For a diseased person, the survival function in this states is the same
+as above, but the survival function for a person without disease (at
+age 0) is (see figure \ref{fig:states}):
+\[
+S(a) = \ptxt{Well}\!(a) + \ptxt{DM}\!(a)
+\]
+In the appendix of the paper \cite{Carstensen.2008c} is an indication
+of how to compute the probability of being in any of the four states
+shown in figure \ref{fig:states}, which I shall repeat here:
+
+In terms of the rates, the probability of being in the ``Well'' box is
+simply the probability of escaping both death (at a rate of $\mu_W(a)$)
+and diabetes (at a rate of $\lambda(a)$):
+\[
+ \ptxt{Well}(a) = \exp\left(\!-\int_0^a\!\!\mu_W(u)+\lambda(u) \right) \dif u
+\]
+The probability of being alive with diabetes at age $a$, is computed given that
+ diabetes occurred at age $s$ ($s<a$) and then integrated over $s$ from $0$
+ to $a$:
+\begin{align*}
+ \ptxt{DM}(a) = \int_0^a\!\! & \ptxt{survive to $s$, DM diagnosed at $s$} \\
+ & \times \ptxt{survive with DM from $s$ to $a$} \dif s \\
+ = \int_0^a\!\! & \lambda(s)
+ \exp\left(\!-\int_0^s\!\!\mu_W(u)+\lambda(u) \dif u \right) \\
+ & \times \exp\left(\!-\int_s^a\!\!\mu_D(u) \dif u \right) \dif s
+\end{align*}
+Sometimes we will use a version where the mortality among diabetes
+patients depend both on age $a$ and duration of diabetes, $d$,
+$\mu_D(a,d)$, in which case we get:
+\begin{align*}
+ \ptxt{DM}(a) = \int_0^a \! & \lambda(s)
+ \exp\left(-\int_0^s\!\mu_W(u)+\lambda(u) \dif u \right) \\
+ & \times \exp\left(-\int_s^a\!\mu_D(u,u-s) \dif u \right) \dif s
+\end{align*}
+because the integration variable $u$ is the age-scale and the second
+integral refers to mortality among persons diagnosed at age $s$, that
+is, with duration $u-s$ at age $u$.
+
+The option of using duration-dependent mortality rates among diseased
+individuals is not implemented yet.
+
+\subsubsection{Conditional survival functions}
+
+Unlike the immune approach, the conditional survival function in the
+more realistic case is not just a ratio of the unconditional to the
+value at the conditioning age, $A$, say. This would amount to
+conditioning on being merely \emph{alive} at age $A$, but what we want
+is to condition on being in the ``Well'' state at age $A$.
+
+The formulae for the conditional probabilities of being either in
+``Well'' or ``DM'', given being in ``Well'' at age $A$ are basically
+replicates of the unconditional, albeit with changes in integration
+limits:
+\begin{align*}
+\ptxt{Well|Well at $A$}(a) &= \exp\left(-\int_A^a \! \mu_W(u)+\lambda(u) \right) \dif u \\
+ \ptxt{DM|Well at $A$}(a) &= \int_A^a \! \lambda(s)
+ \exp\left(-\int_A^s\!\mu_W(u)+\lambda(u) \dif u \right) \\
+ & \qquad \times \exp\left(-\int_s^a\!\mu_D(u,u-s) \dif u \right) \dif s
+\end{align*}
+The calculation of these conditional survival functions is implemented
+but not allowing for duration-dependence. Thus it is only implemented
+assuming $\mu_D(a,d)=\mu_D(a)$.
+
+\chapter{Analyses for DM in Denmark}
+
+The rates we use as basis for the following calculations are derived
+from the NDR, where we have omitted the blood-glucose criteria,
+because there is compelling evidence that these have quite a low
+specificity (particularly in the younger ages among women), and do
+not substantially contribute to the sensitivity.
+
+As noted above the calculations of YLL requires access to
+(age-specific) rates of incidence of DM and mortality for persons with
+and without DM.
+
+\section{Modeling mortality and incidence data}
+
+We read in the dataset of DM and population mortality and incidence, \texttt{DMepi}:
+<<>>=
+data( DMepi )
+@ %
+The dataset \texttt{DMepi} contains diabetes events, deaths and
+person-years for persons without diabetes and deaths and person-years
+for persons with diabetes:
+<<>>=
+str( DMepi )
+head( DMepi )
+@ %
+For each combination of sex, age, period and date of birth in 1 year
+age groups, we have the person-years in the ``Well'' (\texttt{Y.nD})
+and the ``DM'' (\texttt{Y.DM}) states, as well as the number of deaths
+from these (\texttt{D.nD}, \texttt{D.DM}) and the number of incident
+diabetes cases from the ``Well'' state (\texttt{X}).
+
+In order to compute the years of life lost to diabetes and how this
+has changed over time, we fit models for the mortality and incidence
+of both groups (and of course, separately for men and women). The
+models we use will be age-period-cohort models \cite{Carstensen.2007a}
+providing estimated mortality rates for ages 0--99 and dates
+1.1.1996--1.1.2016.
+
+First we transform the age and period variables to reflect the mean
+age and period in each of the Lexis triangles. We also compute the
+total number of deaths and amount of risk time, as we are going to
+model the total mortality as well. Finally we restrict the dataset to
+ages over 30 only:
+<<>>=
+DMepi <- transform( subset( DMepi, A>30 ),
+ D.T = D.nD + D.DM,
+ Y.T = Y.nD + Y.DM )
+head(DMepi)
+@ %
+With the correct age and period coding in the Lexis triangles, we fit
+models for the mortalities and incidences. Note that we for
+comparative purposes also fit a model for the \emph{total} mortality,
+ignoring the
+<<>>=
+# Knots used in all models
+( a.kn <- seq(40,95,,6) )
+( p.kn <- seq(1997,2015,,4) )
+( c.kn <- seq(1910,1976,,6) )
+# Check the number of events between knots
+ae <- xtabs( cbind(D.nD,D.DM,X) ~ cut(A,c(30,a.kn,Inf)) + sex, data=DMepi )
+ftable( addmargins(ae,1), col.vars=3:2 )
+pe <- xtabs( cbind(D.nD,D.DM,X) ~ cut(P,c(1990,p.kn,Inf)) + sex, data=DMepi )
+ftable( addmargins(pe,1), col.vars=3:2 )
+ce <- xtabs( cbind(D.nD,D.DM,X) ~ cut(P-A,c(-Inf,c.kn,Inf)) + sex, data=DMepi )
+ftable( addmargins(ce,1), col.vars=3:2 )
+# Fit an APC-model for all transitions, seperately for men and women
+mW.m <- glm( D.nD ~ -1 + Ns(A ,knots=a.kn,int=TRUE) +
+ Ns( P,knots=p.kn,ref=2005) +
+ Ns(P-A,knots=c.kn,ref=1950),
+ offset = log(Y.nD),
+ family = poisson,
+ data = subset( DMepi, sex=="M" ) )
+mD.m <- update( mW.m, D.DM ~ . , offset=log(Y.DM) )
+mT.m <- update( mW.m, D.T ~ . , offset=log(Y.T ) )
+lW.m <- update( mW.m, X ~ . )
+# Model for women
+mW.f <- update( mW.m, data = subset( DMepi, sex=="F" ) )
+mD.f <- update( mD.m, data = subset( DMepi, sex=="F" ) )
+mT.f <- update( mT.m, data = subset( DMepi, sex=="F" ) )
+lW.f <- update( lW.m, data = subset( DMepi, sex=="F" ) )
+@ %
+
+\section{Residual life time and years lost to DM}
+
+We now collect the estimated years of life lost classified by method
+(immune assumption or not), sex, age and calendar time:
+<<>>=
+a.ref <- 30:90
+p.ref <- 1996:2016
+aYLL <- NArray( list( type = c("Imm","Tot","Sus"),
+ sex = levels( DMepi$sex ),
+ age = a.ref,
+ date = p.ref ) )
+str( aYLL )
+system.time(
+for( ip in p.ref )
+ {
+ nd <- data.frame( A = seq(30,90,0.2)+0.1,
+ P = ip,
+ Y.nD = 1,
+ Y.DM = 1,
+ Y.T = 1 )
+ muW.m <- ci.pred( mW.m, nd )[,1]
+ muD.m <- ci.pred( mD.m, nd )[,1]
+ muT.m <- ci.pred( mT.m, nd )[,1]
+ lam.m <- ci.pred( lW.m, nd )[,1]
+ muW.f <- ci.pred( mW.f, nd )[,1]
+ muD.f <- ci.pred( mD.f, nd )[,1]
+ muT.f <- ci.pred( mT.f, nd )[,1]
+ lam.f <- ci.pred( lW.f, nd )[,1]
+ aYLL["Imm","M",,paste(ip)] <- yll( int=0.2, muW.m, muD.m, lam=NULL,
+ A=a.ref, age.in=30, note=FALSE )[-1]
+ aYLL["Imm","F",,paste(ip)] <- yll( int=0.2, muW.f, muD.f, lam=NULL,
+ A=a.ref, age.in=30, note=FALSE )[-1]
+ aYLL["Tot","M",,paste(ip)] <- yll( int=0.2, muT.m, muD.m, lam=NULL,
+ A=a.ref, age.in=30, note=FALSE )[-1]
+ aYLL["Tot","F",,paste(ip)] <- yll( int=0.2, muT.f, muD.f, lam=NULL,
+ A=a.ref, age.in=30, note=FALSE )[-1]
+ aYLL["Sus","M",,paste(ip)] <- yll( int=0.2, muW.m, muD.m, lam=lam.m,
+ A=a.ref, age.in=30, note=FALSE )[-1]
+ aYLL["Sus","F",,paste(ip)] <- yll( int=0.2, muW.f, muD.f, lam=lam.f,
+ A=a.ref, age.in=30, note=FALSE )[-1]
+ } )
+round( ftable( aYLL[,,seq(1,61,10),], col.vars=c(3,2) ), 1 )
+@ %
+We now have the relevant points for the graph showing YLL to diabetes
+for men and women by age, and calendar year, both under the immunity
+and susceptibility models for the calculation of YLL.
+<<imm,fig=TRUE,width=8,height=5>>=
+plyll <- function(wh){
+par( mfrow=c(1,2), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, bty="n", las=1 )
+
+matplot( a.ref, aYLL[wh,"M",,],
+ type="l", lty=1, col="blue", lwd=1:2,
+ ylim=c(0,12), xlab="Age",
+ ylab="Years lost to DM", yaxs="i" )
+abline(v=50,h=1:10,col=gray(0.7))
+text( 90, 11, "Men", col="blue", adj=1 )
+text( 40, aYLL[wh,"M","40","1996"], "1996", adj=c(0,0), col="blue" )
+text( 43, aYLL[wh,"M","44","2016"], "2016", adj=c(1,1), col="blue" )
+
+matplot( a.ref, aYLL[wh,"F",,],
+ type="l", lty=1, col="red", lwd=1:2,
+ ylim=c(0,12), xlab="Age",
+ ylab="Years lost to DM", yaxs="i" )
+abline(v=50,h=1:10,col=gray(0.7))
+text( 90, 11, "Women", col="red", adj=1 )
+text( 40, aYLL[wh,"F","40","1996"], "1996", adj=c(0,0), col="red" )
+text( 43, aYLL[wh,"F","44","2016"], "2016", adj=c(1,1), col="red" )
+}
+plyll("Imm")
+@ %
+<<tot,fig=TRUE,width=8,height=5>>=
+plyll("Tot")
+@ %
+<<sus,fig=TRUE,width=8,height=5>>=
+plyll("Sus")
+@ %
+\begin{figure}[h]
+ \centering
+ \includegraphics[width=\textwidth]{yll-imm}
+ \caption{Years of life lost to DM: the difference in expected
+ residual life time at different ages between persons with and
+ without diabetes, assuming the persons without diabetes at a given
+ age remain free from diabetes (immunity assumption --- not
+ reasonable). The lines refer to date of evaluation; the top lines
+ refer to 1.1.1996 the bottom ones to 1.1.2016. Blue curves are
+ men, red women.}
+ \label{fig:imm}
+\end{figure}
+
+\begin{figure}[h]
+ \centering
+ \includegraphics[width=\textwidth]{yll-sus}
+ \caption{Years of life lost to DM: the difference in expected
+ residual life time at different ages between persons with and
+ without diabetes, allowing the persons without diabetes at a given
+ to contract diabetes and thus be subject to higher mortality. The
+ lines refer to date of evaluation; the top lines refer to 1.1.1996
+ the bottom ones to 1.1.2016. Blue curves are men, red women.}
+ \label{fig:sus}
+\end{figure}
+
+\begin{figure}[h]
+ \centering
+ \includegraphics[width=\textwidth]{yll-tot}
+ \caption{Years of life lost to DM: the difference in expected
+ residual life time at different ages between persons with and
+ without diabetes. Allowance for susceptibility is approximated by
+ using the total population mortality instead of non-DM
+ mortality. The lines refer to date of evaluation; the top lines
+ refer to 1.1.1996 the bottom ones to 1.1.2016. Blue curves are
+ men, red women.}
+ \label{fig:tot}
+\end{figure}
+
+From figure \ref{fig:sus} we see that for men aged 50 the years lost to
+diabetes has decreased from a bit over 8 to a bit less than 6 years,
+and for women from 8.5 to 5 years; so a greater improvement for women.
+
+\chapter{Practical implementation}
+
+We have devised functions that wraps these formulae up for practical
+use.
+
+\section{Function definitions}
+
+<<echo=FALSE,results=hide,eval=TRUE>>=
+source( "../R/erl.R", keep.source=TRUE )
+@
+When using the functions it is assumed that the functions $\mu_W$,
+$\mu_D$ and $\lambda$ are given as vectors corresponding to
+equidistantly (usually tightly) spaced ages from 0 to $K$ where K is
+the age where everyone can safely be assumed dead.
+
+\texttt{surv1} is a simple function that computes the survival
+function from a vector of mortality rates, and optionally the
+conditional survival given being alive at prespecified ages:
+<<>>=
+surv1
+@ %
+\texttt{erl1} basically just expands the result of \texttt{surv1} with
+a column of expected residual life times:
+<<>>=
+erl1
+@ %
+We also define a function, \texttt{surv2}, that computes the survival
+function for a non-diseased person that may become diseased with rate
+\texttt{lam} and after that die at a rate of \texttt{muD}
+(corresponding to the formulae above). This is the sane way of
+handling years of life lost to a particular illness:
+<<>>=
+surv2
+@ %
+Finally we devised a function using these to compute the expected
+residual lifetime at select ages:
+<<>>=
+erl
+@ %
+\ldots and a wrapper for this if we only want the years of life lost
+returned:
+<<>>=
+yll
+@ %
+
+\bibliographystyle{plain}
+
+\begin{thebibliography}{1}
+
+\bibitem{Carstensen.2007a}
+B~Carstensen.
+\newblock Age-{P}eriod-{C}ohort models for the {L}exis diagram.
+\newblock {\em Statistics in Medicine}, 26(15):3018--3045, July 2007.
+
+\bibitem{Carstensen.2008c}
+B.~Carstensen, J.K. Kristensen, P.~Ottosen, and K.~Borch-Johnsen.
+\newblock The {D}anish {N}ational {D}iabetes {R}egister: {T}rends in incidence,
+ prevalence and mortality.
+\newblock {\em Diabetologia}, 51:2187--2196, 2008.
+
+\end{thebibliography}
+
+\addcontentsline{toc}{chapter}{References}
+
+\end{document}
diff --git a/man/DMepi.Rd b/man/DMepi.Rd
new file mode 100644
index 0000000..e6cd779
--- /dev/null
+++ b/man/DMepi.Rd
@@ -0,0 +1,41 @@
+\name{DMepi}
+\alias{DMepi}
+\docType{data}
+\title{Epidmiological rates for diabetes in Denmark 1996--2015}
+\description{Register based counts and person-uears for incidece of
+ diabetes and mortality with and without diabetes.
+}
+\usage{data("DMepi")}
+\format{
+ A data frame with 4000 observations on the following 8 variables.
+ \describe{
+ \item{\code{sex}}{a factor with levels \code{M} \code{F}}
+ \item{\code{A}}{Age glass 0 -- 99}
+ \item{\code{P}}{Calendar year, 1996-2015}
+ \item{\code{X}}{Number of new diagnoses of diabetes}
+ \item{\code{D.nD}}{Number of deaths among persons without diabetes}
+ \item{\code{Y.nD}}{Person-years among persons without diabetes}
+ \item{\code{D.DM}}{Number of deaths among persons with diabetes}
+ \item{\code{Y.DM}}{Person-years among persons with diabetes}
+ }
+}
+\details{Based on registers of the Danish population. Only included for
+ illustrative purposes. Cannot be used as scientifically validaed data.
+}
+\examples{
+data(DMepi)
+# Total deaths and person-years in the Danish population
+ftable( addmargins( xtabs( cbind( Deaths=D.nD+D.DM,
+ PYrs=Y.nD+Y.DM ) ~ P + sex,
+ data=DMepi ),
+ 2 ),
+ row.vars = 1 )
+# Deaths and person-years in the population of diabetes patients
+round(
+ftable( addmargins( xtabs( cbind( Deaths=D.DM,
+ PYrs=Y.DM ) ~ P + sex,
+ data=DMepi ),
+ 2 ),
+ row.vars = 1 ) )
+}
+\keyword{datasets}
diff --git a/man/Icens.Rd b/man/Icens.Rd
index 285cc30..5b5c313 100644
--- a/man/Icens.Rd
+++ b/man/Icens.Rd
@@ -62,7 +62,7 @@
}
\author{
Martyn Plummer, \email{plummer at iarc.fr},
- Bendix Carstensen, \email{bxc at steno.dk} }
+ Bendix Carstensen, \email{b at bxc.dk} }
\seealso{
\code{\link{fit.add}}
\code{\link{fit.mult}}
diff --git a/man/LCa.fit.Rd b/man/LCa.fit.Rd
index 505dbec..d86cef4 100644
--- a/man/LCa.fit.Rd
+++ b/man/LCa.fit.Rd
@@ -24,14 +24,14 @@ LCa.fit( data, A, P, D, Y,
a.ref, # age reference for the interactions
pi.ref = a.ref, # age reference for the period interaction
ci.ref = a.ref, # age reference for the cohort interaction
- p.ref, # period reference for the intercation
+ p.ref, # period reference for the interaction
c.ref, # cohort reference for the interactions
npar = c(a = 6, # no. knots for main age-effect
- p = 6, # no. knots for peroid-effect
+ p = 6, # no. knots for period-effect
c = 6, # no. knots for cohort-effect
pi = 6, # no. knots for age in the period interaction
ci = 6), # no. knots for age in the cohort interaction
- VC = TRUE, # numerical calculation of the Hessia?
+ VC = TRUE, # numerical calculation of the Hessian?
alpha = 0.05, # 1 minus confidence level
eps = 1e-6, # convergence criterion
maxit = 100, # max. no iterations
@@ -81,9 +81,9 @@ LCa.fit( data, A, P, D, Y,
Reference period for the time-interaction term \code{kp(t)} where \code{kc(c.ref)=0}.
}
\item{model}{
- Character, either \code{APa} which is the classical Lee-Carter model
- for log-rates, other possibilities are \code{ACa}, \code{APCa},
- \code{APaC} or \code{APaCa}, see details.
+ Character, either \code{"APa"} which is the classical Lee-Carter model
+ for log-rates, other possibilities are \code{"ACa"}, \code{"APCa"},
+ \code{"APaC"} or \code{"APaCa"}, see details.
}
\item{npar}{
A (possibly named) vector or list, with either the number of knots or
@@ -113,7 +113,7 @@ LCa.fit( data, A, P, D, Y,
Maximal number of iterations.
}
\item{quiet}{
- Shall I shup up or talk extensively to you about iteration progression etc.?
+ Shall I shut up or talk extensively to you about iteration progression etc.?
}
\item{object}{An \code{LCa} object, see under "Value".}
\item{show.est}{Logical. Should the estimates be printed?}
@@ -139,8 +139,8 @@ LCa.fit( data, A, P, D, Y,
multiplicative age by cohort or even both. Thus the most extensive
model is:
- \deqn{\log(lambda(a,p)) = f(a) + b_p(a)k_p(a) + b_c(a)k_c(a)}%
- {log(lambda(a,p)) = f(a) + b_p(a)k_p(a) + b_c(a)k_c(a)}
+ \deqn{\log(\lambda(a,p)) = f(a) + b_p(a)k_p(a) + b_c(a)k_c(a)}{%
+ log(lambda(a,p)) = f(a) + b_p(a)k_p(a) + b_c(a)k_c(a)}
The naming convention for the models is a capital \code{P} and/or
\code{C} if the effect is in the model followed by a lower case
@@ -156,7 +156,7 @@ LCa.fit( data, A, P, D, Y,
The coefficients and the variance-covariance matrix of these are used
in \code{predict.LCa} for a parametric bootstrap (that is, a simulation from
- a multivariate normal with mean equal to the parameterestiamtes and
+ a multivariate normal with mean equal to the parameter estimates and
variance as the estimated variance-covariance) to get confidence
intervals for the predictions if \code{sim} is \code{TRUE} --- which
it is by default if they are part of the object.
@@ -201,7 +201,7 @@ LCa.fit( data, A, P, D, Y,
\item{vcov}{Variance-covariance matrix of coefficients from both
models, in the same order as in the \code{coef}. Only present if
\code{LCa.fit} were called with \code{VC=TRUE}.}
-\item{knots}{List of vectors of knots used in for the age, perido and
+\item{knots}{List of vectors of knots used in for the age, period and
cohort effects.}
\item{refs}{List of reference points used for the age, period and
cohort terms in the interactions.}
@@ -222,15 +222,22 @@ LCa.fit( data, A, P, D, Y,
\code{newdata}. If \code{LCa.fit} were called with \code{VC=TRUE}
there will be 3 columns, namely prediction (1st column) and c.i.s
based on a simulation of parameters from a multivariate normal with
- mean \code{coef} abd cariance \code{vcov} using the median and
+ mean \code{coef} and variance \code{vcov} using the median and
\code{alpha}/2 quantiles from the \code{sim} simulations. If
\code{LCa.fit} were called with \code{VC=FALSE} there
will be 6 columns, namely estimates and c.i.s from age-time model (\code{mod.at}), and
from the age-interaction model (\code{mod.b}), both using conditional variances.
}
\author{
-Bendix Carstensen, \url{http://BendixCarstensen.com}
+ Bendix Carstensen, \url{http://BendixCarstensen.com}
+
+ This function was conceived during a course on APC models at the Max
+ Planck Institute of Demographic Research (MPIDR,
+ \url{https://www.demogr.mpg.de/en/}) in Rostock in May 2016
+ (\url{http://bendixcarstensen.com/APC/MPIDR-2016/}), and finished
+ during a research stay there, kindly sponsored by the MPIDR.
}
+
\seealso{
\code{\link{apc.fit}},
\code{\link{apc.LCa}},
@@ -238,8 +245,7 @@ Bendix Carstensen, \url{http://BendixCarstensen.com}
\code{\link[demography]{lca}} }
\examples{
library( Epi )
-# Load the male lung cancer data by Lexis triangles and
-# rename age and period as required by LCa.fit
+# Load the testis cancer data by Lexis triangles
data( testisDK )
tc <- subset( testisDK, A>14 & A<60 )
head( tc )
@@ -248,7 +254,7 @@ head( tc )
tc$Y <- tc$Y / 10^5
# Fit the Lee-Carter model with age-period interaction (default)
-LCa.tc <- LCa.fit( tc, model="ACa", a.ref=30, p.ref=1980, quiet=FALSE, eps=10e-5, maxit=100 )
+LCa.tc <- LCa.fit( tc, model="ACa", a.ref=30, p.ref=1980, quiet=FALSE, eps=10e-4, maxit=50 )
LCa.tc
summary( LCa.tc )
diff --git a/man/Lexis.Rd b/man/Lexis.Rd
index c2b3c59..b131c23 100644
--- a/man/Lexis.Rd
+++ b/man/Lexis.Rd
@@ -153,7 +153,9 @@ Lexis( entry = list( per=entry ),
\code{\link{plot.Lexis}},
\code{\link{splitLexis}},
\code{\link{cutLexis}},
+ \code{\link{mcutLexis}},
\code{\link{merge.Lexis}},
+ \code{\link{addCov.Lexis}},
\code{\link{subset.Lexis}},
\code{\link{cbind.Lexis}},
\code{\link{rbind.Lexis}},
diff --git a/man/Ns.Rd b/man/Ns.Rd
index 0de2c34..1753ec1 100644
--- a/man/Ns.Rd
+++ b/man/Ns.Rd
@@ -55,7 +55,7 @@ Ns( x, ref = NULL, df = NULL,
\code{diag(weight)} (an assumption which is checked).
}
\author{
- Bendix Carstensen \email{bxc at steno.dk},
+ Bendix Carstensen \email{b at bxc.dk},
Lars Jorge D\'iaz, Steno Diabetes Center.
}
\note{
diff --git a/man/addCov.Lexis.Rd b/man/addCov.Lexis.Rd
new file mode 100644
index 0000000..86e8c08
--- /dev/null
+++ b/man/addCov.Lexis.Rd
@@ -0,0 +1,133 @@
+\name{addCov.Lexis}
+\alias{addCov.Lexis}
+\title{
+ Add covariates (typically clinical measurements) taken at known times
+ to a Lexis object.
+}
+\description{
+ When follow-up in a multistate model is represented in a
+ \code{\link{Lexis}} object we may want to add information on
+ covariates, for example clinical measurements, obtained at different
+ times. This function cuts the follow-up time (see
+ \code{\link{cutLexis}}) at the times of measurement and carries the
+ measurements forward in time to the next measurement occasion.
+}
+\usage{
+\method{addCov}{Lexis}( Lx,
+ clin,
+ timescale = 1,
+ exnam,
+ tfc = "tfc",
+ addScales = FALSE )
+}
+\arguments{
+ \item{Lx}{
+ A Lexis object with follow-up of a cohort.
+ }
+ \item{clin}{
+ A data frame with the covariates to add (typically clinical
+ measurements). Must contain a variable \code{lex.id} identifying the
+ persons represented in \code{Lx}, as well as a variable with the
+ same name as one of the \code{\link{timeScales}} in \code{Lx},
+ identifying the time at which covariates are measured.
+
+ The times must be unique within each person; if not records with
+ duplicate times are discarded, and a warning issued. This is done
+ using \code{duplicated}, so not very well-defined, it is advisable
+ that you do this by yourself.
+ }
+ \item{timescale}{
+ Numerical or character. Number or name of a timescale in
+ \code{Lx}. The \code{clin} data frame must have a variable of this
+ name indicating the time at which the covariate measurements were
+ taken.
+ }
+ \item{exnam}{
+ Character. Name of the variable in \code{clin} with the examination
+ names (such as \code{wave1}, \code{wave2} etc.). Values may not be
+ repeated within person. Will be carried over to the resulting
+ \code{Lexis} object. If there is no variable of this name in
+ \code{clin} it will be constructed; if argument omitted, a variable
+ called \code{exnam} with values \code{ex.1}, \code{ex.2} etc. will
+ be constructed.
+ }
+ \item{tfc}{
+ Character (\code{t}ime \code{f}rom \code{c}ovariate). Name of the
+ variable in the result which will contain the time since the most
+ recent covariate date. This is not a time scale as it is reset to
+ 0 at each new covariate time. Also note that by this very token,
+ this variable will be meaningless if you \code{splitLexis}
+ \emph{after} using \code{addCov.Lexis}.
+ }
+ \item{addScales}{
+ Logical. Should timescales representing time since each covariate
+ time be added? They will be named \code{paste("tf",exnam)}. Not
+ implemented, argument currently ignored.
+ }
+}
+\value{
+ A \code{Lexis} object representing the same follow-up as \code{Lx},
+ with cuts added at the times of examination, and covariate
+ measurements added for all records representing follow-up after the
+ most recent time of measurement.
+}
+\details{
+ The implementation is clumpy, the function is slow.
+}
+\author{
+ Bendix Carstensen, \email{b at bxc.dk}, \url{http://BendixCarstensen.com}
+}
+\seealso{
+\code{\link{cutLexis}},
+\code{\link{mcutLexis}},
+\code{\link{splitLexis}},
+\code{\link{Lexis}}
+}
+\examples{
+# A small bogus cohort
+xcoh <- structure( list( id = c("A", "B", "C"),
+ birth = c("1952-07-14", "1954-04-01", "1987-06-10"),
+ entry = c("1965-08-04", "1972-09-08", "1991-12-23"),
+ exit = c("1997-06-27", "1995-05-23", "1998-07-24"),
+ fail = c(1, 0, 1) ),
+ .Names = c("id", "birth", "entry", "exit", "fail"),
+ row.names = c("1", "2", "3"),
+ class = "data.frame" )
+
+# Convert the character dates into numerical variables (fractional years)
+xcoh$bt <- cal.yr( xcoh$birth )
+xcoh$en <- cal.yr( xcoh$entry )
+xcoh$ex <- cal.yr( xcoh$exit )
+
+# Define as Lexis object with timescales calendar time and age
+Lcoh <- Lexis( entry = list( per=en ),
+ exit = list( per=ex, age=ex-bt ),
+ exit.status = factor( fail, 0:1, c("Alive","Dead") ),
+ data = xcoh )
+str( Lcoh )
+Lx <- Lcoh[,1:7]
+
+# Data frame with clinical examination data, date of examination in per
+clin <- data.frame( lex.id = c(1,1,3,2),
+ per = c(1977.3,1971.7,1996.2,1990.6),
+ bp = c(120,140,160,157),
+ chol = c(5,7,8,9) )
+Lx
+clin
+
+# Works with time split BEFORE adding clinical data:
+Lb <- addCov.Lexis( splitLexis( Lx,
+ time.scale="age",
+ breaks=seq(0,80,5) ),
+ clin,
+ exnam="clX" )
+Lb
+# With time split AFTER adding clinincal data, variable tfc is largely meaningless:
+La <- splitLexis( addCov.Lexis( Lx,
+ clin,
+ exnam="clX" ),
+ breaks=seq(0,80,5),
+ time.scale="age" )
+La
+}
+\keyword{ survival }
diff --git a/man/cal.yr.Rd b/man/cal.yr.Rd
index 9a1f5f6..ddfa704 100644
--- a/man/cal.yr.Rd
+++ b/man/cal.yr.Rd
@@ -47,7 +47,7 @@
}
\author{
Bendix Carstensen, Steno Diabetes Center \& Dept. of Biostatistics,
- University of Copenhagen, \email{bxc at steno.dk},
+ University of Copenhagen, \email{b at bxc.dk},
\url{http://BendixCarstensen.com}
}
\seealso{
diff --git a/man/ci.cum.Rd b/man/ci.cum.Rd
index f8ecb0f..b38ba61 100644
--- a/man/ci.cum.Rd
+++ b/man/ci.cum.Rd
@@ -38,17 +38,20 @@ ci.cum( obj,
\item{sample}{Should a sample of the original parameters be used to
compute a cumulative rate?}
}
-\details{
- The purpose of this function is to compute cumulative rate based on a
- model for the rates. If the model is a multiplicative model for the
- rates, the purpose of \code{ctr.mat} is to return a vector of rates or
- log-rates when applied to the coefficients of the model. If log-rates
- are returned from the model, the they should be exponentiated before
+\details{
+ The purpose of this function is to the compute cumulative rate
+ (integrated intensity) at a set of points based on a model for the
+ rates.
+ \code{ctr.mat} is a matrix which, when premultiplied to the parameters
+ of the model reurn the (log)rates at a set of increasing time points. If log-rates are
+ returned from the model, the they should be exponentiated before
cumulated, and the variances computed accordingly. Since the primary
use is for log-linear Poisson models the \code{Exp} parameter defaults
to TRUE.
- The \code{ci.Exp} argumen ensures that the confidence intervals for
+ The \code{ci.Exp} argument ensures that the confidence intervals for
+ the cumulaive rates are alwys positive, so that exp(-cum.rate) is
+ always in [0,1].
}
\value{
A matrix with 4 columns: Estimate, lower and upper c.i. and standard
@@ -91,10 +94,14 @@ MM <- ns( tB, knots=c(50,100,200,400,700), intercept=TRUE )
mp <- glm( D ~ MM - 1 + offset(log(Y)),
family=poisson, eps=10^-8, maxit=25 )
+# mp is now a model for the rates along the time scale tB
+
# Contrast matrix to extract effects, i.e. matrix to multiply with the
# coefficients to produce the log-rates: unique rows of MM, in time order.
T.pt <- sort( unique( tB ) )
T.wh <- match( T.pt, tB )
+
+# ctr.mat=MM[T.wh,] selects the rates as evaluated at times T.pt:
Lambda <- ci.cum( mp, ctr.mat=MM[T.wh,], intl=diff(c(0,T.pt)) )
# Put the estimated survival function on top of the KM-estimator
diff --git a/man/ci.lin.Rd b/man/ci.lin.Rd
index 2a4e9e6..f4cb47b 100644
--- a/man/ci.lin.Rd
+++ b/man/ci.lin.Rd
@@ -3,6 +3,7 @@
\alias{ci.mat}
\alias{ci.exp}
\alias{ci.pred}
+\alias{ci.ratio}
\alias{Wald}
\title{
Compute linear functions of parameters with standard errors and
@@ -30,8 +31,13 @@ Wald( obj, H0=0, ... )
ci.mat( alpha = 0.05, df = Inf )
ci.pred( obj, newdata,
Exp = NULL,
+ alpha = 0.05 )
+ci.ratio( r1, r2,
+ se1 = NULL,
+ se2 = NULL,
+ log.tr = !is.null(se1) & !is.null(se2),
alpha = 0.05,
- df = Inf )
+ pval = FALSE )
}
\arguments{
\item{obj}{A model object (of class
@@ -85,8 +91,8 @@ ci.pred( obj, newdata,
function; if FALSE, the prediction on the link scale is returned.}
\item{sample}{Logical or numerical. If \code{TRUE} or numerical a
sample of size \code{as.numeric(sample)} is drawn from the
- multivariate normal with mean equal to the (\code{subet} defined)
- coefficients and variance equal to the estimated variance-covriance
+ multivariate normal with mean equal to the (\code{subset} defined)
+ coefficients and variance equal to the estimated variance-covariance
of these. These are then transformed by \code{ctr.mat} and
returned.}
\item{pval}{Logical. Should a column of P-values be included with the estimates
@@ -96,11 +102,19 @@ ci.pred( obj, newdata,
parameter vector.}
\item{\ldots}{Parameters passed on to \code{ci.lin}.}
\item{newdata}{Data frame of covariates where prediction is made.}
+ \item{r1,r2}{Estimates of rates in two independent groups, with
+ confidence intervals.}
+ \item{se1,se2}{Standard errors of log-rates in the two groups. If
+ given, it is assumed that \code{r1} and \code{r2} represent
+ log-rates.}
+ \item{log.tr}{Logical, if true, it is assumed that \code{r1} and
+ \code{r2} represent log-rates with confidence intervals.}
}
\value{
\code{ci.lin} returns a matrix with number of rows and row names as
\code{ctr.mat}. The columns are Estimate, Std.Err, z, P, 2.5\% and
- 97.5\%. If \code{vcov=TRUE} a list with components \code{est}, the
+ 97.5\% (or according to the value of \code{alpha}). If
+ \code{vcov=TRUE} a list with components \code{est}, the
desired functional of the parameters and \code{vcov}, the variance
covariance matrix of this, is returned but not printed. If
\code{Exp==TRUE} the confidence intervals for the parameters are
@@ -114,11 +128,11 @@ ci.pred( obj, newdata,
\code{Wald} computes a Wald test for a subset of (possibly linear
combinations of) parameters being equal to the vector of null
- values. The selection of the subset of parameters is the same as for
- \code{ci.lin}. Using the \code{ctr.mat} argument makes it possible to
- do a Wald test for equality of parameters. \code{Wald} returns a named
- numerical vector of length 3, with names \code{Chisq}, \code{d.f.} and
- \code{P}.
+ values as given by \code{H0}. The selection of the subset of
+ parameters is the same as for \code{ci.lin}. Using the \code{ctr.mat}
+ argument makes it possible to do a Wald test for equality of
+ parameters. \code{Wald} returns a named numerical vector of length 3,
+ with names \code{Chisq}, \code{d.f.} and \code{P}.
\code{ci.mat} returns a 2 by 3 matrix with rows \code{c(1,0,0)} and
\code{c(0,-1,1)*1.96}, devised to post-multiply to a p by 2 matrix with
@@ -135,6 +149,12 @@ ci.pred( obj, newdata,
scale, and by default transformed by the inverse link, since the most
common use for this is for multiplicative Poisson or binomial models
with either log or logit link.
+
+ \code{ci.ratio} returns the rate-ratio of two independent set of
+ rates given with confidence intervals or s.e.s. If \code{se1} and
+ \code{se2} are given and \code{log.tr=FALSE} it is assumed that
+ \code{r1} and \code{r2} are rates and \code{se1} and \code{se2} are
+ standard errors of the log-rates.
}
\author{
Bendix Carstensen,
@@ -169,6 +189,10 @@ Wald( mm, subset="f", ctr.mat=CM )
# or alternatively
( CM <- rbind( c(1,-1,0,0), c(0,1,-1,0)) )
Wald( mm, subset="f", ctr.mat=CM )
+
+# Confidnece intervas for ratio of rates
+ci.ratio( cbind(10,8,12.5), cbind(5,4,6.25) )
+ci.ratio( cbind(8,12.5), cbind(4,6.25) )
}
\keyword{models}
\keyword{regression}
diff --git a/man/cutLexis.Rd b/man/cutLexis.Rd
index 4cb3871..2f3640e 100644
--- a/man/cutLexis.Rd
+++ b/man/cutLexis.Rd
@@ -17,10 +17,10 @@
cutLexis( data, cut, timescale = 1,
new.state = nlevels(data$lex.Cst)+1,
new.scale = FALSE,
- split.states = FALSE,
- progressive = FALSE,
- precursor.states = NULL,
- count = FALSE)
+ split.states = FALSE,
+ progressive = FALSE,
+ precursor.states = NULL,
+ count = FALSE )
countLexis( data, cut, timescale = 1 )
}
\arguments{
@@ -116,10 +116,11 @@ countLexis( data, cut, timescale = 1 )
needed for as many times as the person with most events. But also it
is immaterial in what order the cutpoints are entered.
}
-\author{Bendix Carstensen, Steno Diabetes Center, \email{bxc at steno.dk},
+\author{Bendix Carstensen, Steno Diabetes Center, \email{b at bxc.dk},
Martyn Plummer, IARC, \email{plummer at iarc.fr}
}
\seealso{
+ \code{\link{mcutLexis}},
\code{\link{splitLexis}},
\code{\link{Lexis}},
\code{\link{summary.Lexis}},
diff --git a/man/erl.Rd b/man/erl.Rd
index f46b0ff..663d1fa 100644
--- a/man/erl.Rd
+++ b/man/erl.Rd
@@ -56,7 +56,7 @@
\item{immune}{
Logical. Should the years of life lost to the disease be computed
using assumptions that non-diseased individuals are immune to the
- disease (\code{lam}=0) and that their mortality is \code{muW}.
+ disease (\code{lam}=0) and that their mortality is yet still \code{muW}.
}
\item{note}{
Logical. Should a warning of silly assumptions be printed?
@@ -67,7 +67,7 @@
}
\details{
The mortality rates given are supposed to refer to the ages
- \code{age.in+i*int/2}, \code{i=1,2,3,...}.
+ \code{age.in+(i-1/2)*int}, \code{i=1,2,3,...}.
The units in which \code{int} is given must correspond to the units in
which the rates \code{mu}, \code{muW}, \code{muD} and \code{lam} are
@@ -84,17 +84,17 @@
starting at \code{age.in}, and the other that the conditioning ages
given in the argument \code{A} will refer to the ages defined by this.
- The \code{immune} argument is \code{TRUE} whenever the disease
- incidence rates are supplied. If set to \code{FALSE}, the years of
- life lost is computed under the assumption that individuals without
- the disease at a given age are immune to the disease in the sense that
- the disease incidence rate is 0, so transitions to the diseased state
- (with presumably higher mortality rates) are assumed not to occur. This
- is a slightly peculiar assumption (but presumably the most used in the
- epidemiological literature) and the resulting object is therefore
- given an attribute, \code{NOTE}, point this out. The default of the
- \code{surv2} function is to take the possibility of disease into
- account in order to potentially rectify this.}
+ The \code{immune} argument is \code{FALSE} whenever the disease
+ incidence rates are supplied. If set to \code{TRUE}, the years of life
+ lost is computed under the assumption that individuals without the
+ disease at a given age are immune to the disease in the sense that the
+ disease incidence rate is 0, so transitions to the diseased state
+ (with presumably higher mortality rates) are assumed not to
+ occur. This is a slightly peculiar assumption (but presumably the most
+ used in the epidemiological literature) and the resulting object is
+ therefore given an attribute, \code{NOTE}, that point this out. The
+ default of the \code{surv2} function is to take the possibility of
+ disease into account in order to potentially rectify this.}
\value{\code{surv1} and \code{surv2} return a matrix whose first column
is the ages at the ends of the
@@ -125,7 +125,7 @@
difference between the columns is added as a
third column, labeled "YLL".
}
-\author{Bendix Carstensen, \email{bxc at steno.dk}
+\author{Bendix Carstensen, \email{b at bxc.dk}
}
\seealso{
\code{\link{ci.cum}}
@@ -144,7 +144,7 @@ Lc <- cutLexis( Lx, cut = Lx$doins-Lx$dob,
precursor.states = "DM" )
summary( Lc )
# Split in small age intervals
-Sc <- splitLexis( Lc, breaks=seq(0,120,1) )
+Sc <- splitLexis( Lc, breaks=seq(0,120,2) )
summary( Sc )
# Overview of object
diff --git a/man/fit.mult.Rd b/man/fit.mult.Rd
index a4e08d9..5a573ed 100644
--- a/man/fit.mult.Rd
+++ b/man/fit.mult.Rd
@@ -46,7 +46,7 @@
}
\author{
Martyn Plummer, \email{plummer at iarc.fr},
- Bendix Carstensen, \email{bxc at steno.dk} }
+ Bendix Carstensen, \email{b at bxc.dk} }
\seealso{
\code{\link{Icens}}
\code{\link{fit.add}}
diff --git a/man/foreign.Lexis.Rd b/man/foreign.Lexis.Rd
index 3bd4aac..7932e8b 100644
--- a/man/foreign.Lexis.Rd
+++ b/man/foreign.Lexis.Rd
@@ -55,7 +55,7 @@ msdata(obj, ...)
class \code{etm}.
}
\author{
-Bendix Carstensen, \email{bxc at steno.dk}, \url{http://BendixCarstensen.com}
+Bendix Carstensen, \email{b at bxc.dk}, \url{http://BendixCarstensen.com}
}
\examples{
data(DMlate)
diff --git a/man/gen.exp.Rd b/man/gen.exp.Rd
index e24fc74..9c6ff3d 100644
--- a/man/gen.exp.Rd
+++ b/man/gen.exp.Rd
@@ -28,13 +28,13 @@ gen.exp(purchase, id = "id", dop = "dop", amt = "amt", dpt = "dpt",
units used for \code{dpt} is assumed to be units of \code{amt} per
units of \code{dop}.}
\item{id}{Character. Name of the id variable in the data frame.}
- \item{dop}{Character. Name of the date of purchase variable in the data frame.}
- \item{amt}{Character. Name of the amount purchased variable in the data frame.}
- \item{dpt}{Character. Name of the dose-per-time variable in the data frame.}
- \item{fu}{Data frame with follow-up period for each person, the person
+ \item{dop}{Character. Name of the \code{d}ate \code{o}f \code{p}urchase variable in the data frame.}
+ \item{amt}{Character. Name of the \code{am}oun\code{t} purchased variable in the data frame.}
+ \item{dpt}{Character. Name of the \code{d}ose-\code{p}er-\code{t}ime variable in the data frame.}
+ \item{fu}{Data frame with \code{f}ollow-\code{u}p period for each person, the person
id variable must have the same name as in the \code{purchase} data frame.}
- \item{doe}{Character. Name of the date of entry variable.}
- \item{dox}{Character. Name of the date of exit variable.}
+ \item{doe}{Character. Name of the \code{d}ate \code{o}f \code{e}ntry variable.}
+ \item{dox}{Character. Name of the \code{d}ate \code{o}f e\code{x}it variable.}
\item{use.dpt}{Logical: should we use information on dose per time.}
\item{breaks}{Numerical vector of dates at which the time since
first exposure, cumulative dose etc. are computed.}
@@ -97,17 +97,17 @@ gen.exp(purchase, id = "id", dop = "dop", amt = "amt", dpt = "dpt",
from possibly the first interval for each person, this will assume values in the set of
the values in \code{breaks}.}
\item{\code{Y}}{the length of interval.}
- \item{\code{tfi}}{time from first initiation of drug.}
- \item{\code{tfc}}{time from latest cessation of drug.}
- \item{\code{cdur}}{cumulative time on the drug.}
- \item{\code{cdos}}{cumulative dose.}
+ \item{\code{tfi}}{\code{t}ime \code{f}rom first \code{i}nitiation of drug.}
+ \item{\code{tfc}}{\code{t}ime \code{f}rom latest \code{c}essation of drug.}
+ \item{\code{cdur}}{\code{c}umulative \code{dur}ation of the drug.}
+ \item{\code{cdos}}{\code{c}umulative \code{dos}e.}
\item{\code{ldos}}{suffixed with one value per element in
\code{lags}, the latter giving the cumulative doses \code{lags} before
\code{dof}.}
}
}
-\author{Bendix Carstensen, \email{bxc at steno.dk}. The development of
+\author{Bendix Carstensen, \email{b at bxc.dk}. The development of
this function was supported partly through a grant from the EFSD
(European Foundation for the Study of Diabetes), ""}
diff --git a/man/lgrep.Rd b/man/lgrep.Rd
new file mode 100644
index 0000000..13a0d0d
--- /dev/null
+++ b/man/lgrep.Rd
@@ -0,0 +1,39 @@
+\name{lgrep}
+\alias{fgrep}
+\alias{ngrep}
+\alias{lgrep}
+\title{
+Convenience versions of grep
+}
+\description{
+Often you want the elements of a vector (or its names or levels) that
+meet a certain pattern. But \code{grep} only gives you the position, so
+these functions are designed to give you that.
+}
+\usage{
+fgrep( pattern, x, ... )
+ngrep( pattern, x, ... )
+lgrep( pattern, x, ... )
+}
+\arguments{
+ \item{pattern}{Pattern searched for.}
+ \item{x}{Object where \code{pattern} is searched. Or in whose \code{names}
+ or \code{levels} attributes \code{pattern} is sought.}
+ \item{...}{Arguments passed on to \code{\link[base]{grep}}.}
+}
+\value{Elements of the input \code{x} (\code{fgrep}) or its names
+ attribute (\code{ngrep}) or levels attribute (\code{lgrep}).
+}
+\author{Bendix Carstensen, \email{b at bxc.dk}, \url{www.bxc.dk} }
+\seealso{\code{\link{grep}}}
+\examples{
+ff <- factor( ll <- paste( sample( letters[1:3], 20, replace=TRUE ),
+ sample( letters[1:3], 20, replace=TRUE ), sep="" ) )
+ff
+fgrep( "a", ff )
+fgrep( "a", ll )
+ngrep( "a", ff )
+lgrep( "a", ff )
+lgrep( "a", ff, invert=TRUE )
+}
+\keyword{ manip }
diff --git a/man/mcutLexis.Rd b/man/mcutLexis.Rd
new file mode 100644
index 0000000..6e85ca9
--- /dev/null
+++ b/man/mcutLexis.Rd
@@ -0,0 +1,96 @@
+\name{mcutLexis}
+\alias{mcutLexis}
+\title{
+Cut follow-up at multiple event dates and keep track of order of events
+}
+\description{
+A generalization of \code{\link{cutLexis}} to the case where different
+events may occur in any order.
+}
+\usage{
+mcutLexis( L0, timescale = 1, wh,
+ new.states = NULL,
+ precursor.states = NULL,
+ seq.states = TRUE,
+ new.scales = NULL,
+ ties.resolve = FALSE )
+}
+\arguments{
+ \item{L0}{A Lexis object.}
+ \item{timescale}{Which time scale do the variables in \code{L0[,wh]}
+ refer to. Can be character or integer.}
+ \item{wh}{Which variables contain the event dates. Character or
+ integer vector}
+ \item{new.states}{Names of the events forming new states. If
+ \code{NULL} equal to the variable names from \code{wh}.}
+ \item{precursor.states}{Which states are precursor states. See
+ \code{\link{cutLexis}} for definition of precursor states.}
+ \item{seq.states}{Should the sequence of events be kept track of? That
+ is, should A-B be considered different from B-A. If \code{FALSE}, the
+ state with combined preceding events A and B will be called A+B.}
+ \item{new.scales}{Should we construct new time scales indicating the
+ time since each of the event occurrences.}
+ \item{ties.resolve}{Should tied event times be resolved by adding
+ random noise to tied event dates. If \code{FALSE} the function will
+ not accept that two events occur at the same time for a person
+ (ties). If \code{TRUE} a random quantity in the range
+ \code{c(-1,1)/100} will be added to all event times in all records
+ with at least one tie. If numeric a random quantity in the range
+ \code{c(-1,1)*ties.resolve} will be added to all event times in all
+ records with at least one tie.}
+}
+\value{A \code{\link{Lexis}} object with extra states created by
+ occurrence of a number of intermediate events.
+}
+\author{
+Bendix Carstensen, \url{http://BendixCarstensen.com}
+}
+\seealso{
+\code{\link{cutLexis}}, \code{\link{Lexis}}, \code{\link{splitLexis}}
+}
+\examples{
+# A dataframe of times
+set.seed(563248)
+dd <- data.frame( id = 1:10,
+ doN = round(runif(10,-30, 0),1),
+ doE = round(runif(10, 0,20),1),
+ doX = round(runif(10, 50,60),1),
+ doD = round(runif(10, 50,60),1),
+ # these are the event times
+ doA = c(NA,20,NA,27,35,NA,52, 5,43,80),
+ doB = c(25,NA,37,40,NA,NA,15,23,36,61) )
+
+# set up a Lexis object with time from entry to death/exit
+Lx <- Lexis( entry = list(time=doE,
+ age=doE-doN),
+ exit = list(time=pmin(doX,doD)),
+ exit.status = factor(doD<doX,labels=c("OK","D")),
+ data = dd )
+summary( Lx )
+
+# cut the follow-up at dates doA and doB
+L2 <- mcutLexis( Lx, "time", wh=c("doA","doB"),
+ new.states = c("A","B"),
+ precursor.states = "OK",
+ seq.states = TRUE,
+ new.scales = c("tfA","tfB") )
+summary( L2 )
+L2
+
+# show the states
+boxes( L2, boxpos=list(x=c(10,50,50,90,50,90),
+ y=c(50,90,50,90,10,10)),
+ show.R=FALSE, show.BE=TRUE )
+
+
+L3 <- mcutLexis( Lx, "time", wh=c("doA","doB"),
+ new.states = c("A","B"),
+ precursor.states = "OK",
+ seq.states = FALSE,
+ new.scales = c("tfA","tfB") )
+summary( L3 )
+boxes( L3, boxpos=list(x=c(10,50,50,90,50),
+ y=c(50,90,50,50,10)),
+ show.R=FALSE, show.BE=TRUE )
+}
+\keyword{survival}
diff --git a/man/ncut.Rd b/man/ncut.Rd
index 47af77e..028d52c 100644
--- a/man/ncut.Rd
+++ b/man/ncut.Rd
@@ -26,7 +26,7 @@ The function uses the base function \code{findInterval}.
A numerical vector of the same length as \code{x}.
}
\author{
- Bendix Carstensen, Steno Diabetes Center, \email{bxc at steno.dk},
+ Bendix Carstensen, Steno Diabetes Center, \email{b at bxc.dk},
\url{http://BendixCarstensen.com}, with essential input
from Martyn Plummer, IARC.
}
diff --git a/man/nice.Rd b/man/nice.Rd
index fef2ea3..c103096 100644
--- a/man/nice.Rd
+++ b/man/nice.Rd
@@ -17,7 +17,7 @@ nice(x, log = F, lpos = c(1, 2, 5), ...)
\item{\dots}{Arguments passed on to \code{pretty} if \code{log}=FALSE}
}
\value{A vector of breakpoints.}
-\author{Bendix Carstensen, \email{bxc at steno.dk}, \url{http://BendixCarstensen.com}}
+\author{Bendix Carstensen, \email{b at bxc.dk}, \url{http://BendixCarstensen.com}}
\seealso{pretty}
\examples{
nice( exp( rnorm( 100 ) ), log=TRUE )
diff --git a/man/plotCIF.Rd b/man/plotCIF.Rd
new file mode 100644
index 0000000..2ec05ad
--- /dev/null
+++ b/man/plotCIF.Rd
@@ -0,0 +1,122 @@
+\name{plotCIF}
+\alias{plotCIF}
+\alias{stackedCIF}
+
+\title{Plotting Aalen-Johansen curves for competing events
+}
+\description{Function \code{plotCIF} plots, for one or more groups,
+the cumulative incidence curves for a selected event out of two or more
+competing events. Function \code{stackedCIF} plots, for one group or population,
+the cumulative incidence curves for two or more competing events such that
+the cumulative incidences are stacked upon each other. The CIFs are
+are estimated by the Aalen-Johansen method.
+}
+\usage{
+## S3 method for class 'survfit'
+ plotCIF( x, event = 1,
+ xlab = "Time",
+ ylab = "Cumulative incidence",
+ ylim = c(0, 1),
+ lty = NULL,
+ col = NULL, ... )
+
+## S3 method for class 'survfit'
+stackedCIF( x, group = 1,
+ colour = NULL,
+ ylim = c(0, 1),
+ xlab = "Time",
+ ylab = "Cumulative incidence", ... )
+}
+\arguments{
+ \item{x}{An object of class \code{\link{survfit}}, the \code{type}
+ of \code{event} in \code{Surv()} being "\code{mstate}";
+ the first level of the event factor represents censoring and the
+ remaining ones the alternative competing events.
+}
+ \item{event}{Determines the event for which the cumulative incidence curve is plotted by \code{plotCIF}.
+}
+\item{group}{An integer showing the selected level of a possible grouping factor
+ appearing in the model formula in \code{survfit} when plotting by \code{stackedCIF}
+}
+\item{col}{A vector specifying the plotting color(s) of the curve(s)
+ for the different groups in \code{\link{plotCIF}}-- default: all "black".
+}
+\item{colour}{A vector indicating the colours to be
+ used for shading the areas pertinent to the separate outcomes in \code{\link{stackedCIF}} --
+ default: all \code{"white"}.
+}
+\item{xlab}{Label for the $x$-axis.
+}
+ \item{ylab}{Label for the $y$-axis.
+}
+ \item{ylim}{Limits of the $y$-axis.
+}
+
+ \item{lty}{A vector specifying the line type(s) of the curve(s) for the different groups --
+ default: all 1 (=solid).
+}
+ \item{\dots}{Further graphical parameters to be passed.
+}
+}
+\details{
+The order in which the curves with \code{\link{stackedCIF}} are piled
+upon each other is the same as the ordering of the values or levels of
+the competing events in the pertinent event variable. The ordering can
+be changed by permuting the levels as desired using function
+\code{Relevel}, after which \code{survfit} is called with the relevelled
+\code{event} variable in \code{Surv()}
+}
+\value{No value is returned but a plot is produced as a side-effect.
+}
+\references{Putter, H., Fiocco, M., Geskus, R.B. (2007).
+Tutorial in biostatistics: competing risks and multi-state models.
+Statistics in Medicine, 26: 2389--2430.
+}
+\author{Esa L{\"a}{\"a}r{\"a}, \email{esa.laara at oulu.fi}
+}
+\note{
+Aalen-Johansen curves for competing events in several groups can also
+be plotted by function \code{\link{plot.survfit}} of the survival
+library as well as by some functions in other packages covering analysis
+of time-to-event data.}
+\seealso{
+\code{\link{survfit}}, \code{\link{plot}}, \code{\link{plot.survfit}}.
+}
+\examples{
+library(survival) # requires version 2.39-4 or later
+head(mgus1)
+# Aalen-Johansen estimates of CIF are plotted by sex for two
+# competing events: (1) progression (pcm), and (2) death, in
+# a cohort of patients with monoclonal gammopathy.
+
+# The data are actually covering transitions from pcm to death, too,
+# for those entering the state of pcm. Such patients have two rows
+# in the data frame, and in their 2nd row the 'start' time is
+# the time to pcm (in days).
+
+# In our analysis we shall only include those time intervals with value 0
+# for variable 'start'. Thus, the relevant follow-up time is represented
+# by variable 'stop' (days). For convenience, days are converted to years.
+
+fitCI <- survfit(Surv(stop/365.25, event, type="mstate") ~ sex,
+ data= subset(mgus1, start==0) )
+par(mfrow=c(1,2))
+plotCIF(fitCI, event = 1, col = c("red", "blue"),
+ main = "Progression", xlab="Time (years)" )
+text( 38, 0.15, "Men", pos = 2)
+text( 38, 0.4, "Women", pos = 2)
+plotCIF(fitCI, event = 2, col = c("red", "blue"),
+ main = "Death", xlab="Time (years)" )
+text( 38, 0.8, "Men", pos = 2)
+text( 38, 0.5, "Women", pos = 2)
+
+par(mfrow=c(1,2))
+stackedCIF(fitCI, group = 1, colour = c("gray80", "gray90"),
+ main = "Women", xlab="Time (years)" )
+text( 36, 0.15, "PCM", pos = 2)
+text( 36, 0.6, "Death", pos = 2)
+stackedCIF(fitCI, group = 2, colour = c("gray80", "gray90"),
+ main = "Men", xlab="Time (years)" )
+text( 39, 0.10, "PCM", pos = 2)
+text( 39, 0.6, "Death", pos = 2)
+}
diff --git a/man/plotEst.Rd b/man/plotEst.Rd
index b9dd732..1a99bc5 100644
--- a/man/plotEst.Rd
+++ b/man/plotEst.Rd
@@ -87,7 +87,7 @@ pointsEst( ests, y = dim(ests)[1]:1, pch = 16, cex = 1, lwd = 2,
}
\author{
Bendix Carstensen,
- \email{bxc at steno.dk},
+ \email{b at bxc.dk},
\url{http://BendixCarstensen.com}}
\seealso{
ci.lin
diff --git a/man/splitLexis.Rd b/man/splitLexis.Rd
index 0a1e8e0..b870baf 100644
--- a/man/splitLexis.Rd
+++ b/man/splitLexis.Rd
@@ -85,5 +85,9 @@ tapply( status(x2,"exit")==1, list( timeBand(x2,"age","left"),
tapply( dur(x2), list( timeBand(x2,"age","left"),
timeBand(x2,"per","left") ), sum )
}
-\seealso{\code{\link{timeBand}}, \code{\link{cutLexis}}, \code{\link{summary.Lexis}}}
+\seealso{
+ \code{\link{timeBand}},
+ \code{\link{cutLexis}},
+ \code{\link{mcutLexis}},
+ \code{\link{summary.Lexis}}}
\keyword{manip}
diff --git a/man/stack.Lexis.Rd b/man/stack.Lexis.Rd
index e79406e..1228f82 100644
--- a/man/stack.Lexis.Rd
+++ b/man/stack.Lexis.Rd
@@ -49,7 +49,7 @@ a factor with levels made up of combinations of the levels of
joined by a "\code{->}".}
\author{
-Bendix Carstensen, \email{bxc at steno.dk}, \url{http://BendixCarstensen.com}
+Bendix Carstensen, \email{b at bxc.dk}, \url{http://BendixCarstensen.com}
}
\examples{
data(DMlate)
diff --git a/man/summary.Lexis.rd b/man/summary.Lexis.rd
index c85beb9..296b790 100644
--- a/man/summary.Lexis.rd
+++ b/man/summary.Lexis.rd
@@ -45,7 +45,9 @@
If the argument \code{Rates} is FALSE (the default), then only the
first component of the list is returned.
}
-\author{Bendix Carstensen, \email{bxc at steno.dk}}
+\author{
+ Bendix Carstensen, \url{http://BendixCarstensen.com}
+ }
\examples{
data( nickel )
# Lung cancer deaths and other deaths are coded 1 and 2
diff --git a/man/transform.Lexis.Rd b/man/transform.Lexis.Rd
index ec4f090..a13bb6b 100644
--- a/man/transform.Lexis.Rd
+++ b/man/transform.Lexis.Rd
@@ -5,6 +5,8 @@
\alias{factorize}
\alias{factorize.Lexis}
\alias{levels.Lexis}
+\alias{order.Lexis}
+\alias{sort.Lexis}
\title{Transform a Lexis (or stacked.Lexis) objects}
\description{
Modify a Lexis object.
@@ -15,6 +17,8 @@
\method{Relevel}{Lexis}( x, states, print = TRUE, \dots )
\method{levels}{Lexis}( x )
\method{factorize}{Lexis}( x, states, print = TRUE, \dots )
+%\method{order}{Lexis}( x, \dots )
+%\method{sort}{Lexis}( x, \dots )
\method{transform}{stacked.Lexis}( `_data`, \dots)
}
\arguments{
@@ -43,7 +47,13 @@
not passed to the function, the returned object have levels of
\code{lex.Cst}, \code{lex.Xst} (and for \code{stacked.Lexis} objects
\code{lex.Tr}) shaved down to the actually occurring values.
- }
+
+ \code{order} returns the order of the rows in a Lexis object to sort
+ it by (\code{lex.id},\code{timeScales[x]}).
+
+ \code{sort} returns the Lexis object sorted by
+ (\code{lex.id},\code{timeScales[x]}).
+}
\value{
A transformed \code{Lexis} object.
diff --git a/src/chinv2.c b/src/chinv2.c
index dce6d6c..c19f353 100644
--- a/src/chinv2.c
+++ b/src/chinv2.c
@@ -13,7 +13,9 @@
** Copied from the survival package by Terry Therneau, version 2.35-7
*/
-void chinv2(double **matrix , int n)
+#include <R_ext/Visibility.h>
+
+void attribute_hidden chinv2(double **matrix , int n)
{
register double temp;
register int i,j,k;
diff --git a/src/cholesky2.c b/src/cholesky2.c
index 81bcf62..6736639 100644
--- a/src/cholesky2.c
+++ b/src/cholesky2.c
@@ -22,7 +22,9 @@
* Copied from the survival package by Terry Therneau, version 2.35-7
*/
-int cholesky2(double **matrix, int n, double toler)
+#include <R_ext/Visibility.h>
+
+int attribute_hidden cholesky2(double **matrix, int n, double toler)
{
double temp;
int i,j,k;
diff --git a/src/chsolve2.c b/src/chsolve2.c
index 97052b4..8c90241 100644
--- a/src/chsolve2.c
+++ b/src/chsolve2.c
@@ -14,7 +14,9 @@
* Copied from the survival package by Terry Therneau, version 2.35-7
*/
-void chsolve2(double **matrix, int n, double *y)
+#include <R_ext/Visibility.h>
+
+void attribute_hidden chsolve2(double **matrix, int n, double *y)
{
register int i,j;
register double temp;
diff --git a/src/clogit.c b/src/clogit.c
index 92f3f8a..216b23a 100644
--- a/src/clogit.c
+++ b/src/clogit.c
@@ -7,9 +7,9 @@
These are specially adapted to work with information matrices that
are not of full rank.
*/
-int cholesky2(double **matrix, int m, double toler);
-void chsolve2(double **matrix, int m, double *y);
-void chinv2(double **matrix, int m);
+extern int cholesky2(double **matrix, int m, double toler);
+extern void chsolve2(double **matrix, int m, double *y);
+extern void chinv2(double **matrix, int m);
/*
Efficient calculation of the conditional log likelihood, along with
diff --git a/src/init.c b/src/init.c
new file mode 100644
index 0000000..07a3015
--- /dev/null
+++ b/src/init.c
@@ -0,0 +1,18 @@
+#include <R.h>
+#include <Rinternals.h>
+#include <stdlib.h> // for NULL
+#include <R_ext/Rdynload.h>
+
+extern SEXP clogit(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP);
+
+static const R_CallMethodDef CallEntries[] = {
+ {"C_clogit", (DL_FUNC) &clogit, 7},
+ {NULL, NULL, 0}
+};
+
+void R_init_Epi(DllInfo *dll)
+{
+ R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
+ R_useDynamicSymbols(dll, FALSE);
+ R_forceSymbols(dll, TRUE);
+}
diff --git a/vignettes/index.html b/vignettes/index.html
index 5d67758..9da1be2 100644
--- a/vignettes/index.html
+++ b/vignettes/index.html
@@ -10,6 +10,9 @@
<li> <a href="simLexis.pdf">
Simulation in multistate models with multiple timescales
</a>, and the corresponding <a href="simLexis.R"> R-code </a>.
+ <li> <a href="yll.pdf">
+ Calculation of years of life lost,
+ </a>, and the corresponding <a href="yll.R"> R-code </a>.
</ul>
Here is the website for the <a href="http://BendixCarstensen.com/Epi/">Epi package</a>.
</html>
diff --git a/vignettes/simLexis.rnw b/vignettes/simLexis.rnw
index aa12586..e5c9abf 100644
--- a/vignettes/simLexis.rnw
+++ b/vignettes/simLexis.rnw
@@ -366,21 +366,27 @@ knots so that the number of events is the same between each pair of
knots, but only half of this beyond each of the boundary knots,
whereas for the timescales \texttt{DMdur} and \texttt{tIns} where we
have observation from a well-defined 0, we put knots at 0 and place the
-remaining knots so that the number of events is the same between them
-and beyond the last:
+remaining knots so that the number of events is the same between each
+pair of knots as well as outside the boundary knots.
<<knots>>=
nk <- 5
-( ai.kn <- with( subset(Si,lex.Xst=="Ins"),
+( ai.kn <- with( subset(Si,lex.Xst=="Ins" & lex.Cst!=lex.Xst ),
quantile( Age+lex.dur , probs=(1:nk-0.5)/nk ) ) )
( ad.kn <- with( subset(Si,lex.Xst=="Dead"),
quantile( Age+lex.dur , probs=(1:nk-0.5)/nk ) ) )
-( di.kn <- with( subset(Si,lex.Xst=="Ins"),
+( di.kn <- with( subset(Si,lex.Xst=="Ins" & lex.Cst!=lex.Xst ),
c(0,quantile( DMdur+lex.dur, probs=(1:(nk-1))/nk ) )) )
( dd.kn <- with( subset(Si,lex.Xst=="Dead"),
c(0,quantile( DMdur+lex.dur, probs=(1:(nk-1))/nk ) )) )
( ti.kn <- with( subset(Si,lex.Xst=="Dead(Ins)"),
c(0,quantile( t.Ins+lex.dur, probs=(1:(nk-1))/nk ) )) )
@ %
+Note that when we tease out the event records for transition to
+\emph{transient} states (in this case ``Ins'', that is
+verb|lex.Xst=="Ins"|), we should add \verb|lex.Cst!=lex.Xst|, to
+include only transition records and avoiding including records of
+sojourn time in the transient state.
+
We then fit Poisson models to transition rates, using the wrapper
\texttt{Ns} from the \texttt{Epi} package to simplify the
specification of the rates:
@@ -1025,7 +1031,7 @@ follows:
This section explains the actually existing code for
\texttt{simLexis}, as it is in the current version of \texttt{Epi}.
<<echo=FALSE,results=hide>>=
-options( keep.source=TRUE )
+source( "../R/simLexis.R", keep.source=TRUE )
@ %
The function \texttt{simLexis} takes a \texttt{Lexis} object as
input. This defines the initial state(s) and times of the start for a
@@ -1088,7 +1094,7 @@ and objects of class \texttt{function} is made.
The dataset on which this prediction is done has
\texttt{length(time.pts)} rows per person.
<<>>=
-Epi:::simX
+simX
@ %
As we see, \texttt{simX} calls \texttt{sim1} which simulates the
transition for one person.
@@ -1101,7 +1107,7 @@ The predicted cumulative intensities are fed, person by person, to
appended to the \texttt{nd} data frame. This way we have simulated
\emph{one} transition (time and state) for each person:
<<>>=
-Epi:::sim1
+sim1
@ %
The \texttt{sim1} function uses \texttt{lint} to do linear interpolation.
@@ -1114,7 +1120,7 @@ a custom made linear interpolator that does the interpolation
exploiting the fact the the \texttt{ci} is non-decreasing and
\texttt{tt} is strictly monotonously increasing:
<<>>=
-Epi:::lint
+lint
@
\subsection{\texttt{get.next}}
@@ -1125,7 +1131,7 @@ defined as time since entry to one of these states be initialized to 0
before a call to \texttt{simX} is made for these persons. This
accomplished by the function \texttt{get.next}:
<<>>=
-Epi:::get.next
+get.next
@
\subsection{\texttt{chop.lex}}
@@ -1135,7 +1141,7 @@ after \emph{each} new entry to a transient state. In order to groom
the output data we use \texttt{chop.lex} to censor all persons at the
same designated time after \emph{initial} entry.
<<>>=
-Epi:::chop.lex
+chop.lex
@
\section{Probabilities from simulated \texttt{Lexis} objects}
diff --git a/vignettes/yll.rnw b/vignettes/yll.rnw
new file mode 100644
index 0000000..5da0094
--- /dev/null
+++ b/vignettes/yll.rnw
@@ -0,0 +1,641 @@
+\SweaveOpts{results=verbatim,keep.source=TRUE,include=FALSE,eps=FALSE}
+%\VignetteIndexEntry{Years of life lost: simLexis}
+\documentclass[a4paper,twoside,12pt]{report}
+
+% ----------------------------------------------------------------------
+% General information for the title page and the page headings
+\newcommand{\Title}{Years of Life Lost (YLL)
+ to disease\\Diabetes in DK as example}
+\newcommand{\Tit}{YLL}
+\newcommand{\Version}{Version 1.2}
+\newcommand{\Dates}{February 2017}
+\newcommand{\Where}{SDC}
+\newcommand{\Homepage}{\url{http://bendixcarstensen.com/Epi}}
+\newcommand{\Faculty}{\begin{tabular}{rl}
+Bendix Carstensen
+ & Steno Diabetes Center, Gentofte, Denmark\\
+ & {\small \& Department of Biostatistics,
+ University of Copenhagen} \\
+ & \texttt{b at bxc.dk}\\
+ & \url{http://BendixCarstensen.com} \\[1em]
+ \end{tabular}}
+% Packages
+\usepackage[utf8]{inputenc}
+\usepackage[T1]{fontenc}
+\usepackage[english]{babel}
+\usepackage[font=it,labelfont=normalfont]{caption}
+\usepackage[colorlinks,urlcolor=blue,linkcolor=red,,citecolor=Maroon]{hyperref}
+\usepackage[ae,hyper]{Rd}
+\usepackage[dvipsnames]{xcolor}
+\usepackage[super]{nth}
+\usepackage[noae]{Sweave}
+\usepackage{makeidx,floatflt,amsmath,amsfonts,amsbsy,enumitem,dcolumn,needspace}
+\usepackage{ifthen,calc,eso-pic,everyshi}
+\usepackage{booktabs,longtable,rotating,graphicx,subfig}
+\usepackage{pdfpages,verbatim,fancyhdr,datetime,%
+afterpage}
+\usepackage[abspath]{currfile}
+% \usepackage{times}
+\renewcommand{\textfraction}{0.0}
+\renewcommand{\topfraction}{1.0}
+\renewcommand{\bottomfraction}{1.0}
+\renewcommand{\floatpagefraction}{0.9}
+% \usepackage{mslapa}
+\definecolor{blaa}{RGB}{99,99,255}
+\DeclareGraphicsExtensions{.png,.pdf,.jpg}
+% Make the Sweave output nicer
+\DefineVerbatimEnvironment{Sinput}{Verbatim}{fontsize=\small,fontshape=sl,formatcom=\color{Blue}}
+\DefineVerbatimEnvironment{Soutput}{Verbatim}{fontsize=\small,formatcom=\color{Maroon},xleftmargin=0em}
+\DefineVerbatimEnvironment{Scode}{Verbatim}{fontsize=\small}
+\fvset{listparameters={\setlength{\topsep}{-0.1ex}}}
+\renewenvironment{Schunk}%
+{\renewcommand{\baselinestretch}{0.85} \vspace{\topsep}}%
+{\renewcommand{\baselinestretch}{1.00} \vspace{\topsep}}
+\providecommand{\ptxt}[1]{\Pp\left\{\text{#1}\right\}}
+\providecommand{\dif}{{\,\mathrm d}}
+\DeclareMathOperator{\YLL}{YLL}
+\DeclareMathOperator{\Pp}{P}
+
+%----------------------------------------------------------------------
+% Set up layout of pages
+\oddsidemargin 1mm
+\evensidemargin 1mm
+\topmargin -10mm
+\headheight 8mm
+\headsep 5mm
+\textheight 240mm
+\textwidth 165mm
+%\footheight 5mm
+\footskip 15mm
+\renewcommand{\topfraction}{0.9}
+\renewcommand{\bottomfraction}{0.9}
+\renewcommand{\textfraction}{0.1}
+\renewcommand{\floatpagefraction}{0.9}
+\renewcommand{\headrulewidth}{0.1pt}
+\setcounter{secnumdepth}{4}
+% \setcounter{tocdepth}{2}
+
+%----------------------------------------------------------------------
+% How to insert a figure in a .rnw file
+\newcommand{\rwpre}{./graph/gr}
+\newcommand{\insfig}[3]{
+\begin{figure}[h]
+ \centering
+ \includegraphics[width=#2\textwidth]{\rwpre-#1}
+ \caption{#3}
+ \label{fig:#1}
+% \afterpage{\clearpage}
+\end{figure}}
+
+%----------------------------------------------------------------------
+% Here is the document starting with the titlepage
+\begin{document}
+
+%----------------------------------------------------------------------
+% The title page
+\setcounter{page}{1}
+\pagenumbering{roman}
+\pagestyle{plain}
+\thispagestyle{empty}
+% \vspace*{0.05\textheight}
+\flushright
+% The blank below here is necessary in order not to muck up the
+% linespacing in title if it has more than 2 lines
+{\Huge \bfseries \Title
+
+}\ \\[-1.5ex]
+\noindent\textcolor{blaa}{\rule[-1ex]{\textwidth}{5pt}}\\[2.5ex]
+\large
+\Where \\
+\Dates \\
+\Homepage \\
+\Version \\[1em]
+\normalsize
+Compiled \today,\ \currenttime\\
+from: \texttt{\currfileabspath}\\[1em]
+% \input{wordcount}
+\normalsize
+\vfill
+\Faculty
+% End of titlepage
+% \newpage
+
+%----------------------------------------------------------------------
+% Table of contents
+\tableofcontents
+
+%----------------------------------------------------------------------
+% General text layout
+\raggedright
+\parindent 1em
+\parskip 0ex
+\cleardoublepage
+
+%----------------------------------------------------------------------
+% General page style
+\pagenumbering{arabic}
+\setcounter{page}{1}
+\pagestyle{fancy}
+\renewcommand{\chaptermark}[1]{\markboth{\textsl{#1}}{}}
+\renewcommand{\sectionmark}[1]{\markright{\thesection\ \textsl{#1}}{}}
+\fancyhead[EL]{\bf \thepage \quad \rm \leftmark}
+\fancyhead[ER,OL]{\sl \Tit}
+\fancyhead[OR]{\rm \rightmark \quad \bf \thepage}
+\fancyfoot{}
+
+<<echo=FALSE>>=
+options( width=90,
+ SweaveHooks=list( fig=function()
+ par(mar=c(3,3,1,1),mgp=c(3,1,0)/1.6,las=1,bty="n") ) )
+@ %
+\renewcommand{\rwpre}{./yll}
+
+%----------------------------------------------------------------------
+% Here comes the substance part
+\chapter{Theory and technicalities}
+
+This vignette for the \texttt{Epi} package describes the
+probabilistic/demographic background for and technical implementation
+of the \texttt{erl} and \texttt{yll} functions that computes the
+expected residual life time and years of life lost in an illness-death
+model.
+
+\section{Years of life lost (YLL)}
+
+\ldots to diabetes or any other disease for that matter.
+
+The general concept in calculation of ``years lost to\ldots'' is the
+comparison of the expected lifetime between two groups of persons; one
+with and one without disease (in this example DM). The expected
+lifetime is the area under the survival curve, so basically the
+exercise requires that two survival curves that are deemed relevant be
+available.
+
+The years of life lost is therefore just the area between the survival
+curves for those ``Well'', $S_W(t)$, and for those ``Diseased'',
+$S_D(t)$:
+\[
+ \YLL = \int_0^\infty S_W(t) - S_D(t) \dif t
+\]
+The time $t$ could of course be age, but it could also be ``time after
+age 50'' and the survival curves compared would then be survival
+curves \emph{conditional} on survival till age 50, and the YLL would
+be the years of life lost for a 50-year old person with diabetes.
+
+If we are referring to the expected lifetime we will more precisely use
+the label expected residual lifetime, ERL.
+
+\section{Constructing the survival curves}
+
+YLL can be computed in two different ways, depending on the way the
+survival curve and hence the expected lifetime of a person
+\emph{without} diabetes is computed:
+\begin{itemize}
+\item Assume that the ``Well'' persons are \emph{immune} to disease
+ --- using only the non-DM mortality rates throughout for calculation
+ of expected life time.
+\item Assume that the ``Well'' persons \emph{can} acquire the disease and
+ thereby see an increased mortality, thus involving all three rates
+ shown in figure \ref{fig:states}.
+\end{itemize}
+The former gives a higher YLL because the comparison is to persons
+assumed immune to DM (and yet with the same mortality as non-immune
+prior to diagnosis), the latter gives a more realistic picture of the
+comparison of group of persons with and without diabetes at a given
+age that can be interpreted in the real world.
+
+The differences can be illustrated by figure \ref{fig:states}; the
+immune approach corresponds to an assumption of $\lambda(t)=0$ in the
+calculation of the survival curve for a person in the ``Well'' state.
+
+Calculation of the survival of a diseased person already in the ``DM''
+state is unaffected by assumptions about $\lambda$.
+
+<<states,fig=TRUE,height=5,width=5,echo=FALSE>>=
+library( Epi )
+TM <- matrix(NA,4,4)
+rownames(TM) <-
+colnames(TM) <- c("Well","DM","Dead","Dead(DM)")
+TM[1,2:3] <- TM[2,4] <- 1
+zz <- boxes( TM, boxpos=list(x=c(20,80,20,80),y=c(80,80,20,20)), wm=1.5, hm=4 )
+@
+<<states,fig=TRUE,height=4,width=5,echo=FALSE>>=
+zz$Arrowtext <- c( expression(lambda),
+ expression(mu[W]),
+ expression(mu[D][M]) )
+boxes( zz )
+@ %
+\insfig{states}{0.7}{Illness-death model describing diabetes incidence
+ and -mortality.}
+
+\subsection{Total mortality --- a shortcut?}
+
+A practical crude shortcut could be to compare the ERL in the diabetic
+population to the ERL for the \emph{entire} population (that is use
+the total mortality ignoring diabetes status).
+
+Note however that this approach also counts the mortality of persons
+that acquired the disease earlier, thus making the comparison
+population on average more ill than the population we aim at, namely
+those well at a given time, which only then become more gradually ill.
+
+How large these effects are can however be empirically explored, as we
+shall do later.
+
+\subsection{Disease duration}
+
+In the exposition above there is no explicit provision for the effect of
+disease duration, but if we were able to devise mortality rates for
+any combination of age and duration, this could be taken into account.
+
+There are however severe limitations in this as we in principle would
+want to have duration effects as long as the age-effects --- in
+principle for all $(a,d)$ where $d\leq A$, where $A$ is the age at
+which we condition. So even if we were only to compute ERL from
+age, say, 40 we would still need duration effects up to 60 years
+(namely to age 100).
+
+The incorporation of duration effects is in principle trivial from a
+computational point of view, but we would be forced to entertain
+models predicting duration effects way beyond what is actually
+observed disease duration in any practical case.
+
+\subsection{Computing integrals}
+
+The practical calculations of survival curves, ERL and YLL involves
+calculation of (cumulative) integrals of rates and functions of these
+as we shall see below. This is easy if we have a closed form
+expression of the function, so its value may be computed at any time
+point --- this will be the case if we model rates by smooth parametric
+functions.
+
+Computing the (cumulative) integral of a function is done as follows:
+\begin{itemize}
+\item Compute the value of the function (mortality rate for example)
+ at the midpoints of a sequence of narrow equidistant intervals ---
+ for example one- or three month intervals of age, say.
+\item Take the cumulative sum of these values multiplied by the
+ interval length --- this will be a very close approximation to the
+ cumulative integral evaluated at the end of each interval.
+\item If the intervals are really small (like 1/100 year), the
+ distinction between the value at the middle and at the end of each
+ interval becomes irrelevant.
+\end{itemize}
+Note that in the above it is assumed that the rates are given in units
+corresponding to the interval length --- or more precisely, as the
+cumulative rates over the interval.
+
+\section{Survival functions in the illness-death model}
+
+The survival functions for persons in the ``Well'' state can be
+computed under two fundamentally different scenarios, depending on
+whether persons in the ``Well'' state are assumed to be immune to the
+disease ($\lambda(a)=0$) or not.
+
+\subsection{Immune approach}
+
+In this case both survival functions for person in the two states are
+the usual simple transformation of the cumulative mortality rates:
+\[
+ S_W(a) = \exp\left(-\int_0^a\!\!\mu_W(u) \dif u \right), \qquad
+ S_D(a) = \exp\left(-\int_0^a\!\!\mu_D(u) \dif u \right)
+\]
+
+\subsubsection{Conditional survival functions}
+
+If we want the \emph{conditional} survival functions given survival to
+age $A$, say, they are just:
+\[
+ S_W(a|A) = S_W(a)/S_W(A), \qquad S_D(a|A) = S_D(a)/S_D(A)
+\]
+
+\subsection{Non-immune approach}
+
+For a diseased person, the survival function in this states is the same
+as above, but the survival function for a person without disease (at
+age 0) is (see figure \ref{fig:states}):
+\[
+S(a) = \ptxt{Well}\!(a) + \ptxt{DM}\!(a)
+\]
+In the appendix of the paper \cite{Carstensen.2008c} is an indication
+of how to compute the probability of being in any of the four states
+shown in figure \ref{fig:states}, which I shall repeat here:
+
+In terms of the rates, the probability of being in the ``Well'' box is
+simply the probability of escaping both death (at a rate of $\mu_W(a)$)
+and diabetes (at a rate of $\lambda(a)$):
+\[
+ \ptxt{Well}(a) = \exp\left(\!-\int_0^a\!\!\mu_W(u)+\lambda(u) \right) \dif u
+\]
+The probability of being alive with diabetes at age $a$, is computed given that
+ diabetes occurred at age $s$ ($s<a$) and then integrated over $s$ from $0$
+ to $a$:
+\begin{align*}
+ \ptxt{DM}(a) = \int_0^a\!\! & \ptxt{survive to $s$, DM diagnosed at $s$} \\
+ & \times \ptxt{survive with DM from $s$ to $a$} \dif s \\
+ = \int_0^a\!\! & \lambda(s)
+ \exp\left(\!-\int_0^s\!\!\mu_W(u)+\lambda(u) \dif u \right) \\
+ & \times \exp\left(\!-\int_s^a\!\!\mu_D(u) \dif u \right) \dif s
+\end{align*}
+Sometimes we will use a version where the mortality among diabetes
+patients depend both on age $a$ and duration of diabetes, $d$,
+$\mu_D(a,d)$, in which case we get:
+\begin{align*}
+ \ptxt{DM}(a) = \int_0^a \! & \lambda(s)
+ \exp\left(-\int_0^s\!\mu_W(u)+\lambda(u) \dif u \right) \\
+ & \times \exp\left(-\int_s^a\!\mu_D(u,u-s) \dif u \right) \dif s
+\end{align*}
+because the integration variable $u$ is the age-scale and the second
+integral refers to mortality among persons diagnosed at age $s$, that
+is, with duration $u-s$ at age $u$.
+
+The option of using duration-dependent mortality rates among diseased
+individuals is not implemented yet.
+
+\subsubsection{Conditional survival functions}
+
+Unlike the immune approach, the conditional survival function in the
+more realistic case is not just a ratio of the unconditional to the
+value at the conditioning age, $A$, say. This would amount to
+conditioning on being merely \emph{alive} at age $A$, but what we want
+is to condition on being in the ``Well'' state at age $A$.
+
+The formulae for the conditional probabilities of being either in
+``Well'' or ``DM'', given being in ``Well'' at age $A$ are basically
+replicates of the unconditional, albeit with changes in integration
+limits:
+\begin{align*}
+\ptxt{Well|Well at $A$}(a) &= \exp\left(-\int_A^a \! \mu_W(u)+\lambda(u) \right) \dif u \\
+ \ptxt{DM|Well at $A$}(a) &= \int_A^a \! \lambda(s)
+ \exp\left(-\int_A^s\!\mu_W(u)+\lambda(u) \dif u \right) \\
+ & \qquad \times \exp\left(-\int_s^a\!\mu_D(u,u-s) \dif u \right) \dif s
+\end{align*}
+The calculation of these conditional survival functions is implemented
+but not allowing for duration-dependence. Thus it is only implemented
+assuming $\mu_D(a,d)=\mu_D(a)$.
+
+\chapter{Analyses for DM in Denmark}
+
+The rates we use as basis for the following calculations are derived
+from the NDR, where we have omitted the blood-glucose criteria,
+because there is compelling evidence that these have quite a low
+specificity (particularly in the younger ages among women), and do
+not substantially contribute to the sensitivity.
+
+As noted above the calculations of YLL requires access to
+(age-specific) rates of incidence of DM and mortality for persons with
+and without DM.
+
+\section{Modeling mortality and incidence data}
+
+We read in the dataset of DM and population mortality and incidence, \texttt{DMepi}:
+<<>>=
+data( DMepi )
+@ %
+The dataset \texttt{DMepi} contains diabetes events, deaths and
+person-years for persons without diabetes and deaths and person-years
+for persons with diabetes:
+<<>>=
+str( DMepi )
+head( DMepi )
+@ %
+For each combination of sex, age, period and date of birth in 1 year
+age groups, we have the person-years in the ``Well'' (\texttt{Y.nD})
+and the ``DM'' (\texttt{Y.DM}) states, as well as the number of deaths
+from these (\texttt{D.nD}, \texttt{D.DM}) and the number of incident
+diabetes cases from the ``Well'' state (\texttt{X}).
+
+In order to compute the years of life lost to diabetes and how this
+has changed over time, we fit models for the mortality and incidence
+of both groups (and of course, separately for men and women). The
+models we use will be age-period-cohort models \cite{Carstensen.2007a}
+providing estimated mortality rates for ages 0--99 and dates
+1.1.1996--1.1.2016.
+
+First we transform the age and period variables to reflect the mean
+age and period in each of the Lexis triangles. We also compute the
+total number of deaths and amount of risk time, as we are going to
+model the total mortality as well. Finally we restrict the dataset to
+ages over 30 only:
+<<>>=
+DMepi <- transform( subset( DMepi, A>30 ),
+ D.T = D.nD + D.DM,
+ Y.T = Y.nD + Y.DM )
+head(DMepi)
+@ %
+With the correct age and period coding in the Lexis triangles, we fit
+models for the mortalities and incidences. Note that we for
+comparative purposes also fit a model for the \emph{total} mortality,
+ignoring the
+<<>>=
+# Knots used in all models
+( a.kn <- seq(40,95,,6) )
+( p.kn <- seq(1997,2015,,4) )
+( c.kn <- seq(1910,1976,,6) )
+# Check the number of events between knots
+ae <- xtabs( cbind(D.nD,D.DM,X) ~ cut(A,c(30,a.kn,Inf)) + sex, data=DMepi )
+ftable( addmargins(ae,1), col.vars=3:2 )
+pe <- xtabs( cbind(D.nD,D.DM,X) ~ cut(P,c(1990,p.kn,Inf)) + sex, data=DMepi )
+ftable( addmargins(pe,1), col.vars=3:2 )
+ce <- xtabs( cbind(D.nD,D.DM,X) ~ cut(P-A,c(-Inf,c.kn,Inf)) + sex, data=DMepi )
+ftable( addmargins(ce,1), col.vars=3:2 )
+# Fit an APC-model for all transitions, seperately for men and women
+mW.m <- glm( D.nD ~ -1 + Ns(A ,knots=a.kn,int=TRUE) +
+ Ns( P,knots=p.kn,ref=2005) +
+ Ns(P-A,knots=c.kn,ref=1950),
+ offset = log(Y.nD),
+ family = poisson,
+ data = subset( DMepi, sex=="M" ) )
+mD.m <- update( mW.m, D.DM ~ . , offset=log(Y.DM) )
+mT.m <- update( mW.m, D.T ~ . , offset=log(Y.T ) )
+lW.m <- update( mW.m, X ~ . )
+# Model for women
+mW.f <- update( mW.m, data = subset( DMepi, sex=="F" ) )
+mD.f <- update( mD.m, data = subset( DMepi, sex=="F" ) )
+mT.f <- update( mT.m, data = subset( DMepi, sex=="F" ) )
+lW.f <- update( lW.m, data = subset( DMepi, sex=="F" ) )
+@ %
+
+\section{Residual life time and years lost to DM}
+
+We now collect the estimated years of life lost classified by method
+(immune assumption or not), sex, age and calendar time:
+<<>>=
+a.ref <- 30:90
+p.ref <- 1996:2016
+aYLL <- NArray( list( type = c("Imm","Tot","Sus"),
+ sex = levels( DMepi$sex ),
+ age = a.ref,
+ date = p.ref ) )
+str( aYLL )
+system.time(
+for( ip in p.ref )
+ {
+ nd <- data.frame( A = seq(30,90,0.2)+0.1,
+ P = ip,
+ Y.nD = 1,
+ Y.DM = 1,
+ Y.T = 1 )
+ muW.m <- ci.pred( mW.m, nd )[,1]
+ muD.m <- ci.pred( mD.m, nd )[,1]
+ muT.m <- ci.pred( mT.m, nd )[,1]
+ lam.m <- ci.pred( lW.m, nd )[,1]
+ muW.f <- ci.pred( mW.f, nd )[,1]
+ muD.f <- ci.pred( mD.f, nd )[,1]
+ muT.f <- ci.pred( mT.f, nd )[,1]
+ lam.f <- ci.pred( lW.f, nd )[,1]
+ aYLL["Imm","M",,paste(ip)] <- yll( int=0.2, muW.m, muD.m, lam=NULL,
+ A=a.ref, age.in=30, note=FALSE )[-1]
+ aYLL["Imm","F",,paste(ip)] <- yll( int=0.2, muW.f, muD.f, lam=NULL,
+ A=a.ref, age.in=30, note=FALSE )[-1]
+ aYLL["Tot","M",,paste(ip)] <- yll( int=0.2, muT.m, muD.m, lam=NULL,
+ A=a.ref, age.in=30, note=FALSE )[-1]
+ aYLL["Tot","F",,paste(ip)] <- yll( int=0.2, muT.f, muD.f, lam=NULL,
+ A=a.ref, age.in=30, note=FALSE )[-1]
+ aYLL["Sus","M",,paste(ip)] <- yll( int=0.2, muW.m, muD.m, lam=lam.m,
+ A=a.ref, age.in=30, note=FALSE )[-1]
+ aYLL["Sus","F",,paste(ip)] <- yll( int=0.2, muW.f, muD.f, lam=lam.f,
+ A=a.ref, age.in=30, note=FALSE )[-1]
+ } )
+round( ftable( aYLL[,,seq(1,61,10),], col.vars=c(3,2) ), 1 )
+@ %
+We now have the relevant points for the graph showing YLL to diabetes
+for men and women by age, and calendar year, both under the immunity
+and susceptibility models for the calculation of YLL.
+<<imm,fig=TRUE,width=8,height=5>>=
+plyll <- function(wh){
+par( mfrow=c(1,2), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, bty="n", las=1 )
+
+matplot( a.ref, aYLL[wh,"M",,],
+ type="l", lty=1, col="blue", lwd=1:2,
+ ylim=c(0,12), xlab="Age",
+ ylab="Years lost to DM", yaxs="i" )
+abline(v=50,h=1:10,col=gray(0.7))
+text( 90, 11, "Men", col="blue", adj=1 )
+text( 40, aYLL[wh,"M","40","1996"], "1996", adj=c(0,0), col="blue" )
+text( 43, aYLL[wh,"M","44","2016"], "2016", adj=c(1,1), col="blue" )
+
+matplot( a.ref, aYLL[wh,"F",,],
+ type="l", lty=1, col="red", lwd=1:2,
+ ylim=c(0,12), xlab="Age",
+ ylab="Years lost to DM", yaxs="i" )
+abline(v=50,h=1:10,col=gray(0.7))
+text( 90, 11, "Women", col="red", adj=1 )
+text( 40, aYLL[wh,"F","40","1996"], "1996", adj=c(0,0), col="red" )
+text( 43, aYLL[wh,"F","44","2016"], "2016", adj=c(1,1), col="red" )
+}
+plyll("Imm")
+@ %
+<<tot,fig=TRUE,width=8,height=5>>=
+plyll("Tot")
+@ %
+<<sus,fig=TRUE,width=8,height=5>>=
+plyll("Sus")
+@ %
+\begin{figure}[h]
+ \centering
+ \includegraphics[width=\textwidth]{yll-imm}
+ \caption{Years of life lost to DM: the difference in expected
+ residual life time at different ages between persons with and
+ without diabetes, assuming the persons without diabetes at a given
+ age remain free from diabetes (immunity assumption --- not
+ reasonable). The lines refer to date of evaluation; the top lines
+ refer to 1.1.1996 the bottom ones to 1.1.2016. Blue curves are
+ men, red women.}
+ \label{fig:imm}
+\end{figure}
+
+\begin{figure}[h]
+ \centering
+ \includegraphics[width=\textwidth]{yll-sus}
+ \caption{Years of life lost to DM: the difference in expected
+ residual life time at different ages between persons with and
+ without diabetes, allowing the persons without diabetes at a given
+ to contract diabetes and thus be subject to higher mortality. The
+ lines refer to date of evaluation; the top lines refer to 1.1.1996
+ the bottom ones to 1.1.2016. Blue curves are men, red women.}
+ \label{fig:sus}
+\end{figure}
+
+\begin{figure}[h]
+ \centering
+ \includegraphics[width=\textwidth]{yll-tot}
+ \caption{Years of life lost to DM: the difference in expected
+ residual life time at different ages between persons with and
+ without diabetes. Allowance for susceptibility is approximated by
+ using the total population mortality instead of non-DM
+ mortality. The lines refer to date of evaluation; the top lines
+ refer to 1.1.1996 the bottom ones to 1.1.2016. Blue curves are
+ men, red women.}
+ \label{fig:tot}
+\end{figure}
+
+From figure \ref{fig:sus} we see that for men aged 50 the years lost to
+diabetes has decreased from a bit over 8 to a bit less than 6 years,
+and for women from 8.5 to 5 years; so a greater improvement for women.
+
+\chapter{Practical implementation}
+
+We have devised functions that wraps these formulae up for practical
+use.
+
+\section{Function definitions}
+
+<<echo=FALSE,results=hide,eval=TRUE>>=
+source( "../R/erl.R", keep.source=TRUE )
+@
+When using the functions it is assumed that the functions $\mu_W$,
+$\mu_D$ and $\lambda$ are given as vectors corresponding to
+equidistantly (usually tightly) spaced ages from 0 to $K$ where K is
+the age where everyone can safely be assumed dead.
+
+\texttt{surv1} is a simple function that computes the survival
+function from a vector of mortality rates, and optionally the
+conditional survival given being alive at prespecified ages:
+<<>>=
+surv1
+@ %
+\texttt{erl1} basically just expands the result of \texttt{surv1} with
+a column of expected residual life times:
+<<>>=
+erl1
+@ %
+We also define a function, \texttt{surv2}, that computes the survival
+function for a non-diseased person that may become diseased with rate
+\texttt{lam} and after that die at a rate of \texttt{muD}
+(corresponding to the formulae above). This is the sane way of
+handling years of life lost to a particular illness:
+<<>>=
+surv2
+@ %
+Finally we devised a function using these to compute the expected
+residual lifetime at select ages:
+<<>>=
+erl
+@ %
+\ldots and a wrapper for this if we only want the years of life lost
+returned:
+<<>>=
+yll
+@ %
+
+\bibliographystyle{plain}
+
+\begin{thebibliography}{1}
+
+\bibitem{Carstensen.2007a}
+B~Carstensen.
+\newblock Age-{P}eriod-{C}ohort models for the {L}exis diagram.
+\newblock {\em Statistics in Medicine}, 26(15):3018--3045, July 2007.
+
+\bibitem{Carstensen.2008c}
+B.~Carstensen, J.K. Kristensen, P.~Ottosen, and K.~Borch-Johnsen.
+\newblock The {D}anish {N}ational {D}iabetes {R}egister: {T}rends in incidence,
+ prevalence and mortality.
+\newblock {\em Diabetologia}, 51:2187--2196, 2008.
+
+\end{thebibliography}
+
+\addcontentsline{toc}{chapter}{References}
+
+\end{document}
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/r-cran-epi.git
More information about the debian-med-commit
mailing list