Binary files /tmp/tmpaTZoIU/2Rrh_8wnSC/r-cran-epi-2.0/build/vignette.rds and /tmp/tmpaTZoIU/0lMC5lKUYD/r-cran-epi-2.7/build/vignette.rds differ diff -Nru r-cran-epi-2.0/CHANGES r-cran-epi-2.7/CHANGES --- r-cran-epi-2.0/CHANGES 2016-01-06 18:14:41.000000000 +0000 +++ r-cran-epi-2.7/CHANGES 2016-10-02 08:56:12.000000000 +0000 @@ -1,3 +1,57 @@ +Changes in 2.6 + +o Added function erl computing Expected Residual Lifetime in an + illness-death model added, together with companions surv1, surv2, + erl1 and yll (Years of Life Lost). + +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. + +o rm.tr added; removes transitions from a Lexis object + +o boxes.MS no longer isssues a warning when show.BE is set to TRUE. + +o LCa.fit rewritten and expanded to encompass both age-period and + age-cohort multiplicative interactions. + +o APC.LCa added, fits all possible Lee-Carter type models and + APC-models. boxes.APC.LCa plots the relationship between models + including the residual deviances, and optionally places boxes to + provide overview of best fitting models. + +Changes in 2.4 + +o Ns updated with the possibility of clamping effects to have 0 slope + beyond the outer knots, see argument "fixsl". + +Changes in 2.3 + +o Cplot (usually called from rateplot) now checks if age- and + period-groupings are of the same length, and tells you if they are + not, instead of just plotting (almost) nothing. + +o Grooming of LCa functions, and in particular the documentation. + +o Update of apc.fit so that also Y^2/D (the observed information about + the rate) and Y (person-time) is allowed as weight when defining the + inner product inducing orthogonality between linear and non-linear + effects. + +Changes in 2.2 + +o LCa.fit added: Fits Lee-Carter models with smooth age and time + effects to rate-data. print, summary, plot and predict methods also + supplied. + +Changes in 2.1 + +o The show.BE="nz" feature in boxes.Lexis is now documented + +o "[.Lexis" is now exported and works... + Changes in 2.0 o cbind, rbind and "[" methods for Lexis objects have been added diff -Nru r-cran-epi-2.0/debian/changelog r-cran-epi-2.7/debian/changelog --- r-cran-epi-2.0/debian/changelog 2016-04-27 19:17:46.000000000 +0000 +++ r-cran-epi-2.7/debian/changelog 2016-10-28 22:01:47.000000000 +0000 @@ -1,3 +1,12 @@ +r-cran-epi (2.7-1) unstable; urgency=medium + + * New upstream version + * Convert to dh-r + * Canonical homepage for CRAN + * New Build-Depends: r-cran-numderiv, r-cran-data.table + + -- Andreas Tille Sat, 29 Oct 2016 00:01:47 +0200 + r-cran-epi (2.0-2) unstable; urgency=medium * Fix missing Depends: r-cran-plyr diff -Nru r-cran-epi-2.0/debian/control r-cran-epi-2.7/debian/control --- r-cran-epi-2.0/debian/control 2016-04-27 19:17:41.000000000 +0000 +++ r-cran-epi-2.7/debian/control 2016-10-28 22:01:47.000000000 +0000 @@ -4,25 +4,26 @@ Section: gnu-r Priority: optional Build-Depends: debhelper (>= 9), - cdbs, + dh-r, r-base-dev, r-cran-mass, r-cran-cmprsk, r-cran-etm, - r-cran-plyr + r-cran-plyr, + r-cran-numderiv, + r-cran-data.table Standards-Version: 3.9.8 Vcs-Browser: https://anonscm.debian.org/viewvc/debian-med/trunk/packages/R/r-cran-epi/trunk/ Vcs-Svn: svn://anonscm.debian.org/debian-med/trunk/packages/R/r-cran-epi/trunk/ -Homepage: http://staff.pubhealth.ku.dk/~bxc/Epi/ +Homepage: https://cran.r-project.org/package=Epi Package: r-cran-epi Architecture: any Depends: ${shlibs:Depends}, - ${R:Depends}, - r-cran-mass, - r-cran-cmprsk, - r-cran-etm, - r-cran-plyr + ${misc:Depends}, + ${R:Depends} +Recommends: ${R:Recommends} +Suggests: ${R:Suggests} Description: GNU R epidemiological analysis Functions for demographic and epidemiological analysis in the Lexis diagram, i.e. register and cohort follow-up data, including interval censored data and diff -Nru r-cran-epi-2.0/debian/copyright r-cran-epi-2.7/debian/copyright --- r-cran-epi-2.0/debian/copyright 2015-09-24 09:12:07.000000000 +0000 +++ r-cran-epi-2.7/debian/copyright 2016-10-28 22:01:03.000000000 +0000 @@ -1,14 +1,14 @@ -Format: http://www.debian.org/doc/packaging-manuals/copyright-format/1.0/ +Format: https://www.debian.org/doc/packaging-manuals/copyright-format/1.0/ Upstream-Name: Epi Upstream-Contact: Bendix Carstensen -Source: http://cran.r-project.org/src/contrib/ +Source: https://cran.r-project.org/package=Epi Files: * -Copyright: 2008-2012 Bendix Carstensen, Martyn Plummer, Esa Laara, Michael Hills et. al. +Copyright: 2008-2016 Bendix Carstensen, Martyn Plummer, Esa Laara, Michael Hills et. al. License: GPL-2+ Files: debian/* -Copyright: 2008-2015 Andreas Tille +Copyright: 2008-2016 Andreas Tille License: GPL-2+ License: GPL-2+ diff -Nru r-cran-epi-2.0/debian/rules r-cran-epi-2.7/debian/rules --- r-cran-epi-2.0/debian/rules 2014-06-27 14:03:42.000000000 +0000 +++ r-cran-epi-2.7/debian/rules 2016-10-28 22:01:27.000000000 +0000 @@ -1,8 +1,6 @@ #!/usr/bin/make -f -# -*- makefile -*- -# debian/rules file for the Debian/GNU Linux r-cran-epi package -# Copyright 2008 by Andreas Tille - -include /usr/share/R/debian/r-cran.mk DEB_COMPRESS_EXCLUDE_ALL := .pdf + +%: + dh $@ --buildsystem R diff -Nru r-cran-epi-2.0/DESCRIPTION r-cran-epi-2.7/DESCRIPTION --- r-cran-epi-2.0/DESCRIPTION 2016-01-06 19:29:29.000000000 +0000 +++ r-cran-epi-2.7/DESCRIPTION 2016-10-04 12:35:49.000000000 +0000 @@ -1,6 +1,6 @@ Package: Epi -Version: 2.0 -Date: 2016-01-06 +Version: 2.7 +Date: 2016-10-04 Title: A Package for Statistical Analysis in Epidemiology Authors@R: c(person("Bendix", "Carstensen", role = c("aut", "cre"), email = "bxc@steno.dk"), @@ -9,24 +9,26 @@ person("Esa", "Laara", role = "ctb"), person("Michael", "Hills", role = "ctb")) Depends: R (>= 3.0.0), utils -Imports: cmprsk, etm, splines, MASS, survival, plyr +Imports: cmprsk, etm, splines, MASS, survival, plyr, Matrix, numDeriv, + data.table Suggests: mstate, nlme, lme4 Description: Functions for demographic and epidemiological analysis in the Lexis diagram, i.e. register and cohort follow-up data, in particular representation, manipulation and simulation of multistate data - the Lexis suite of functions, which includes interfaces to - mstate, etm and cmprsk packages. - Also contains functions for Age-Period-Cohort modeling and a - function for interval censored data and some useful functions for - tabulation and plotting, as well some epidemiological datasets. + '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. License: GPL-2 URL: http://BendixCarstensen.com/Epi/ NeedsCompilation: yes -Packaged: 2016-01-06 18:15:40 UTC; bendix +Packaged: 2016-10-03 19:59:16 UTC; bendix Author: Bendix Carstensen [aut, cre], Martyn Plummer [aut], Esa Laara [ctb], Michael Hills [ctb] Maintainer: Bendix Carstensen Repository: CRAN -Date/Publication: 2016-01-06 20:29:29 +Date/Publication: 2016-10-04 14:35:49 diff -Nru r-cran-epi-2.0/inst/CITATION r-cran-epi-2.7/inst/CITATION --- r-cran-epi-2.0/inst/CITATION 2015-05-27 15:27:08.000000000 +0000 +++ r-cran-epi-2.7/inst/CITATION 2016-10-01 16:53:04.000000000 +0000 @@ -13,14 +13,14 @@ as.person("Michael Hills")), year = year, note = note, - url = "http://CRAN.R-project.org/package=Epi", + url = "https://CRAN.R-project.org/package=Epi", textVersion = paste("Bendix Carstensen, Martyn Plummer, Esa Laara, Michael Hills", sprintf("(%s).", year), "Epi: A Package for Statistical Analysis in Epidemiology.", paste(note, ".", sep = ""), - "URL http://CRAN.R-project.org/package=Epi") + "URL https://CRAN.R-project.org/package=Epi") ) citEntry(entry = "Article", Binary files /tmp/tmpaTZoIU/2Rrh_8wnSC/r-cran-epi-2.0/inst/doc/Follow-up.pdf and /tmp/tmpaTZoIU/0lMC5lKUYD/r-cran-epi-2.7/inst/doc/Follow-up.pdf differ diff -Nru r-cran-epi-2.0/inst/doc/Follow-up.R r-cran-epi-2.7/inst/doc/Follow-up.R --- r-cran-epi-2.0/inst/doc/Follow-up.R 2016-01-04 13:51:22.000000000 +0000 +++ r-cran-epi-2.7/inst/doc/Follow-up.R 2016-10-03 19:59:16.000000000 +0000 @@ -1,159 +1,159 @@ -### R code from vignette source 'Follow-up.rnw' - -################################################### -### code chunk number 1: Follow-up.rnw:65-67 -################################################### -library(Epi) -print( sessionInfo(), l=F ) - - -################################################### -### code chunk number 2: Follow-up.rnw:146-153 -################################################### -data( nickel ) -nicL <- Lexis( entry = list( per=agein+dob, - age=agein, - tfh=agein-age1st ), - exit = list( age=ageout ), - exit.status = ( icd %in% c(162,163) )*1, - data = nickel ) - - -################################################### -### code chunk number 3: Follow-up.rnw:163-166 -################################################### -str( nickel ) -str( nicL ) -head( nicL ) - - -################################################### -### code chunk number 4: Follow-up.rnw:175-176 -################################################### -summary( nicL ) - - -################################################### -### code chunk number 5: nicL1 -################################################### -plot( nicL ) - - -################################################### -### code chunk number 6: nicL2 -################################################### -par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 ) -plot( nicL, 1:2, lwd=1, col=c("blue","red")[(nicL$exp>0)+1], - grid=TRUE, lty.grid=1, col.grid=gray(0.7), - xlim=1900+c(0,90), xaxs="i", - ylim= 10+c(0,90), yaxs="i", las=1 ) -points( nicL, 1:2, pch=c(NA,3)[nicL$lex.Xst+1], - col="lightgray", lwd=3, cex=1.5 ) -points( nicL, 1:2, pch=c(NA,3)[nicL$lex.Xst+1], - col=c("blue","red")[(nicL$exp>0)+1], lwd=1, cex=1.5 ) - - -################################################### -### code chunk number 7: Follow-up.rnw:229-232 -################################################### -nicS1 <- splitLexis( nicL, "age", breaks=seq(0,100,10) ) -summary( nicL ) -summary( nicS1 ) - - -################################################### -### code chunk number 8: Follow-up.rnw:239-240 -################################################### -round( subset( nicS1, id %in% 8:10 ), 2 ) - - -################################################### -### code chunk number 9: Follow-up.rnw:245-247 -################################################### -nicS2 <- splitLexis( nicS1, "tfh", breaks=c(0,1,5,10,20,30,100) ) -round( subset( nicS2, id %in% 8:10 ), 2 ) - - -################################################### -### code chunk number 10: Follow-up.rnw:253-258 -################################################### -timeBand( nicS2, "age", "middle" )[1:20] -# For nice printing and column labelling use the data.frame() function: -data.frame( nicS2[,c("id","lex.id","per","age","tfh","lex.dur")], - mid.age=timeBand( nicS2, "age", "middle" ), - mid.tfh=timeBand( nicS2, "tfh", "middle" ) )[1:20,] - - -################################################### -### code chunk number 11: Follow-up.rnw:276-281 -################################################### -subset( nicL, id %in% 8:10 ) -agehi <- nicL$age1st + 50 / nicL$exposure -nicC <- cutLexis( data=nicL, cut=agehi, timescale="age", - new.state=2, precursor.states=0 ) -subset( nicC, id %in% 8:10 ) - - -################################################### -### code chunk number 12: Follow-up.rnw:287-292 -################################################### -subset( nicS2, id %in% 8:10 ) -agehi <- nicS2$age1st + 50 / nicS2$exposure -nicS2C <- cutLexis( data=nicS2, cut=agehi, timescale="age", - new.state=2, precursor.states=0 ) -subset( nicS2C, id %in% 8:10 ) - - -################################################### -### code chunk number 13: Follow-up.rnw:312-321 -################################################### -data( nickel ) -nicL <- Lexis( entry = list( per=agein+dob, - age=agein, - tfh=agein-age1st ), - exit = list( age=ageout ), - exit.status = ( icd > 0 ) + ( icd %in% c(162,163) ), - data = nickel ) -summary( nicL ) -subset( nicL, id %in% 8:10 ) - - -################################################### -### code chunk number 14: Follow-up.rnw:325-333 -################################################### -nicL <- Lexis( entry = list( per=agein+dob, - age=agein, - tfh=agein-age1st ), - exit = list( age=ageout ), - exit.status = ( icd > 0 ) + ( icd %in% c(162,163) ), - data = nickel, - states = c("Alive","D.oth","D.lung") ) -summary( nicL ) - - -################################################### -### code chunk number 15: Follow-up.rnw:345-353 -################################################### -nicL$agehi <- nicL$age1st + 50 / nicL$exposure -nicC <- cutLexis( data = nicL, - cut = nicL$agehi, - timescale = "age", - new.state = "HiExp", - precursor.states = "Alive" ) -subset( nicC, id %in% 8:10 ) -summary( nicC, scale=1000 ) - - -################################################### -### code chunk number 16: Follow-up.rnw:371-379 -################################################### -nicC <- cutLexis( data = nicL, - cut = nicL$agehi, - timescale = "age", - new.state = "Hi", - split.states=TRUE, new.scale=TRUE, - precursor.states = "Alive" ) -subset( nicC, id %in% 8:10 ) -summary( nicC, scale=1000 ) - - +### R code from vignette source 'Follow-up.rnw' + +################################################### +### code chunk number 1: Follow-up.rnw:65-67 +################################################### +library(Epi) +print( sessionInfo(), l=F ) + + +################################################### +### code chunk number 2: Follow-up.rnw:146-153 +################################################### +data( nickel ) +nicL <- Lexis( entry = list( per=agein+dob, + age=agein, + tfh=agein-age1st ), + exit = list( age=ageout ), + exit.status = ( icd %in% c(162,163) )*1, + data = nickel ) + + +################################################### +### code chunk number 3: Follow-up.rnw:163-166 +################################################### +str( nickel ) +str( nicL ) +head( nicL ) + + +################################################### +### code chunk number 4: Follow-up.rnw:175-176 +################################################### +summary( nicL ) + + +################################################### +### code chunk number 5: nicL1 +################################################### +plot( nicL ) + + +################################################### +### code chunk number 6: nicL2 +################################################### +par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 ) +plot( nicL, 1:2, lwd=1, col=c("blue","red")[(nicL$exp>0)+1], + grid=TRUE, lty.grid=1, col.grid=gray(0.7), + xlim=1900+c(0,90), xaxs="i", + ylim= 10+c(0,90), yaxs="i", las=1 ) +points( nicL, 1:2, pch=c(NA,3)[nicL$lex.Xst+1], + col="lightgray", lwd=3, cex=1.5 ) +points( nicL, 1:2, pch=c(NA,3)[nicL$lex.Xst+1], + col=c("blue","red")[(nicL$exp>0)+1], lwd=1, cex=1.5 ) + + +################################################### +### code chunk number 7: Follow-up.rnw:229-232 +################################################### +nicS1 <- splitLexis( nicL, "age", breaks=seq(0,100,10) ) +summary( nicL ) +summary( nicS1 ) + + +################################################### +### code chunk number 8: Follow-up.rnw:239-240 +################################################### +round( subset( nicS1, id %in% 8:10 ), 2 ) + + +################################################### +### code chunk number 9: Follow-up.rnw:245-247 +################################################### +nicS2 <- splitLexis( nicS1, "tfh", breaks=c(0,1,5,10,20,30,100) ) +round( subset( nicS2, id %in% 8:10 ), 2 ) + + +################################################### +### code chunk number 10: Follow-up.rnw:253-258 +################################################### +timeBand( nicS2, "age", "middle" )[1:20] +# For nice printing and column labelling use the data.frame() function: +data.frame( nicS2[,c("id","lex.id","per","age","tfh","lex.dur")], + mid.age=timeBand( nicS2, "age", "middle" ), + mid.tfh=timeBand( nicS2, "tfh", "middle" ) )[1:20,] + + +################################################### +### code chunk number 11: Follow-up.rnw:276-281 +################################################### +subset( nicL, id %in% 8:10 ) +agehi <- nicL$age1st + 50 / nicL$exposure +nicC <- cutLexis( data=nicL, cut=agehi, timescale="age", + new.state=2, precursor.states=0 ) +subset( nicC, id %in% 8:10 ) + + +################################################### +### code chunk number 12: Follow-up.rnw:287-292 +################################################### +subset( nicS2, id %in% 8:10 ) +agehi <- nicS2$age1st + 50 / nicS2$exposure +nicS2C <- cutLexis( data=nicS2, cut=agehi, timescale="age", + new.state=2, precursor.states=0 ) +subset( nicS2C, id %in% 8:10 ) + + +################################################### +### code chunk number 13: Follow-up.rnw:312-321 +################################################### +data( nickel ) +nicL <- Lexis( entry = list( per=agein+dob, + age=agein, + tfh=agein-age1st ), + exit = list( age=ageout ), + exit.status = ( icd > 0 ) + ( icd %in% c(162,163) ), + data = nickel ) +summary( nicL ) +subset( nicL, id %in% 8:10 ) + + +################################################### +### code chunk number 14: Follow-up.rnw:325-333 +################################################### +nicL <- Lexis( entry = list( per=agein+dob, + age=agein, + tfh=agein-age1st ), + exit = list( age=ageout ), + exit.status = ( icd > 0 ) + ( icd %in% c(162,163) ), + data = nickel, + states = c("Alive","D.oth","D.lung") ) +summary( nicL ) + + +################################################### +### code chunk number 15: Follow-up.rnw:345-353 +################################################### +nicL$agehi <- nicL$age1st + 50 / nicL$exposure +nicC <- cutLexis( data = nicL, + cut = nicL$agehi, + timescale = "age", + new.state = "HiExp", + precursor.states = "Alive" ) +subset( nicC, id %in% 8:10 ) +summary( nicC, scale=1000 ) + + +################################################### +### code chunk number 16: Follow-up.rnw:371-379 +################################################### +nicC <- cutLexis( data = nicL, + cut = nicL$agehi, + timescale = "age", + new.state = "Hi", + split.states=TRUE, new.scale=TRUE, + precursor.states = "Alive" ) +subset( nicC, id %in% 8:10 ) +summary( nicC, scale=1000 ) + + diff -Nru r-cran-epi-2.0/inst/doc/index.html r-cran-epi-2.7/inst/doc/index.html --- r-cran-epi-2.0/inst/doc/index.html 2016-01-01 15:21:26.000000000 +0000 +++ r-cran-epi-2.7/inst/doc/index.html 2016-10-01 10:08:21.000000000 +0000 @@ -7,9 +7,9 @@
  • Follow-up data with the Epi package , and the corresponding R-code . -
  • +
  • Simulation in multistate models with multiple timescales - , and the corresponding R-code . + , and the corresponding R-code . Here is the website for the Epi package. Binary files /tmp/tmpaTZoIU/2Rrh_8wnSC/r-cran-epi-2.0/inst/doc/sim-Lexis.pdf and /tmp/tmpaTZoIU/0lMC5lKUYD/r-cran-epi-2.7/inst/doc/sim-Lexis.pdf differ Binary files /tmp/tmpaTZoIU/2Rrh_8wnSC/r-cran-epi-2.0/inst/doc/simLexis.pdf and /tmp/tmpaTZoIU/0lMC5lKUYD/r-cran-epi-2.7/inst/doc/simLexis.pdf differ diff -Nru r-cran-epi-2.0/inst/doc/sim-Lexis.R r-cran-epi-2.7/inst/doc/sim-Lexis.R --- r-cran-epi-2.0/inst/doc/sim-Lexis.R 2016-01-04 13:51:22.000000000 +0000 +++ r-cran-epi-2.7/inst/doc/sim-Lexis.R 1970-01-01 00:00:00.000000000 +0000 @@ -1,458 +0,0 @@ -### R code from vignette source 'sim-Lexis.rnw' -### Encoding: ISO8859-1 - -################################################### -### code chunk number 1: start -################################################### -options( width=90 ) -library( Epi ) -print( sessionInfo(), l=F ) - - -################################################### -### code chunk number 2: Lexis -################################################### -data(DMlate) -dml <- Lexis( entry = list(Per=dodm, Age=dodm-dobth, DMdur=0 ), - exit = list(Per=dox), - exit.status = factor(!is.na(dodth),labels=c("DM","Dead")), - data = DMlate ) - - -################################################### -### code chunk number 3: cut -################################################### -dmi <- cutLexis( dml, cut = dml$doins, - pre = "DM", - new.state = "Ins", - new.scale = "t.Ins", - split.states = TRUE ) -summary( dmi ) -str(dmi) - - -################################################### -### code chunk number 4: boxes -################################################### -boxes( dmi, boxpos = list(x=c(20,20,80,80), - y=c(80,20,80,20)), - scale.R = 1000, show.BE = TRUE ) - - -################################################### -### code chunk number 5: split -################################################### -Si <- splitLexis( dmi, 0:30/2, "DMdur" ) -dim( Si ) -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"), - 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"), - 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 ) )) ) - - -################################################### -### code chunk number 7: Poisson -################################################### -library( splines ) -DM.Ins <- glm( (lex.Xst=="Ins") ~ Ns( Age , knots=ai.kn ) + - Ns( DMdur, knots=di.kn ) + - I(Per-2000) + sex, - family=poisson, offset=log(lex.dur), - data = subset(Si,lex.Cst=="DM") ) -DM.Dead <- glm( (lex.Xst=="Dead") ~ Ns( Age , knots=ad.kn ) + - Ns( DMdur, knots=dd.kn ) + - I(Per-2000) + sex, - family=poisson, offset=log(lex.dur), - data = subset(Si,lex.Cst=="DM") ) -Ins.Dead <- glm( (lex.Xst=="Dead(Ins)") ~ Ns( Age , knots=ad.kn ) + - Ns( DMdur, knots=dd.kn ) + - Ns( t.Ins, knots=ti.kn ) + - I(Per-2000) + sex, - family=poisson, offset=log(lex.dur), - data = subset(Si,lex.Cst=="Ins") ) - - -################################################### -### code chunk number 8: prop-haz -################################################### -with( Si, table(lex.Cst) ) -All.Dead <- glm( (lex.Xst %in% c("Dead(Ins)","Dead")) ~ - Ns( Age , knots=ad.kn ) + - Ns( DMdur, knots=dd.kn ) + - lex.Cst + - I(Per-2000) + sex, - family=poisson, offset=log(lex.dur), - data = Si ) -round( ci.exp( All.Dead ), 3 ) - - -################################################### -### code chunk number 9: get-dev -################################################### -what <- c("null.deviance","df.null","deviance","df.residual") -( rD <- unlist( DM.Dead[what] ) ) -( rI <- unlist( Ins.Dead[what] ) ) -( rA <- unlist( All.Dead[what] ) ) -round( c( dd <- rA-(rI+rD), "pVal"=1-pchisq(dd[3],dd[4]+1) ), 3 ) - - -################################################### -### code chunk number 10: pr-array -################################################### -pr.rates <- NArray( list( DMdur = seq(0,12,0.1), - DMage = 4:7*10, - r.Ins = c(NA,0,2,5), - model = c("DM/Ins","All"), - what = c("rate","lo","hi") ) ) -str( pr.rates ) - - -################################################### -### code chunk number 11: sim-Lexis.rnw:481-482 -################################################### -ci.pred - - -################################################### -### code chunk number 12: make-pred -################################################### -nd <- data.frame( DMdur = as.numeric( dimnames(pr.rates)[[1]] ), - lex.Cst = factor( 1, levels=1:4, - labels=levels(Si$lex.Cst) ), - sex = factor( 1, levels=1:2, labels=c("M","F")), - lex.dur = 1000 ) -for( ia in dimnames(pr.rates)[[2]] ) - { -dnew <- transform( nd, Age = as.numeric(ia)+DMdur, - Per = 1998+DMdur ) -pr.rates[,ia,1,"DM/Ins",] <- ci.pred( DM.Dead, newdata = dnew ) -pr.rates[,ia,1,"All" ,] <- ci.pred( All.Dead, newdata = dnew ) -for( ii in dimnames(pr.rates)[[3]][-1] ) - { -dnew = transform( dnew, lex.Cst = factor( 2, levels=1:4, - labels=levels(Si$lex.Cst) ), - t.Ins = ifelse( (DMdur-as.numeric(ii)) >= 0, - DMdur-as.numeric(ii), NA ) ) -pr.rates[,ia, ii ,"DM/Ins",] <- ci.pred( Ins.Dead, newdata = dnew ) -pr.rates[,ia, ii ,"All" ,] <- ci.pred( All.Dead, newdata = dnew ) - } - } - - -################################################### -### code chunk number 13: mort-int -################################################### -par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, las=1 ) -plot( NA, xlim=c(40,82), ylim=c(5,300), bty="n", - log="y", xlab="Age", ylab="Mortality rate per 1000 PY" ) -abline( v=seq(40,80,5), h=outer(1:9,10^(0:2),"*"), col=gray(0.8) ) -for( aa in 4:7*10 ) for( ii in 1:4 ) - matlines( aa+as.numeric(dimnames(pr.rates)[[1]]), - cbind( pr.rates[,paste(aa),ii,"DM/Ins",], - pr.rates[,paste(aa),ii,"All" ,] ), - type="l", lty=1, lwd=c(3,1,1), - col=rep(c("red","limegreen"),each=3) ) - - -################################################### -### code chunk number 14: Tr -################################################### -Tr <- list( "DM" = list( "Ins" = DM.Ins, - "Dead" = DM.Dead ), - "Ins" = list( "Dead(Ins)" = Ins.Dead ) ) - - -################################################### -### code chunk number 15: make-ini -################################################### -str( Si[NULL,1:9] ) -ini <- subset(Si,FALSE,select=1:9) -str( ini ) -ini <- subset(Si,select=1:9)[NULL,] -str( ini ) - - -################################################### -### code chunk number 16: ini-fill -################################################### -ini[1:2,"lex.id"] <- 1:2 -ini[1:2,"lex.Cst"] <- "DM" -ini[1:2,"Per"] <- 1995 -ini[1:2,"Age"] <- 60 -ini[1:2,"DMdur"] <- 5 -ini[1:2,"sex"] <- c("M","F") -ini - - -################################################### -### code chunk number 17: simL -################################################### -Nsim <- 5000 -system.time( simL <- simLexis( Tr, - ini, - t.range = 12, - N = Nsim ) ) - - -################################################### -### code chunk number 18: sum-simL -################################################### -summary( simL, by="sex" ) - - -################################################### -### code chunk number 19: Tr.p-simP -################################################### -Tr.p <- list( "DM" = list( "Ins" = DM.Ins, - "Dead" = All.Dead ), - "Ins" = list( "Dead(Ins)" = All.Dead ) ) -system.time( simP <- simLexis( Tr.p, - ini, - t.range = 12, - N = Nsim ) ) -summary( simP, by="sex" ) - - -################################################### -### code chunk number 20: Cox-dur -################################################### -library( survival ) -Cox.Dead <- coxph( Surv( DMdur, DMdur+lex.dur, - lex.Xst %in% c("Dead(Ins)","Dead")) ~ - Ns( Age-DMdur, knots=ad.kn ) + - I(lex.Cst=="Ins") + - I(Per-2000) + sex, - data = Si ) -round( ci.exp( Cox.Dead ), 3 ) -round( ci.exp( All.Dead ), 3 ) - - -################################################### -### code chunk number 21: TR.c -################################################### -Tr.c <- list( "DM" = list( "Ins" = Tr$DM$Ins, - "Dead" = Cox.Dead ), - "Ins" = list( "Dead(Ins)" = Cox.Dead ) ) -system.time( simC <- simLexis( Tr.c, - ini, - t.range = 12, - N = Nsim ) ) -summary( simC, by="sex" ) - - -################################################### -### code chunk number 22: nState -################################################### -system.time( -nSt <- nState( subset(simL,sex=="M"), - at=seq(0,11,0.2), from=1995, time.scale="Per" ) ) -nSt[1:10,] - - -################################################### -### code chunk number 23: pstate0 -################################################### -pM <- pState( nSt, perm=c(1,2,4,3) ) -head( pM ) -par( mfrow=c(1,2), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 ) -plot( pM ) -plot( pM, border="black", col="transparent", lwd=3 ) -text( rep(as.numeric(rownames(pM)[nrow(pM)-1]),ncol(pM)), - pM[nrow(pM),]-diff(c(0,pM[nrow(pM),]))/5, - colnames( pM ), adj=1 ) -box( col="white", lwd=3 ) -box() - - -################################################### -### code chunk number 24: pstatex -################################################### -clr <- c("limegreen","orange") -# expand with a lighter version of the two chosen colors -clx <- c( clr, rgb( t( col2rgb( clr[2:1] )*2 + rep(255,3) ) / 3, max=255 ) ) -par( mfrow=c(1,2), las=1, mar=c(3,3,4,2), mgp=c(3,1,0)/1.6 ) -# Men -plot( pM, col=clx ) -lines( as.numeric(rownames(pM)), pM[,2], lwd=3 ) -mtext( "60 year old male, diagnosed 1990, aged 55", side=3, line=2.5, adj=0, col=gray(0.6) ) -mtext( "Survival curve", side=3, line=1.5, adj=0 ) -mtext( "DM, no insulin DM, Insulin", side=3, line=0.5, adj=0, col=clr[1] ) -mtext( "DM, no insulin", side=3, line=0.5, adj=0, col=clr[2] ) -axis( side=4 ) -axis( side=4, at=1:19/20, labels=FALSE ) -axis( side=4, at=1:99/100, labels=FALSE, tcl=-0.3 ) -# Women -pF <- pState( nState( subset(simL,sex=="F"), - at=seq(0,11,0.2), - from=1995, - time.scale="Per" ), - perm=c(1,2,4,3) ) -plot( pF, col=clx ) -lines( as.numeric(rownames(pF)), pF[,2], lwd=3 ) -mtext( "60 year old female, diagnosed 1990, aged 55", side=3, line=2.5, adj=0, col=gray(0.6) ) -mtext( "Survival curve", side=3, line=1.5, adj=0 ) -mtext( "DM, no insulin DM, Insulin", side=3, line=0.5, adj=0, col=clr[1] ) -mtext( "DM, no insulin", side=3, line=0.5, adj=0, col=clr[2] ) -axis( side=4 ) -axis( side=4, at=1:19/20, labels=FALSE ) -axis( side=4, at=1:99/100, labels=FALSE, tcl=-0.3 ) - - -################################################### -### code chunk number 25: pstatey -################################################### -par( mfrow=c(1,2), las=1, mar=c(3,3,4,2), mgp=c(3,1,0)/1.6 ) -# Men -pM <- pState( nState( subset(simL,sex=="M"), - at=seq(0,11,0.2), - from=60, - time.scale="Age" ), - perm=c(1,2,4,3) ) -plot( pM, col=clx, xlab="Age" ) -lines( as.numeric(rownames(pM)), pM[,2], lwd=3 ) -mtext( "60 year old male, diagnosed 1990, aged 55", side=3, line=2.5, adj=0, col=gray(0.6) ) -mtext( "Survival curve", side=3, line=1.5, adj=0 ) -mtext( "DM, no insulin DM, Insulin", side=3, line=0.5, adj=0, col=clr[1] ) -mtext( "DM, no insulin", side=3, line=0.5, adj=0, col=clr[2] ) -axis( side=4 ) -axis( side=4, at=1:19/20, labels=FALSE ) -axis( side=4, at=1:19/20, labels=FALSE, tcl=-0.4 ) -axis( side=4, at=1:99/100, labels=FALSE, tcl=-0.3 ) -# Women -pF <- pState( nState( subset(simL,sex=="F"), - at=seq(0,11,0.2), - from=60, - time.scale="Age" ), - perm=c(1,2,4,3) ) -plot( pF, col=clx, xlab="Age" ) -lines( as.numeric(rownames(pF)), pF[,2], lwd=3 ) -mtext( "60 year old female, diagnosed 1990, aged 55", side=3, line=2.5, adj=0, col=gray(0.6) ) -mtext( "Survival curve", side=3, line=1.5, adj=0 ) -mtext( "DM, no insulin DM, Insulin", side=3, line=0.5, adj=0, col=clr[1] ) -mtext( "DM, no insulin", side=3, line=0.5, adj=0, col=clr[2] ) -axis( side=4 ) -axis( side=4, at=1:9/10, labels=FALSE ) -axis( side=4, at=1:19/20, labels=FALSE, tcl=-0.4 ) -axis( side=4, at=1:99/100, labels=FALSE, tcl=-0.3 ) - - -################################################### -### code chunk number 26: comp-0 -################################################### -PrM <- pState( nState( subset(simP,sex=="M"), - at=seq(0,11,0.2), - from=60, - time.scale="Age" ), - perm=c(1,2,4,3) ) -PrF <- pState( nState( subset(simP,sex=="F"), - at=seq(0,11,0.2), - from=60, - time.scale="Age" ), - perm=c(1,2,4,3) ) -CoxM <- pState( nState( subset(simC,sex=="M"), - at=seq(0,11,0.2), - from=60, - time.scale="Age" ), - perm=c(1,2,4,3) ) -CoxF <- pState( nState( subset(simC,sex=="F"), - at=seq(0,11,0.2), - from=60, - time.scale="Age" ), - perm=c(1,2,4,3) ) - -par( mfrow=c(1,2), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 ) - plot( pM, border="black", col="transparent", lwd=3 ) -lines( PrM, border="blue" , col="transparent", lwd=3 ) -lines( CoxM, border="red" , col="transparent", lwd=3 ) -text( 60.5, 0.05, "M" ) -box( lwd=3 ) - - plot( pF, border="black", col="transparent", lwd=3 ) -lines( PrF, border="blue" , col="transparent", lwd=3 ) -lines( CoxF, border="red" , col="transparent", lwd=3 ) -text( 60.5, 0.05, "F" ) -box( lwd=3 ) - - -################################################### -### code chunk number 27: sim-Lexis.rnw:987-988 -################################################### -options( keep.source=TRUE ) - - -################################################### -### code chunk number 28: sim-Lexis.rnw:1004-1007 -################################################### -cbind( -attr( ini, "time.scale" ), -attr( ini, "time.since" ) ) - - -################################################### -### code chunk number 29: sim-Lexis.rnw:1032-1033 -################################################### -simLexis - - -################################################### -### code chunk number 30: sim-Lexis.rnw:1050-1051 -################################################### -Epi:::simX - - -################################################### -### code chunk number 31: sim-Lexis.rnw:1063-1064 -################################################### -Epi:::sim1 - - -################################################### -### code chunk number 32: sim-Lexis.rnw:1076-1077 -################################################### -Epi:::lint - - -################################################### -### code chunk number 33: sim-Lexis.rnw:1087-1088 -################################################### -Epi:::get.next - - -################################################### -### code chunk number 34: sim-Lexis.rnw:1097-1098 -################################################### -Epi:::chop.lex - - -################################################### -### code chunk number 35: sim-Lexis.rnw:1115-1116 -################################################### -nState - - -################################################### -### code chunk number 36: sim-Lexis.rnw:1125-1126 -################################################### -pState - - -################################################### -### code chunk number 37: sim-Lexis.rnw:1130-1132 -################################################### -plot.pState -lines.pState - - diff -Nru r-cran-epi-2.0/inst/doc/simLexis.R r-cran-epi-2.7/inst/doc/simLexis.R --- r-cran-epi-2.0/inst/doc/simLexis.R 2016-01-06 18:15:40.000000000 +0000 +++ r-cran-epi-2.7/inst/doc/simLexis.R 2016-10-03 19:59:16.000000000 +0000 @@ -121,7 +121,7 @@ ################################################### -### code chunk number 11: simLexis.rnw:510-511 +### code chunk number 11: simLexis.rnw:509-510 ################################################### ci.pred @@ -389,13 +389,13 @@ ################################################### -### code chunk number 27: simLexis.rnw:1028-1029 +### code chunk number 27: simLexis.rnw:1027-1028 ################################################### options( keep.source=TRUE ) ################################################### -### code chunk number 28: simLexis.rnw:1045-1048 +### code chunk number 28: simLexis.rnw:1044-1047 ################################################### cbind( attr( ini, "time.scale" ), @@ -403,55 +403,55 @@ ################################################### -### code chunk number 29: simLexis.rnw:1073-1074 +### code chunk number 29: simLexis.rnw:1072-1073 ################################################### simLexis ################################################### -### code chunk number 30: simLexis.rnw:1091-1092 +### code chunk number 30: simLexis.rnw:1090-1091 ################################################### Epi:::simX ################################################### -### code chunk number 31: simLexis.rnw:1104-1105 +### code chunk number 31: simLexis.rnw:1103-1104 ################################################### Epi:::sim1 ################################################### -### code chunk number 32: simLexis.rnw:1117-1118 +### code chunk number 32: simLexis.rnw:1116-1117 ################################################### Epi:::lint ################################################### -### code chunk number 33: simLexis.rnw:1128-1129 +### code chunk number 33: simLexis.rnw:1127-1128 ################################################### Epi:::get.next ################################################### -### code chunk number 34: simLexis.rnw:1138-1139 +### code chunk number 34: simLexis.rnw:1137-1138 ################################################### Epi:::chop.lex ################################################### -### code chunk number 35: simLexis.rnw:1156-1157 +### code chunk number 35: simLexis.rnw:1155-1156 ################################################### nState ################################################### -### code chunk number 36: simLexis.rnw:1166-1167 +### code chunk number 36: simLexis.rnw:1165-1166 ################################################### pState ################################################### -### code chunk number 37: simLexis.rnw:1171-1173 +### code chunk number 37: simLexis.rnw:1170-1172 ################################################### plot.pState lines.pState diff -Nru r-cran-epi-2.0/inst/doc/simLexis.rnw r-cran-epi-2.7/inst/doc/simLexis.rnw --- r-cran-epi-2.0/inst/doc/simLexis.rnw 2016-01-06 18:15:40.000000000 +0000 +++ r-cran-epi-2.7/inst/doc/simLexis.rnw 2016-10-03 19:59:16.000000000 +0000 @@ -31,7 +31,6 @@ \usepackage[ae,hyper]{Rd} \usepackage[dvipsnames]{xcolor} \usepackage[super]{nth} -\usepackage[retainorgcmds]{IEEEtrantools} \usepackage{makeidx,Sweave,floatflt,amsmath,amsfonts,amsbsy,enumitem,dcolumn,needspace} \usepackage{ifthen,calc,eso-pic,everyshi} \usepackage{booktabs,longtable,rotating,graphicx} diff -Nru r-cran-epi-2.0/man/apc.fit.Rd r-cran-epi-2.7/man/apc.fit.Rd --- r-cran-epi-2.0/man/apc.fit.Rd 2015-05-27 15:27:03.000000000 +0000 +++ r-cran-epi-2.7/man/apc.fit.Rd 2016-05-25 06:50:31.000000000 +0000 @@ -6,7 +6,7 @@ \description{ Fits the classical five models to tabulated rate data (cases, person-years) classified by two of age, period, cohort: - Age, Age-drift, Age-Period, Age-Cohort and Age-period. There are no + Age, Age-drift, Age-Period, Age-Cohort and Age-Period-Cohort. There are no assumptions about the age, period or cohort classes being of the same length, or that tabulation should be only by two of the variables. Only requires that mean age and period for each tabulation unit is given. @@ -21,7 +21,7 @@ ref.p, dist = c("poisson","binomial"), model = c("ns","bs","ls","factor"), - dr.extr = c("weighted","Holford"), + dr.extr = "weighted", parm = c("ACP","APC","AdCP","AdPC","Ad-P-C","Ad-C-P","AC-P","AP-C"), npar = c( A=5, P=5, C=5 ), scale = 1, @@ -32,7 +32,7 @@ \item{data}{Data frame with (at least) variables, \code{A} (age), \code{P} (period), \code{D} (cases, deaths) and \code{Y} (person-years). Cohort (date of birth) is computed as \code{P-A}. - If thsi argument is given the arguments \code{A}, \code{P}, + If this argument is given the arguments \code{A}, \code{P}, \code{D} and \code{Y} are ignored.} \item{A}{Age; numerical vector with mean age at diagnosis for each unit.} \item{P}{Period; numerical vector with mean date of diagnosis for each @@ -42,11 +42,11 @@ data, see the \code{dist} argument.} \item{ref.c}{Reference cohort, numerical. Defaults to median date of birth among cases. If used with \code{parm="AdCP"} or \code{parm="AdPC"}, - the resdiual cohort effects will be 1 at \code{ref.c}} + the residual cohort effects will be 1 at \code{ref.c}} \item{ref.p}{Reference period, numerical. Defaults to median date of diagnosis among cases.} - \item{dist}{Distribution (or more precisely: Likelihood) used for modelling. - if a binomial model us ised, \code{Y} is assuemd to be the + \item{dist}{Distribution (or more precisely: Likelihood) used for modeling. + if a binomial model us used, \code{Y} is assumed to be the denominator; \code{"binomial"} gives a binomial model with logit link.} \item{model}{Type of model fitted: @@ -61,17 +61,33 @@ is ignored in this case. } } - \item{dr.extr}{Character. How the drift parameter should be extracted from - the age-period-cohort model. \code{"weighted"} (default) lets the - weighted average (by marginal no. cases, \code{D}) of the estimated - period and cohort effects have 0 slope. \code{"Holford"} uses the - naive average over all values for the estimated effects, - disregarding the no. cases.} + \item{dr.extr}{Character or numeric. + How the drift parameter should be extracted from the + age-period-cohort model. Specifies the inner product used for + definition of orthogonality of the period / cohort effects to the + linear effects --- in terms of a diagonal matrix. + + \code{"weighted"} (or "t") (default) uses the no. cases, \code{D}, + corresponding to the observed information about the log-rate + (usually termed "theta", hence the "\code{t}"). \code{"r"} or \code{"l"} uses + \code{Y*Y/D} corresponding to the observed information about the + rate (usually termed "lambda", hence the "\code{l}"). \code{"y"} + uses the person-years as the weight in the inner product. If given + "\code{n}" (Naive) (well, in fact any other character value) will + induce the use of the standard inner product putting equal weight on + all units in the dataset. + + If given as a numeric vector this is used as the diagonal of the + matrix inducing the inner product. + + The setting of this parameter has no effect on the fit of the model, + only on the parametrization. + } \item{parm}{Character. Indicates the parametrization of the effects. - The first four refer to the ML-fit of the Age-Period-Cohort model, the last four - give Age-effects from a smaller model and residuals relative to - this. If one of the latter is chosen, the argument \code{dr.extr} - is ignored. Possible values for \code{parm} are: + The first four refer to the ML-fit of the Age-Period-Cohort model, + the last four give Age-effects from a smaller model and residuals + relative to this. If one of the latter is chosen, the argument + \code{dr.extr} is ignored. Possible values for \code{parm} are: \itemize{ \item \code{"ACP"}: ML-estimates. Age-effects as rates for the reference cohort. Cohort effects as RR relative to the reference @@ -81,14 +97,14 @@ period. Cohort effects constrained to be 0 on average with 0 slope. \item \code{"AdCP"}: ML-estimates. Age-effects as rates for the reference cohort. Cohort and period effects constrained to be 0 on - average with 0 slope. These effects do not multiply to the fitted - rates, the drift is missing and needs to be included to produce - the fitted values. + average with 0 slope. In this case returned effects do not + multiply to the fitted rates, the drift is missing and needs to be + included to produce the fitted values. \item \code{"AdPC"}: ML-estimates. Age-effects as rates for the reference period. Cohort and period effects constrained to be 0 on - average with 0 slope. These effects do not multiply to the fitted - rates, the drift is missing and needs to be included to produce - the fitted values. + average with 0 slope. In this case returned effects do not + multiply to the fitted rates, the drift is missing and needs to be + included to produce the fitted values. \item \code{"Ad-C-P"}: Age effects are rates for the reference cohort in the Age-drift model (cohort drift). Cohort effects are from the model with cohort alone, using log(fitted values) from the Age-drift @@ -112,14 +128,14 @@ } } \item{npar}{The number of parameters/knots to use for each of the terms in the model. If it is vector of length 3, the numbers are taken as the - no. of knots for Age, Period and Cohort, respctively. Unless it has + no. of knots for Age, Period and Cohort, respectively. Unless it has a names attribute with vales "A", "P" and "C" in which case these will be used. The knots chosen are the quantiles - \code{(1:nk+0.1)/(nk+0.2)} of the events (i.e. of \code{rep(A,D)}) + \code{(1:nk-0.5)/nk} of the events (i.e. of \code{rep(A,D)}) \code{npar} may also be a named list of three numerical vectors with names "A", "P" and "C", in which case these taken as the knots for - the age, period and cohort effect, the first and last element in + the age, period and cohort effect, the smallest and largest element in each vector are used as the boundary knots.} \item{alpha}{The significance level. Estimates are given with (1-\code{alpha}) confidence limits.} @@ -132,19 +148,21 @@ \code{\link{apc.lines}}) --- a list with components: \item{Type}{Text describing the model and parametrization returned} \item{Model}{The model object(s) on which the parametrization is based.} - \item{Age}{Matrix with 4 colums: \code{A.pt} with the ages (equals + \item{Age}{Matrix with 4 columns: \code{A.pt} with the ages (equals \code{unique(A)}) and three columns giving the estimated rates with c.i.s.} - \item{Per}{Matrix with 4 colums: \code{P.pt} with the dates of + \item{Per}{Matrix with 4 columns: \code{P.pt} with the dates of diagnosis (equals \code{unique(P)}) and three columns giving the estimated RRs with c.i.s.} - \item{Coh}{Matrix with 4 colums: \code{C.pt} with the dates of birth + \item{Coh}{Matrix with 4 columns: \code{C.pt} with the dates of birth (equals \code{unique(P-A)}) and three columns giving the estimated RRs with c.i.s.} - \item{Drift}{A 3 column matrix with drift-estimates and c.i.s: The first row is - the ML-estimate of the drift (as defined by \code{drift}), the - second row is the estimate from the Age-drift model. For the - sequential parametrizations, only the latter is given.} + \item{Drift}{A 3 column matrix with drift-estimates and c.i.s: The + first row is the ML-estimate of the drift (as defined by + \code{drift}), the second row is the estimate from the Age-drift + model. The first row name indicates which type of inner product were + used for projections. For the sequential parametrizations, only the + latter is given.} \item{Ref}{Numerical vector of length 2 with reference period and cohort. If ref.p or ref.c was not supplied the corresponding element is NA.} \item{AOV}{Analysis of deviance table comparing the five classical @@ -153,21 +171,34 @@ \item{Knots}{If \code{model} is one of \code{"ns"} or \code{"bs"}, a list with three components: \code{Age}, \code{Per}, \code{Coh}, each one a vector of knots. The max and the min are the boundary knots.} +} +\details{ + Each record in the input data frame represents a subset of a Lexis + diagram. The subsets need not be of equal length on the age and + period axes, in fact there are no restrictions on the shape of + these; they could be Lexis triangels for example. The requirement is + that \code{A} and \code{P} are coded with the mean age and calendar + time of observation in the subset. This is essential since \code{A} + and \code{P} are used as quantitative variables in the models. + + This is a different approach relative to the vast majority of the + uses of APC-models in the literature where a factor model is used for + age, perido and cohort effects. The latter can be obtained by using + \code{model="factor"}. } \references{ The considerations behind the parametrizations used in this function - are given in details in a preprint from Department of Biostatistics - in Copenhagen: - "Demography and epidemiology: Age-Period-Cohort models in the computer age", - \url{http://biostat.ku.dk/reports/2006/ResearchReport06-1.pdf/}, later - published as: - B. Carstensen: Age-period-cohort models for the Lexis diagram. + are given in details in: + B. Carstensen: Age-Period-Cohort models for the Lexis diagram. Statistics in Medicine, 10; 26(15):3018-45, 2007. + + Various links to course material etc. is available through \url{http://BendixCarstensen.com/APC} } \author{ Bendix Carstensen, \url{http://BendixCarstensen.com} } \seealso{ + \code{\link{LCa.fit}}, \code{\link{apc.frame}}, \code{\link{apc.lines}}, \code{\link{apc.plot}}. diff -Nru r-cran-epi-2.0/man/apc.LCa.Rd r-cran-epi-2.7/man/apc.LCa.Rd --- r-cran-epi-2.0/man/apc.LCa.Rd 1970-01-01 00:00:00.000000000 +0000 +++ r-cran-epi-2.7/man/apc.LCa.Rd 2016-05-25 07:43:03.000000000 +0000 @@ -0,0 +1,89 @@ +\name{apc.LCa} +\alias{apc.LCa} +\alias{show.apc.LCa} +\title{Fit Age-Period-Cohort models and Lee-Carter models with effects + modeled by natural splines. +} +\description{ +\code{apc.LCa} fits an Age-Period-Cohort model and sub-models (using +\code{\link{apc.fit}}) as well as Lee-Carter models (using +\code{\link{LCa.fit}}). \code{boxes.apc.LCa} plots the models in little +boxes with their residual deviance with arrows showing their +relationships. +} +\usage{ +apc.LCa( data, + keep.models = FALSE, + ... ) +show.apc.LCa( x, + dev.scale = TRUE, + top = "Ad", ... ) +} +\arguments{ + \item{data}{A data frame that must have columns \code{A}, \code{P}, + \code{D} and \code{Y}, see e.g. \code{\link{apc.fit}} + } + \item{keep.models}{Logical. Should the \code{apc} object and the 5 + \code{LCa} objects be returned too? + } + \item{...}{Further parameters passed on to \code{\link{LCa.fit}} or \code{\link{boxes.matrix}}. + } + \item{x}{The result from a call to \code{apc.LCa}.} + \item{dev.scale}{Should the vertical position of the boxes with the + models be scales relative to the deviance between the Age-drift + model and the extended Lee-Carter model?} + \item{top}{The model presented at the top of the plot of boxes + (together with any other model with larger deviance) when + vertical position is scaled by deviances. Only "Ad", "AP", "AC", + "APa" or "ACa" will make sense.} +} +\details{The function \code{apc.LCa} fits all 9 models (well, 10) available as + extension and sub-models of the APC-model and compares them by + returning deviance and residual df. +} +\value{A 9 by 2 matrix classified by model and deviance/df; optionally + (if \code{models=TRUE}) a list with the matrix as \code{dev}, \code{apc}, an + \code{apc} object (from \code{\link{apc.fit}}), and \code{LCa}, a list + with 5 \code{LCa} objects (from \code{\link{LCa.fit}}). +} +\author{ +Bendix Carstensen, \url{http://BendixCarstensen.com} +} +\seealso{ +\code{ \link{apc.fit}}, \code{\link{LCa.fit} } +} +\examples{ +library( Epi ) + +# Danish lung cancer incidence in 5x5x5 Lexis triangles +data( lungDK ) +lc <- subset( lungDK, Ax>40 )[,c("Ax","Px","D","Y")] +names( lc )[1:2] <- c("A","P") +head( lc ) + +al <- apc.LCa( lc, npar=c(9,6,6,6,10), keep.models=TRUE, maxit=500, eps=10e-3 ) +show.apc.LCa( al, dev=FALSE ) +show.apc.LCa( al, top="AP" ) +show.apc.LCa( al, top="APa" ) +show.apc.LCa( al, top="ACa" ) + +# Danish mortality data +\dontrun{ +data( M.dk ) +mdk <- subset( M.dk, sex==1 )[,c("A","P","D","Y")] +head( mdk ) + +al <- apc.LCa( mdk, npar=c(15,15,20,6,6), maxit=50, eps=10e-3, + quiet=FALSE, VC=FALSE ) +show.apc.LCa( al, dev=FALSE ) +show.apc.LCa( al, dev=TRUE ) +show.apc.LCa( al, top="AP" ) + +# Fit a reasonable model to Danish mortality data and plot results +mACa <- LCa.fit( mdk, model="ACa", npar=c(15,15,20,6,6), c.ref=1930, + a.ref=70, quiet=FALSE, maxit=250 ) +par( mfrow=c(1,3) ) +plot( mACa )} +} +\keyword{regression} +\keyword{models} diff -Nru r-cran-epi-2.0/man/boxes.MS.Rd r-cran-epi-2.7/man/boxes.MS.Rd --- r-cran-epi-2.0/man/boxes.MS.Rd 2015-05-27 15:27:03.000000000 +0000 +++ r-cran-epi-2.7/man/boxes.MS.Rd 2016-02-15 16:05:15.000000000 +0000 @@ -121,7 +121,8 @@ \item{digits.Y}{How many digits after the decimal point should be used for the person-years.} \item{show.BE}{Logical. Should number of persons beginning - resp. ending follow up in each state be shown?} + resp. ending follow up in each state be shown? If given as charcater + "nz" or "noz" the numbers will be shown, but zeros omitted.} \item{BE.sep}{Character vector of length 4, used for annotation of the number of persons beginning and ending in each state: 1st elemet precedes no. beginning, 2nd trails it, 3rd precedes the no. ending diff -Nru r-cran-epi-2.0/man/ci.lin.Rd r-cran-epi-2.7/man/ci.lin.Rd --- r-cran-epi-2.0/man/ci.lin.Rd 2015-08-03 07:04:22.000000000 +0000 +++ r-cran-epi-2.7/man/ci.lin.Rd 2016-10-01 09:59:29.000000000 +0000 @@ -140,7 +140,6 @@ Bendix Carstensen, \url{BendixCarstensen.com} & Michael Hills - \url{http://www.mhills.pwp.blueyonder.co.uk/} } \seealso{See also \code{\link{ci.cum}} for a function computing cumulative sums of (functions of) parameter estimates.} diff -Nru r-cran-epi-2.0/man/erl.Rd r-cran-epi-2.7/man/erl.Rd --- r-cran-epi-2.0/man/erl.Rd 1970-01-01 00:00:00.000000000 +0000 +++ r-cran-epi-2.7/man/erl.Rd 2016-10-01 11:25:33.000000000 +0000 @@ -0,0 +1,203 @@ +\name{erl} +\alias{surv1} +\alias{surv2} +\alias{erl1} +\alias{erl} +\alias{yll} +\title{Compute survival functions from rates and expected residual + lifetime in an illness-death model as well as years of life lost to disease. +} +\description{ + These functions compute survival functions from a set of mortality and + disease incidence rates in an illness-death model. Expected residual + life time can be computed under various scenarios by the \code{erl} + function, and areas between survival functions can be computed under + various scenarios by the \code{yll} function. Rates are assumed + supplied for equidistant intervals of length \code{int}. + } +\usage{ + surv1( int, mu , age.in = 0, A = NULL ) + erl1( int, mu , age.in = 0 ) + surv2( int, muW, muD, lam, age.in = 0, A = NULL ) + erl( int, muW, muD, lam=NULL, age.in = 0, A = NULL, + immune = is.null(lam), yll=TRUE, note=TRUE ) + yll( int, muW, muD, lam=NULL, age.in = 0, A = NULL, + immune = is.null(lam), note=TRUE ) + } +\arguments{ + \item{int}{ + Scalar. Length of intervals that rates refer to. + } + \item{mu}{ + Numeric vector of mortality rates at midpoints of intervals of length \code{int} + } + \item{muW}{ + Numeric vector of mortality rates among persons in the "Well" state at + midpoints of intervals of length \code{int}. Left endpoint of first + interval is \code{age.in}. + } + \item{muD}{ + Numeric vector of mortality rates among persons in the "Diseased" state + at midpoints of intervals of length \code{int}. Left endpoint of first + interval is \code{age.in}. + } + \item{lam}{ + Numeric vector of disease incidence rates among persons in the "Well" state + at midpoints of intervals of length \code{int}. Left endpoint of first + interval is \code{age.in}. + } + \item{age.in}{ + Scalar indicating the age at the left endpoint of the first interval. + } + \item{A}{ + Numeric vector of conditioning ages for calculation of survival + functions. + } + \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}. + } + \item{note}{ + Logical. Should a warning of silly assumptions be printed? + } + \item{yll}{ + Logical. Should years of life lost be included in the result? + } +} +\details{ + The mortality rates given are supposed to refer to the ages + \code{age.in+i*int/2}, \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 + given. Thus if \code{int} is given in years, the rates must be given + in the unit of events per year. + + The ages in which the survival curves are computed are from + \code{age.in} and then at the end of \code{length(muW)} + (\code{length(mu)}) intervals each of length \code{int}. + + The \code{age.in} argument is merely a device to account for rates + only available from a given age. It has two effects, one is that + labeling of the interval endpoint is offset by this quantity, thus + 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.} + +\value{\code{surv1} and \code{surv2} return a matrix whose first column + is the ages at the ends of the + intervals, thus with \code{length(mu)+1} rows. The following columns + are the survival functions (since \code{age.in}), and conditional on + survival till ages as indicated in \code{A}, thus a matrix with + \code{length(A)+2} columns. Columns are labeled with the actual + conditioning ages; if \code{A} contains values that are not among the + endpoints of the intervals used, the nearest smaller interval border + is used as conditioning age, and columns are named accordingly. + + \code{surv1} returns the survival function for a simple model with one + type of death, occurring at intensity \code{mu}. + + \code{surv2} returns the survival function for a person in the "Well" + state of an illness-death model, taking into account that the person + may move to the "Diseased" state, thus requiring all three transition + rates to be specified. The conditional survival functions are + conditional on being in the "Well" state at ages given in \code{A}. + + \code{erl1} returns a three column matrix with columns \code{age}, + \code{surv} (survival function) and \code{erl} (expected residual life + time) with \code{length(mu)+1} rows. + + \code{erl} returns a two column matrix, columns labeled "Well" and + "Dis", and with row-labels \code{A}. The entries are the expected + residual life times given survival to \code{A}. If \code{yll=TRUE} the + difference between the columns is added as a + third column, labeled "YLL". + } +\author{Bendix Carstensen, \email{bxc@steno.dk} + } +\seealso{ +\code{\link{ci.cum}} +} +\examples{ +library( Epi ) +data( DMlate ) +# Naive Lexis object +Lx <- Lexis( entry = list( age = dodm-dobth ), + exit = list( age = dox -dobth ), + exit.status = factor( !is.na(dodth), labels=c("DM","Dead") ), + data = DMlate ) +# Cut follow-up at insulin inception +Lc <- cutLexis( Lx, cut = Lx$doins-Lx$dob, + new.state = "DM/ins", + precursor.states = "DM" ) +summary( Lc ) +# Split in small age intervals +Sc <- splitLexis( Lc, breaks=seq(0,120,1) ) +summary( Sc ) + +# Overview of object +boxes( Sc, boxpos=TRUE, show.BE=TRUE, scale.R=100 ) + +# Knots for splines +a.kn <- 2:9*10 + +# Mortality among DM +mW <- glm( lex.Xst=="Dead" ~ Ns( age, knots=a.kn ), + offset = log(lex.dur), + family = poisson, + data = subset(Sc,lex.Cst=="DM") ) + +# Mortality among insulin treated +mI <- update( mW, data = subset(Sc,lex.Cst=="DM/ins") ) + +# Total motality +mT <- update( mW, data = Sc ) + +# Incidence of insulin inception +lI <- update( mW, lex.Xst=="DM/ins" ~ . ) + +# From these we can now derive the fitted rates in intervals of 1 year's +# length. In real applications you would use much smaller interval like +# 1 month: +# int <- 1/12 +int <- 1 + +# Prediction frame to return rates in units of cases per 1 year +# - we start at age 40 since rates of insulin inception are largely +# indeterminate before age 40 +nd <- data.frame( age = seq( 40+int, 110, int ) - int/2, + lex.dur = 1 ) +muW <- predict( mW, newdata = nd, type = "response" ) +muD <- predict( mI, newdata = nd, type = "response" ) +lam <- predict( lI, newdata = nd, type = "response" ) + +# Compute the survival function, and the conditional from ages 50 resp. 70 +s1 <- surv1( int, muD, age.in=40, A=c(50,70) ) +round( s1, 3 ) + +s2 <- surv2( int, muW, muD, lam, age.in=40, A=c(50,70) ) +round( s2, 3 ) + +# How much is YLL overrated by ignoring insulin incidence? +round( YLL <- cbind( +yll( int, muW, muD, lam, A = 41:90, age.in = 40 ), +yll( int, muW, muD, lam, A = 41:90, age.in = 40, immune=TRUE ) ), 2 )[seq(1,51,10),] + +par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, bty="n", las=1 ) +matplot( 40:90, YLL, + type="l", lty=1, lwd=3, + ylim=c(0,10), yaxs="i", xlab="Age" ) +} +\keyword{survival} diff -Nru r-cran-epi-2.0/man/foreign.Lexis.Rd r-cran-epi-2.7/man/foreign.Lexis.Rd --- r-cran-epi-2.0/man/foreign.Lexis.Rd 2015-05-27 15:27:03.000000000 +0000 +++ r-cran-epi-2.7/man/foreign.Lexis.Rd 2016-05-25 07:08:03.000000000 +0000 @@ -86,6 +86,6 @@ } \seealso{ \code{\link{stack.Lexis}}, - \code{\link[mstate:msprep]{msdata}}, + \code{\link[mstate:msprep]{msprep}}, \code{\link[etm:etm]{etm}} } \keyword{survival} diff -Nru r-cran-epi-2.0/man/gen.exp.Rd r-cran-epi-2.7/man/gen.exp.Rd --- r-cran-epi-2.0/man/gen.exp.Rd 2015-05-27 15:27:03.000000000 +0000 +++ r-cran-epi-2.7/man/gen.exp.Rd 2016-08-12 05:28:11.000000000 +0000 @@ -21,25 +21,26 @@ } \arguments{ \item{purchase}{Data frame with columns \code{id}-person id, - \code{dop}-date of purchase, \code{amt}-amount purchased, and - optionally \code{dpt}-defined daily dose, that is how much is assumed - to be ingested per unit time. The time unit used here is assumed to be -the same as that used in \code{dop}, so despite the name it is not -necessarily measured per day.} - \item{id}{Name of the id variable in the data frame.} - \item{dop}{Name of the date of purchase variable in the data frame.} - \item{amt}{Name of the amount purchased variable in the data frame.} - \item{dpt}{Name of the dose-per-time variable in the data frame.} + \code{dop} - \code{d}ate \code{o}f \code{p}urchase, \code{amt} - + \code{am}oun\code{t} purchased, and optionally \code{dpt} - + (\code{d}ose \code{p}er \code{t}ime) ("defined daily dose", DDD, + that is), how much is assume to be ingested per unit time. The + 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 id variable must have the same name as in the \code{purchase} data frame.} - \item{doe}{Name of the date of entry variable.} - \item{dox}{Name of the date of exit variable.} - \item{use.dpt}{Logical, should we use information on dose per time.} - \item{breaks}{Numerical vector of time points where the time since exposure and the - cumulative dose are computed.} + \item{doe}{Character. Name of the date of entry variable.} + \item{dox}{Character. Name of the date of exit 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.} \item{lags}{Numerical vector of lag-times used in computing lagged cumulative doses.} - \item{push.max}{How much can purchases maximally be pushed forward in + \item{push.max}{Numerical. How much can purchases maximally be pushed forward in time. See details.} \item{pred.win}{The length of the window used for constructing the average dose per time used to compute the duration of the last purchase} @@ -52,29 +53,40 @@ If \code{use.dpt} is \code{TRUE} then the dose per time information is used to compute the exposure interval associated with each purchase. -Exposure intervals are stacked, that is each interval is put -after any previous. This means that the start of exposure to a given -purchase can be pushed into the future. The parameter \code{push.max} -indicates the maximally tolerated push. If this is reached by a person, -the assumption is that some of the purchased drug is not counted in the -exposure calculations. - -The \code{dpt} can either be a constant, basically translating the -purchased amount into exposure time the same way for all persons, or it -can be a vector with different treatment intensities for each purchase. -In any case the cumulative dose is computed taking this into account. + Exposure intervals are stacked, that is each interval is put after any + previous. This means that the start of exposure to a given purchase + can be pushed into the future. The parameter \code{push.max} indicates + the maximally tolerated push. If this is reached by a person, the + assumption is that some of the purchased drug is not counted in the + exposure calculations. + + The \code{dpt} can either be a constant, basically translating the + purchased amount into exposure time the same way for all persons, or + it can be a vector with different treatment intensities for each + purchase. In any case the cumulative dose is computed taking this + into account. If \code{use.dpt} is \code{FALSE} then the exposure from one purchase is assumed to stretch over the time to the next purchase, so we are effectively assuming different rates of dose per time between any two adjacent purchases. Moreover, with this approach, periods of - non-exposure does not exist. - -The intention of this function is to generate covariates for a -particular drug for the entire follow-up of each person. The reason that -the follow-up prior to drug purchase and post-exposure is included is -that the covariates must be defined for these periods too, in order to -be useful for analysis of disease outcomes. + non-exposure does not exist. Formally this approach conditions on the + future, because the rate of consumption (the accumulation of + cumulative exposure) is computed based on knowledge of when next + purchase is made. + + The intention of this function is to generate covariates for a + particular drug for the entire follow-up of each person. The reason + that the follow-up prior to drug purchase and post-exposure is + included is that the covariates must be defined for these periods too, + in order to be useful for analysis of disease outcomes. + + This function is described in terme of calendar time as underlying + time scale, because this will normally be the time scale for drug + purchases and for entry and exit for persons. In principle the + variables termed as dates might equally well refer to say the age + scale, but this would then have to be true \emph{both} for the + purchase data and the follow-up data. } \value{A data frame with one record per follow-up interval between @@ -94,7 +106,11 @@ \code{dof}.} } } -\author{Bendix Carstensen, \email{bxc@steno.dk}} + +\author{Bendix Carstensen, \email{bxc@steno.dk}. The development of +this function was supported partly through a grant from the EFSD +(European Foundation for the Study of Diabetes), ""} + \seealso{\code{\link{Lexis}}, \code{\link{splitLexis}}} \examples{ # Construct a simple data frame of purchases for 3 persons diff -Nru r-cran-epi-2.0/man/LCa.fit.Rd r-cran-epi-2.7/man/LCa.fit.Rd --- r-cran-epi-2.0/man/LCa.fit.Rd 1970-01-01 00:00:00.000000000 +0000 +++ r-cran-epi-2.7/man/LCa.fit.Rd 2016-10-01 11:23:35.000000000 +0000 @@ -0,0 +1,277 @@ +\name{LCa.fit} +\alias{LCa.fit} +\alias{print.LCa} +\alias{summary.LCa} +\alias{plot.LCa} +\alias{predict.LCa} +\title{ + Fit Lee-Carter-type models for rates to arbitrarily shaped observations + of rates in a Lexis diagram. +} +\description{ + The Lee-Carter model is originally defined as a model for rates + observed in A-sets (age by period) of a Lexis diagram, as + log(rate(x,t)) = a(x) + b(x)k(t), using one parameter per age(x) and + period(t). This function uses natural splines for a(), b() and k(), + placing knots for each effect such that the number of events is the + same between knots. Also fits the continuous time counterparts of all + models supported by the \code{\link[ilc]{lca.rh}} function from the + \code{ilc} package (see details). +} +\usage{ +LCa.fit( data, A, P, D, Y, + model = "APa", # or one of "ACa", "APaC", "APCa" or "APaCa" + 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 + c.ref, # cohort reference for the interactions + npar = c(a = 6, # no. knots for main age-effect + p = 6, # no. knots for peroid-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? + alpha = 0.05, # 1 minus confidence level + eps = 1e-6, # convergence criterion + maxit = 100, # max. no iterations + quiet = TRUE ) # cut the crap +\method{print}{LCa}( x, ... ) +\method{summary}{LCa}( object, show.est=FALSE, ... ) +\method{plot}{LCa}( x, ... ) +\method{predict}{LCa}( object, newdata, + alpha = 0.05, + level = 1-alpha, + sim = ( "vcov" \%in\% names(object) ), + ... ) + } +\arguments{ + \item{data}{ + A data frame. Must have columns \code{A}(age), \code{P}(period, that is + calendar time), \code{D}(no. of events) and \code{Y}(person-time, + exposure). Alternatively these four quantities can be given as + separate vectors: + } + \item{A}{ + Vector of ages (midpoint of observation). + } + \item{P}{ + Vector of period (midpoint of observation). + } + \item{D}{ + Vector of no. of events. + } + \item{Y}{ + Vector of person-time. Demographers would say "exposure", bewildering epidemiologists. + } + \item{a.ref}{ + Reference age for the age-interaction term(s) \code{pi(x)} and/or + \code{pi(x)}, where \code{pi(a.ref)=1} and \code{ci(a.ref)=1}. + } + \item{pi.ref}{ + Same, but specifically for the interaction with period. + } + \item{ci.ref}{ + Same, but specifically for the interaction with cohort. + } + \item{p.ref}{ + Reference period for the time-interaction term \code{kp(t)} where \code{kp(p.ref)=0}. + } + \item{c.ref}{ + 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. + } + \item{npar}{ + A (possibly named) vector or list, with either the number of knots or + the actual vectors of knots for each term. If unnamed, components are + taken to be in the order (a,b,t), if the model is "APaCa" in the order + (a,p,c,pi,ci). If a vector, the three integers indicate the number of + knots for each term; these will be placed so that there is an equal + number of events (\code{D}) between each, and half as many below the + first and above the last knot. If \code{npar} is a list of scalars the + behavior is the same. If \code{npar} is a list of vectors, these are + taken as the knots for the natural splines. + } + \item{VC}{ + Logical. Should the variance-covariance matrix of the parameters be + computed by numerical differentiation? See details. + } + \item{alpha}{ + 1 minus the confidence level used when calculating + confidence intervals for estimates in \code{LCa.fit} and for + predictions by \code{predict.LCa}. + } + \item{eps}{ + Convergence criterion for the deviance, we use the the relative + difference between deviance from the two models fitted. + } + \item{maxit}{ + Maximal number of iterations. + } + \item{quiet}{ + Shall I shup 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?} + \item{x}{An \code{LCa} object, see under "Value".} + \item{newdata}{Prediction data frame, must have columns \code{A} and + \code{P}. Any \code{Y} column is ignored, predictions are given in + units of the \code{Y} supplied for the call that generated the + \code{LCa} object.} + \item{level}{Confidence level.} + \item{sim}{Logical or numeric. If \code{TRUE}, prediction c.i.s will be + based on 1000 simulations from the posterior parameters. If numeric, + it will be based on that number of simulations.} + \item{...}{Additional parameters passed on to the method.} +} +\details{ + The Lee-Carter model is non-linear in age and time so does not fit + in the classical glm-Poisson framework. But for fixed \code{b(x)} it + is a glm, and also for fixed \code{a(x)}, \code{k(t)}. The function + alternately fits the two versions until the same fit is produced (same + deviance). + + The multiplicative age by period term could equally well have been a + 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)} + + 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 + \code{a} if there is an interaction with age. Thus there are 5 different + models that can be fitted: \code{APa}, \code{ACa}, \code{APaC} \code{APCa} + and \code{APaCa}. + + The standard errors of the parameters from the two model fits are however + wrong; they are conditional on some of terms having a fixed value. And + the symbolic calculation of the Hessian is a nightmare, so this is done + numerically using the \code{hessian} function from the \code{numDeriv} + package if \code{VC=TRUE}. + + 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 + 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. + + The \code{plot} for \code{LCa} objects merely produces between 3 and 5 + panels showing each of the terms in the model. These are mainly for + preliminary inspection; real reporting of the effects should use + proper relative scaling of the effects.} + +\value{ + \code{LCa.fit} returns an object of class \code{LCa} (smooth + effects \code{L}ee-\code{Ca}rter model); it is a list with the + following components: +\item{model}{Character, either \code{APa}, \code{ACa}, \code{APaC}, + \code{APCa} or \code{APaCa}, indicating the variable(s) interacting + with age.} +\item{ax}{3-column matrix of age-effects, c.i. from the age-time + model. Rownames are the unique occurring ages in the + dataset. Estimates are rates.} +\item{pi}{3-column matrix of age-period interaction effects, c.i. from the age + model. Rownames are the actually occurring ages in the + dataset. Estimates are multipliers of the log-RRs in \code{kt}, + centered at 1 at \code{ci.ref}.} +\item{kp}{3-column matrix of period-effects, with c.i.s from the + age-time model. Rownames are the actually occurring times in the + dataset. Estimates are rate-ratios centered at 1 at \code{p.ref}.} +\item{ci}{3-column matrix of age-cohort interaction effects, c.i. from the age + model. Rownames are the actually occurring ages in the + dataset. Estimates are multipliers of the log-RRs in \code{kt}, + centered at 1 at \code{ci.ref}.} +\item{kc}{3-column matrix of period-effects, with c.i.s from the age-time + model. Rownames are the actually occurring times in the + dataset. Estimates are rate-ratios centered at 1 at \code{p.ref}.} +\item{mod.at}{\code{glm} object with the final age-time model. Gives + the same fit as the \code{mod.b} model.} +\item{mod.b}{\code{glm} object with the final age model. Gives + the same fit as the \code{mod.at} model.} +\item{coef}{All coefficients from both models, in the order \code{ax}, + \code{kp}, \code{kc}, \code{pi}, \code{ci}. Only present if + \code{LCa.fit} were called with + \code{VC=TRUE} (the default).} +\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 + cohort effects.} +\item{refs}{List of reference points used for the age, period and + cohort terms in the interactions.} +\item{deviance}{Deviance of the model} +\item{df.residual}{Residual degrees of freedom} +\item{iter}{Number of iterations used to reach convergence.} + + \code{plot.LCa} plots the estimated effects in separate panels, + using a log-scale for the baseline rates (\code{ax}) and the time-RR + (\code{kt}). For the \code{APaCa} model 5 panels are plotted. + + \code{summary.LCa} returns (invisibly) a matrix with the parameters + from the models and a column of the conditional se.s and of the se.s + derived from the numerically computed Hessian (if \code{LCa.fit} were + called with \code{VC=TRUE}.) + + \code{predict.LCa} returns a matrix with one row per row in + \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 + \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} +} +\seealso{ + \code{\link{apc.fit}}, + \code{\link{apc.LCa}}, + \code{\link[ilc]{lca.rh}}, + \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 +data( testisDK ) +tc <- subset( testisDK, A>14 & A<60 ) +head( tc ) + +# We want to see rates per 100,000 PY +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 +summary( LCa.tc ) + +# Inspect what we got +names( LCa.tc ) + +# show the estimated effects +par( mfrow=c(1,3) ) +plot( LCa.tc ) + +# Prediction data frame for ages 15 to 60 for three time points: +nd <- data.frame( A=15:60 ) + +p50 <- predict.LCa( LCa.tc, newdata=cbind(nd,P=1950), sim=10000 ) +p70 <- predict.LCa( LCa.tc, newdata=cbind(nd,P=1970), sim=10000 ) +p90 <- predict.LCa( LCa.tc, newdata=cbind(nd,P=1990), sim=10000 ) + +# Inspect the curves from the parametric bootstrap (simulation): +par( mfrow=c(1,1) ) +matplot( nd$A, cbind(p50,p70,p90), type="l", lwd=c(6,3,3), lty=1, + col=rep( c("red","green","blue"), each=3 ), log="y", + ylab="Testis cancer incidence per 100,000 PY, 1970, 80, 90", xlab="Age" ) +} +\keyword{models} +\keyword{regression} diff -Nru r-cran-epi-2.0/man/lls.Rd r-cran-epi-2.7/man/lls.Rd --- r-cran-epi-2.0/man/lls.Rd 2015-05-27 15:27:03.000000000 +0000 +++ r-cran-epi-2.7/man/lls.Rd 2016-10-01 08:20:42.000000000 +0000 @@ -44,8 +44,13 @@ m1 <- glm( y ~ x, family=binomial ) M <- matrix( 1:20, 4, 5 ) .M <- M +dfr <- data.frame(x,y) +attach( dfr ) lls() +search() clear() +search() lls() +lls(all=TRUE) } \keyword{attributes} \ No newline at end of file diff -Nru r-cran-epi-2.0/man/Ns.Rd r-cran-epi-2.7/man/Ns.Rd --- r-cran-epi-2.0/man/Ns.Rd 2016-01-01 14:20:50.000000000 +0000 +++ r-cran-epi-2.7/man/Ns.Rd 2016-04-21 13:16:04.000000000 +0000 @@ -2,8 +2,8 @@ \alias{Ns} \title{ Natural splines - (cubic splines linear beyond outermost knots) with - convenient specification of knots and possibility of centering and - detrending. + convenient specification of knots and possibility of centering, + detrending and clamping. } \description{ This function is partly for convenient specification of natural splines @@ -11,13 +11,15 @@ and the largest of the supplied knots as boundary knots. It also has the option of centering the effects provided at a chosen reference point as well as projecting the columns on the orthogonal space to - that spanned by the intercept and the linear effect of the variable. + that spanned by the intercept and the linear effect of the variable, + and finally fixing slopes beyond boundary knots (clamping). } \usage{ Ns( x, ref = NULL, df = NULL, knots = NULL, intercept = FALSE, Boundary.knots = NULL, + fixsl = c(FALSE,FALSE), detrend = FALSE ) } \arguments{ @@ -31,7 +33,11 @@ \item{intercept}{Should the intercept be included in the resulting basis? Ignored if any of \code{ref} or \code{detrend} is given.} \item{Boundary.knots}{The boundary knots beyond which the spline is - linear.} + linear. Defaults to the minimum and maximum of \code{knots}.} + \item{fixsl}{Specification of whether slopes beyond outer knots should + be fixed to 0. \code{FALSE} correponds to no restriction; a curve with 0 + slope beyond the upper knot is obtained using + \code{c(FALSE,TRUE)}. Ignored if \code{!(detrend==FALSE)}.} \item{detrend}{If \code{TRUE}, the columns of the spline basis will be projected to the orthogonal of \code{cbind(1,x)}. Optionally \code{detrend} can be given as a vector of non-negative numbers used @@ -49,16 +55,18 @@ \code{diag(weight)} (an assumption which is checked). } \author{ - Bendix Carstensen + Bendix Carstensen \email{bxc@steno.dk}, + Lars Jorge D\'iaz, Steno Diabetes Center. } \note{ The need for this function is primarily from analysis of rates in epidemiology and demography, where the dataset are time-split records of follow-up, and the range of data therefore rarely is of any - interest (let alone meaningful). In Poisson modeling of rates - based on time-split records one should aim at having the same number - of \emph{events} between knots, rather than the same number of - observations. + interest (let alone meaningful). + + In Poisson modeling of rates based on time-split records one should + aim at having the same number of \emph{events} between knots, rather + than the same number of observations. } \examples{ require(splines) @@ -91,18 +99,49 @@ ( a.kn <- with( subset(dms,lex.Xst=="Dead"), quantile( Age+lex.dur, probs=(1:n.kn-0.5)/n.kn ) ) ) n.kn <- 4 -( p.kn <- with( subset(dms,lex.Xst=="Dead"), +( p.kn <- with( subset( dms, lex.Xst=="Dead" ), quantile( Per+lex.dur, probs=(1:n.kn-0.5)/n.kn ) ) ) m1 <- glm( lex.Xst=="Dead" ~ Ns( Age, kn=a.kn ) + Ns( Per, kn=p.kn, ref=2000 ), - offset = log( lex.dur ), family=poisson, data=dms ) + offset = log( lex.dur ), + family = poisson, + data = dms ) # Plot estimated age-mortality curve for the year 2005 and knots chosen: -nd <- data.frame(Age=40:90,Per=2005,lex.dur=1000) +nd <- data.frame( Age=seq(40,100,0.1), Per=2005, lex.dur=1000 ) par( mfrow=c(1,2) ) matplot( nd$Age, ci.pred( m1, newdata=nd ), type="l", lwd=c(3,1,1), lty=1, col="black", log="y", - ylab="Mortality rates per 1000 PY", xlab="Age (years)", las=1 ) + ylab="Mortality rates per 1000 PY", xlab="Age (years)", las=1, ylim=c(1,1000) ) +rug( a.kn, lwd=2 ) + +# Clamped Age effect to the right of rightmost knot. +m1.c <- glm( lex.Xst=="Dead" ~ Ns( Age, kn=a.kn, fixsl=c(FALSE,TRUE) ) + + Ns( Per, kn=p.kn, ref=2000 ), + offset = log( lex.dur ), + family = poisson, + data = dms ) + +# Plot estimated age-mortality curve for the year 2005 and knots chosen. +matplot( nd$Age, ci.pred( m1.c, newdata=nd ), + type="l", lwd=c(3,1,1), lty=1, col="black", log="y", + ylab="Mortality rates per 1000 PY", xlab="Age (years)", las=1, ylim=c(1,1000) ) +rug( a.kn, lwd=2 ) + +par( mfrow=c(1,1) ) + +# Including a linear Age effect of 0.05 to the right of rightmost knot. +m1.l <- glm( lex.Xst=="Dead" ~ Ns( Age, kn=a.kn, fixsl=c(FALSE,TRUE) ) + + Ns( Per, kn=p.kn, ref=2000 ), + offset = log( lex.dur ) + pmax( Age, max( a.kn ) ) * 0.05, + family = poisson, + data = dms ) + +# Plot estimated age-mortality curve for the year 2005 and knots chosen. +nd <- data.frame(Age=40:100,Per=2005,lex.dur=1000) +matplot( nd$Age, ci.pred( m1.l, newdata=nd ), + type="l", lwd=c(3,1,1), lty=1, col="black", log="y", + ylab="Mortality rates per 1000 PY", xlab="Age (years)", las=1, ylim=c(1,1000) ) rug( a.kn, lwd=2 ) } \keyword{regression} diff -Nru r-cran-epi-2.0/man/Relevel.Rd r-cran-epi-2.7/man/Relevel.Rd --- r-cran-epi-2.0/man/Relevel.Rd 2015-05-27 15:27:03.000000000 +0000 +++ r-cran-epi-2.7/man/Relevel.Rd 2016-07-20 12:40:11.000000000 +0000 @@ -24,6 +24,7 @@ An unordered factor, where levels of \code{x} have been reordered and/or collapsed. } +\author{Bendix Carstensen, \url{BendixCarstensen.com}.} \seealso{\code{\link{Relevel.Lexis}}} \examples{ ff <- factor( sample( letters[1:5], 100, replace=TRUE ) ) diff -Nru r-cran-epi-2.0/man/rm.tr.Rd r-cran-epi-2.7/man/rm.tr.Rd --- r-cran-epi-2.0/man/rm.tr.Rd 1970-01-01 00:00:00.000000000 +0000 +++ r-cran-epi-2.7/man/rm.tr.Rd 2016-10-01 08:54:08.000000000 +0000 @@ -0,0 +1,67 @@ +\name{rm.tr} +\alias{rm.tr} +\title{ +Remove transitions from a Lexis object. +} +\description{ +Sometimes certain transitions are not of interest. This function removes +these and assigns the risk time in the target state of the transitions +to the originating state. +} +\usage{ +rm.tr(obj, from, to) +} +\arguments{ + \item{obj}{ +A \code{Lexis} object. +} + \item{from}{ +Character; name of the state from which the transition to be purged +originates. Must be a valid state name for \code{obj}. +} + \item{to}{ +Character; name of the state to which the transition to be purged +targets. Must be a valid state name for \code{obj}. +} +} +\details{ +The function removes all transitions from \code{from} to \code{to}, and +assigns all risk time in the \code{to} state after the transition +(\code{lex.dur}) to the \code{from} state. This is only done for risk +time in \code{to} occurring directly after \code{from}. Risk time in +\code{to} occurring after a transition from states different from +\code{from} is not affected. Transitions from \code{to} to another +state, \code{other}, say, will be changed to transitions from +\code{from} to \code{other}. +} +\value{ +A \code{\link{Lexis}} object with the indicated transition removed. +} +\author{ +Bendix Carstensen, \url{BendixCarstensen.com}. +} +\seealso{ +\code{\link{Relevel}} +} +\examples{ +data(DMlate) +dml <- Lexis( entry = list(Per=dodm, Age=dodm-dobth, DMdur=0 ), + exit = list(Per=dox), + exit.status = factor(!is.na(dodth),labels=c("DM","Dead")), + data = DMlate ) + +# A small subset for illustration +dml <- subset( dml, lex.id \%in\% c(13,15,20,28,40) ) + +# Cut the follow-up at start of insulin therapy +dmi <- cutLexis( dml, cut = dml$doins, + pre = "DM", + new.state = "Ins" )[,1:10] + +# How does it look? +dmi + +# Remove all transitions DM -> Ins +rm.tr( dmi, "DM", "Ins" ) +} +\keyword{manip} diff -Nru r-cran-epi-2.0/man/subset.Lexis.Rd r-cran-epi-2.7/man/subset.Lexis.Rd --- r-cran-epi-2.0/man/subset.Lexis.Rd 2016-01-04 13:51:22.000000000 +0000 +++ r-cran-epi-2.7/man/subset.Lexis.Rd 2016-02-17 15:45:16.000000000 +0000 @@ -21,7 +21,7 @@ } \details{ The subset method for \code{Lexis} objects works exactly as the method - for data frames. So does the "[" method. The special methods are neeed in + for data frames. So does the "[" method. The special methods are needed in order to propagate the Lexis-specific attributes. The method for \code{stacked.Lexis} objects also shrinks the set of diff -Nru r-cran-epi-2.0/man/summary.Lexis.rd r-cran-epi-2.7/man/summary.Lexis.rd --- r-cran-epi-2.0/man/summary.Lexis.rd 2015-05-27 15:27:03.000000000 +0000 +++ r-cran-epi-2.7/man/summary.Lexis.rd 2016-10-01 09:32:27.000000000 +0000 @@ -9,7 +9,8 @@ (\code{lex.Cst} and \code{lex.Xst}), as well the risk time in each state. } \usage{ - \method{summary}{Lexis}( object, simplify=TRUE, scale=1, by=NULL, Rates=FALSE, ... ) + \method{summary}{Lexis}( object, simplify=TRUE, scale=1, by=NULL, + Rates=FALSE, timeScales=FALSE, ... ) \method{print}{summary.Lexis}( x, ..., digits=2 ) } \arguments{ @@ -25,6 +26,8 @@ values of this will be used to stratify the summary.} \item{Rates}{Should a component with transition rates be returned (and printed) too?} +\item{timeScales}{Should the names of the timescales and the indication + of since which entry also be given?} \item{x}{A \code{summary.Lexis} object.} \item{digits}{How many digits should be used for printing?} \item{ ... }{Other parameters - ignored} diff -Nru r-cran-epi-2.0/MD5 r-cran-epi-2.7/MD5 --- r-cran-epi-2.0/MD5 2016-01-06 19:29:29.000000000 +0000 +++ r-cran-epi-2.7/MD5 2016-10-04 12:35:50.000000000 +0000 @@ -1,32 +1,34 @@ -ff50971550f0c63fd0986919ea174c08 *CHANGES -b88050b7250d1b5eb37d62c2d67ca04b *DESCRIPTION -2783172f6bf5e6c89f86771d3e8760b4 *NAMESPACE +203c582ea761fd50ca3bb7549f87f920 *CHANGES +b812b1d33d655d933f3c57e067a687f9 *DESCRIPTION +98ff4454ab2bc144d016e0cf6a00b8a0 *NAMESPACE fac3c7f01ab0bd6930cf3cbea20e6bcc *R/Aplot.R -75d172080863bd6e3b2286a0c94b0d34 *R/Cplot.R +a031997be7388601dd1df6a28fc55328 *R/Cplot.R 6a1c4f4cbbf509d1f3c40b9f0adb6399 *R/Icens.R +ecdc4193570616dd8b42ca451941fb29 *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 -1d3f55ffad7833ae4615bdde4c7308c2 *R/Ns.r +9a46235d94dabfbf52b324083f573492 *R/Ns.r e9f582b5dc17d07bf2df96de74441963 *R/Pplot.R b1c0f04443499205645ad578a3dd8573 *R/ROC.R 0a7bbd37e75920b54fd0fbe00a86511d *R/Relevel.R 0e3e3d23b8ab79dcf3046828190462b2 *R/Termplot.R 28a0ec0c8584e22759b516ff69e5bd78 *R/Wald.R -563f0b2d90a3776765207cc61a31e134 *R/apc.fit.R +b1f0000430d6e15f59fe949ed8f1beb1 *R/apc.LCa.R +847f61901b2c9123367be3ebef30abe6 *R/apc.fit.R 3c53b3d667196c204fab811d2231ccd5 *R/apc.frame.R 40d04e6ce1304d51df6911621fe28b52 *R/apc.plot.R 8ce5a4c989947cb9d93becdecfde642c *R/as.Date.cal.yr.R -c231a1067a4cc2fc9af0dc4b0d2c6792 *R/boxes.MS.R -11d3966379eacdc164f2d50c4da8bb70 *R/cal.yr.R +f4e2a039dda786e4cad0b2961e3dd599 *R/boxes.MS.R +44e5f888817c94b8a190eae4aa964a2c *R/cal.yr.R 6988e76eb64f004a8c5029c23f8603af *R/ccwc.R d5171eb23f0b431a7cc89b69d03308a6 *R/ci.cum.R 1e69aa59a2367037c3a054bfa1a2759e *R/ci.lin.R 6f80c7b2b2120068cd16c043fbc29db3 *R/ci.mat.R 886f87cbd08266fa281e2ed7b3b3c253 *R/ci.pd.R -acab8b478be5a44a36561e65d70f2b78 *R/clear.R +7637da11fdfb3612df5d73a64ee044e4 *R/clear.R 89f2fe2a79bba0565efc02fd4aaeca9a *R/clogistic.R d365bc4739ebb104b0318c1c9ee226fd *R/contr.2nd.R bcd4142d3bf3e4517397387fe43fdf95 *R/contr.cum.R @@ -37,6 +39,7 @@ 14270f90ed00f4f043286d50e2de9f26 *R/detrend.R a86250cce524f9fa2215c0c0cac4c6e7 *R/effx.match.r 9b46438b1c6fdc926ac1a253a85096aa *R/effx.r +fac08dcafcfddfe8e9857f2546b5f0ea *R/erl.R 822e02862de39625e5ed9ef6f84bdf39 *R/expand.data.r 0466dba42d10766e0b9cd220ad39dcf4 *R/factorize.R 67ce5a28408f42c7b65084195897b995 *R/fit.add.r @@ -46,7 +49,7 @@ 23a68a19f8f8e70c5664eb5110b3e807 *R/foreign.R da7b5d41fe7876ef6bab01839b613939 *R/ftrend.R f0f84fc26294b1cb3016a9e23974c545 *R/gen.exp.R -6158aa76619a4ee8b0f6e1d5535fdbe7 *R/lexis.R +52ccaf61645561d9c050e5d5665fbdcd *R/lexis.R 9c575049aa6c3841e7e92f563a6fc450 *R/lls.R 79d6e0b02a2c4e98fbbfcd936b21d377 *R/mh.R 9e48e47769874a34f63a2e2d4e6ffffc *R/ncut.r @@ -57,15 +60,16 @@ cca51283f637927702fbe217367011eb *R/print.floated.R 961858aa39320463dffa6154c70c87f6 *R/projection.ip.r 0c2a2e493c190c4e74587b5127e9bef6 *R/rateplot.R +251d9ec19ef0ca158582f78d70255c23 *R/rm.tr.R 804e11d9d564a860fce275258f40db3d *R/simLexis.R 358b6fa60093d8bbee969a76a88797c0 *R/splitLexis.R e08261ab194f4aca8502ee249afd309e *R/stack.Lexis.R fb143d532001643d3f74fbd2e114e976 *R/stattable.R 889282ea37a5a9d9bf7155d790cd9036 *R/summary.Icens.r -ce8c23fd028541b00472f792246e39de *R/summary.Lexis.r +8ceccd0ffa9694670f38c34b65d5be10 *R/summary.Lexis.r 88983ba8314a37910993d469a982d800 *R/twoby2.R 497daef89f1ea705cd2c755132bec69d *R/zzz.R -917367e2fbe5a2d1f271222936fcb8a6 *build/vignette.rds +19da589577365d1c7f87fa6460dfefdc *build/vignette.rds b9d3112271b9545896c98ec6406f5f02 *data/B.dk.rda 666e336a1aefb9c4298abb83a81a3300 *data/DMconv.rda e200a06eee9305b2fa2c59af5982d62b *data/DMlate.rda @@ -90,20 +94,19 @@ 4640a31484934cccc1d626104af1365f *data/occup.rda 404cc2c826ccdfb5b5edbe9834561e02 *data/testisDK.rda 1953ddd7d750051c5682fa08b11cb777 *data/thoro.rda -e0f7d0bbeb4f050e41a55dd412548816 *inst/CITATION -acae1b4099b6e046c314f51aaa0327ce *inst/doc/Follow-up.R -d9c94a943d90f4e94e08c8386779143c *inst/doc/Follow-up.pdf +9bd39178a387f935acc54519109c9f4a *inst/CITATION +e61ded63f726562d6f3172838e4f7337 *inst/doc/Follow-up.R +f6e45b6606e4c70ae545b651bc3ea59b *inst/doc/Follow-up.pdf a6978ca7d45a578a4fb047fb92eeed73 *inst/doc/Follow-up.rnw -aad0285b496c6f1d8385a41ca441158c *inst/doc/index.html -82484f2ce4b2c1dfcde893b123ee9a9e *inst/doc/sim-Lexis.R -0f26b30e8b1f919e01203435d2aeac35 *inst/doc/sim-Lexis.pdf -d9264ec0daa372a8921ae4ce12257eba *inst/doc/simLexis.R -2e393710d9af312888df236d7efe1d23 *inst/doc/simLexis.pdf -0e75c4f99139b36ab0d52f28ee80e3fe *inst/doc/simLexis.rnw +fff6643f55321f70c23cb36b9ee587e0 *inst/doc/index.html +a9ec600c3c54c2b2e0a252f599bda85a *inst/doc/simLexis.R +6ccba2ccef7004dd6f2c129f0e80afce *inst/doc/simLexis.pdf +c38c581cd2ac1acdcce9e60e56dc37ca *inst/doc/simLexis.rnw 317bb6ebeaec5e7376b609d7e382c166 *man/B.dk.Rd 3abc5c20e62c874d70dcb2688ad2cce5 *man/DMconv.Rd 362200b221355a7bcbcd1c7c8df5243f *man/DMlate.Rd 38f9ae52c47f70d01f46eb437fca0f16 *man/Icens.Rd +39d7b1fc6a514f8a6fd7d5013893a4c5 *man/LCa.fit.Rd 1ab9326a322e35e9855be902cc2a5da8 *man/Lexis.Rd c97dabafce248e4d47c5b18b6f04efe0 *man/Lexis.diagram.Rd 2e71ca58213788c8a358708e776e2112 *man/Lexis.lines.Rd @@ -112,24 +115,25 @@ 6ae7494a426801b9944c37190c09b3ee *man/N.dk.Rd 0983e4031974eea1a30e67db90310c25 *man/N2Y.Rd 424aecfdd664b3123ad4815880aa3761 *man/NArray.Rd -f8cd42bd595dbab38f87e5834aeb092a *man/Ns.Rd +5d1d5dcea53632bf7ea6f18adb174014 *man/Ns.Rd 335b75739b8a25293690da804cd5e522 *man/ROC.Rd -5fa418e7a462df9892e7781ef5027412 *man/Relevel.Rd +ca827755eae8ee8286860a715cbd1744 *man/Relevel.Rd bb2162d557dfb7685edac64ac3213b6c *man/S.typh.Rd 9351888a37017068f2c9ac8b843c846d *man/Termplot.Rd b4879df831de32e9e298b7b3d90389a9 *man/Y.dk.Rd -d99521bbb38bd98eac3acebf9b83582c *man/apc.fit.Rd +62c36534ab2b6043f68306dd33407891 *man/apc.LCa.Rd +d5976837ba976d5969aaaedb98b215eb *man/apc.fit.Rd 8e42054ad19c34f2bf77210cfafc560e *man/apc.frame.Rd 92b663829915395cdabf3e5d37cdd2b4 *man/bdendo.Rd 968ba8cef65153e9c9c521dd1541123e *man/bdendo11.Rd d74a4f6e66ce91c73cd2817bb690446a *man/births.Rd 8dc11840302a388e4f26892335d5d888 *man/blcaIT.Rd -63c1de8406186f3a692a7091a0bba2a7 *man/boxes.MS.Rd +6b9f6bc7ca862808b89f6439914b8a0a *man/boxes.MS.Rd 18d942d375c5d5f6899291bf9cd1fe11 *man/brv.Rd a1a09c9591d54235cc026e44ca4ddad8 *man/cal.yr.Rd effba2c17ade00f0afd10184a073d34c *man/ccwc.Rd c47967c97d376b54706caacce489a9a2 *man/ci.cum.Rd -909d557dcade1d5adb6ceda22712b00e *man/ci.lin.Rd +fd00c162c3b545d86546615dfc24750f *man/ci.lin.Rd e7477402aa4746d653a0c72b2239d185 *man/ci.pd.Rd 128c0b1d00a69e41d131a6f74a1bae23 *man/clogistic.Rd cf5e5a6e18e07b7bb8d70dc0417e6166 *man/contr.cum.Rd @@ -139,20 +143,21 @@ 386d0a4aa3515cf754ef458e1f32089c *man/diet.Rd 64dd6599bb830e98f4717b1a34ee12ab *man/effx.Rd 3c3cc3b1ecabece6b7709c1b147bab4e *man/effx.match.Rd +139ed9b0278b94d417e9464472c3ef5c *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 bca8b8c3452180e63d4d22885b04e3ce *man/float.Rd -a00c7d3fce9cf363ccb950996276da76 *man/foreign.Lexis.Rd +c96cdd9ce1f87335dd8eb454341f44ba *man/foreign.Lexis.Rd e5a57a203cb91e9680dedc9023779e59 *man/ftrend.Rd -0b34f125847915dde69f22e1b3759a3c *man/gen.exp.Rd +9a1a04510b887f39bab1c2db4b6d2cfd *man/gen.exp.Rd cbbc2a23f902d83ba2527b27b5ed3adb *man/gmortDK.Rd d9f1d31e109d6058a2160ca0aa49bf81 *man/hivDK.Rd 9c3b059921d728298ccff73298b4dd31 *man/lep.Rd ec0e27968d4ee35f390161e49abd903c *man/lines.apc.Rd -c17abc408d69d11bb660d92788056f77 *man/lls.Rd +78125ba49d54905fb727765dc8997796 *man/lls.Rd 661ec298ffcd0f39cf55c5d033e04fd1 *man/lungDK.Rd 3c462c4ce7ca1ce4b8c6fef50b3c7623 *man/merge.Lexis.Rd 4c26c9edf5d64b0af3e63ad791b2a280 *man/merge.data.frame.Rd @@ -171,14 +176,15 @@ 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 eaad8ae22bd68a18750a94352f65f9c3 *man/start.Lexis.Rd 4f43436fc863747dad4fcdd566948039 *man/stattable.Rd 73cf4597fcb4a9d019bda47431b095be *man/stattable.funs.Rd -19c524bd545956a7ca54b22838edb660 *man/subset.Lexis.Rd -b89f3f488011bd32983b467915dd4813 *man/summary.Lexis.rd +f10fc89dadf571e07f3b13923660cff6 *man/subset.Lexis.Rd +2639c132d73f56d51bca04d8d5e22b66 *man/summary.Lexis.rd 14711d23ae0b89d9d22b9f3b07211093 *man/testisDK.Rd e3446824b09ea4ad8e8a0ddfc1987eef *man/thoro.Rd 6ea9cf64475187fa2f402275c40573a7 *man/time.band.Rd @@ -190,5 +196,5 @@ 33be867717f019d03821e58b662a94c4 *src/chsolve2.c ce9e73f4028463cf8b84a7783f282ac3 *src/clogit.c a6978ca7d45a578a4fb047fb92eeed73 *vignettes/Follow-up.rnw -aad0285b496c6f1d8385a41ca441158c *vignettes/index.html -0e75c4f99139b36ab0d52f28ee80e3fe *vignettes/simLexis.rnw +fff6643f55321f70c23cb36b9ee587e0 *vignettes/index.html +c38c581cd2ac1acdcce9e60e56dc37ca *vignettes/simLexis.rnw diff -Nru r-cran-epi-2.0/NAMESPACE r-cran-epi-2.7/NAMESPACE --- r-cran-epi-2.0/NAMESPACE 2016-01-04 14:00:42.000000000 +0000 +++ r-cran-epi-2.7/NAMESPACE 2016-10-03 10:15:49.000000000 +0000 @@ -27,6 +27,7 @@ contr.diff, detrend, dur, + erl, surv1, surv2, erl1, yll, effx, effx.match, float, @@ -38,9 +39,16 @@ fit.add, fit.mult, ftrend, + show.apc.LCa, + apc.LCa, + LCa.fit, + print.LCa, + summary.LCa, + predict.LCa, + plot.LCa, Lexis.diagram, Lexis.lines, - Life.lines, + Life.lines, Lexis, merge.Lexis, plot.Lexis, @@ -48,6 +56,7 @@ lines.Lexis, PY.ann.Lexis, subset.Lexis, + "[.Lexis", cbind.Lexis, rbind.Lexis, summary.Lexis, @@ -85,6 +94,7 @@ boxarr, boxes, factorize, + rm.tr, PY.ann, N2Y, tmat, @@ -121,8 +131,11 @@ importFrom( plyr, rbind.fill ) importFrom( cmprsk, crr) importFrom( etm, etm ) +# importFrom( mstate, msdata ) importFrom( MASS, mvrnorm ) -importFrom( survival, clogit) +importFrom( survival, clogit ) +importFrom( numDeriv, hessian ) +importFrom( Matrix, nearPD ) importFrom("grDevices", "gray", "rainbow") importFrom("graphics", "abline", "arrows", "axis", "box", "layout", "lines", "locator", "matlines", "matplot", "matpoints", @@ -137,7 +150,6 @@ "poisson", "predict", "qnorm", "qt", "quantile", "runif", "termplot", "update", "vcov", "weighted.mean", "xtabs") importFrom("utils", "flush.console", "str") -## importFrom( mstate, msdata ) # register S3 methods S3method( print, Icens) @@ -153,7 +165,7 @@ S3method( merge, Lexis) S3method( subset, Lexis) S3method( subset, stacked.Lexis) -# S3method( "[", Lexis) +S3method( "[", Lexis) S3method( cbind, Lexis) S3method( rbind, Lexis) S3method( summary, Lexis) @@ -172,6 +184,10 @@ S3method( boxes, MS) S3method( msdata, Lexis) S3method( etm, Lexis) +S3method( print, LCa) +S3method( summary, LCa) +S3method( predict, LCa) +S3method( plot, LCa) S3method( merge, data.frame) S3method( print, stat.table) S3method( print, clogistic) diff -Nru r-cran-epi-2.0/R/apc.fit.R r-cran-epi-2.7/R/apc.fit.R --- r-cran-epi-2.0/R/apc.fit.R 2015-05-28 08:04:38.000000000 +0000 +++ r-cran-epi-2.7/R/apc.fit.R 2016-05-24 22:35:38.000000000 +0000 @@ -8,17 +8,17 @@ ref.p, dist = c("poisson","binomial"), model = c("ns","bs","ls","factor"), - dr.extr = c("weighted","Holford"), + dr.extr = "weighted", parm = c("ACP","APC","AdCP","AdPC","Ad-P-C","Ad-C-P","AC-P","AP-C"), npar = c( A=5, P=5, C=5 ), scale = 1, alpha = 0.05, print.AOV = TRUE ) { -dist <- match.arg(dist) + dist <- match.arg(dist) model <- match.arg(model) -drtyp <- match.arg(dr.extr) -parm <- toupper(match.arg(parm)) +drtyp <- deparse(substitute(dr.extr)) + parm <- toupper(match.arg(parm)) if(!missing(data)) { if (length(match(c("A", "P", "D", "Y"), names(data))) != 4) @@ -115,16 +115,23 @@ Rc <- MC[abs(P - A - c0) == min(abs(P - A - c0)), , drop = FALSE][1, ] } if (model == "ns") { - knl <- is.list(npar) - MA <- if (knl) Ns(A, knots = npar[["A"]] ) - else Ns(A, knots = quantile( rep(A,D), - probs=(0:npar["A"]+0.1)/(npar["A"]+0.2) ) ) - MP <- if (knl) Ns(P, knots = npar[["P"]] ) - else Ns(P, knots = quantile( rep(P,D), - probs=(0:npar["P"]+0.1)/(npar["P"]+0.2) ) ) - MC <- if (knl) Ns(P - A, knots = npar[["C"]] ) - else Ns(P-A, knots = quantile( rep(P-A,D), - probs=(0:npar["C"]+0.1)/(npar["C"]+0.2) ) ) + # is npar a list + knl <- is.list( npar ) + # if scalar expand + if( !knl & length(npar)==1 ) npar <- rep( npar, 3 ) + # if no names, provide them + if( is.null(names(npar)) ) names(npar) <- c("A","P","C") + # if names too long or wrong case, rectify + names( npar ) <- toupper( substr(names(npar),1,1) ) + MA <- if (knl) Ns( A, knots = npar[["A"]] ) + else Ns( A, knots = quantile( rep(A,D), + probs=(1:npar["A"]-0.5)/npar["A"] ) ) + MP <- if (knl) Ns(P , knots = npar[["P"]] ) + else Ns(P , knots = quantile( rep(P,D), + probs=(1:npar["P"]-0.5)/npar["P"] ) ) + MC <- if (knl) Ns(P-A, knots = npar[["C"]] ) + else Ns(P-A, knots = quantile( rep(P-A,D), + probs=(1:npar["C"]-0.5)/npar["C"] ) ) Rp <- ns(p0, knots = attr(MP,"knots"), Boundary.knots = attr(MP,"Boundary.knots")) Rc <- ns(c0, knots = attr(MC,"knots"), @@ -147,7 +154,7 @@ Boundary.knots = npar[["P"]][ c(1,nk[2])], degree = deg) else bs(P, df = npar[["P"]], degree = deg) MC <- if (knl) bs(P - A, knots = npar[["C"]][-c(1, nk[3])], - Boundary.knots = npar[["C"]][c(1, nk[3])], degree = deg) + Boundary.knots = npar[["C"]][ c(1, nk[3])], degree = deg) else bs(P - A, df = npar[["C"]], degree = deg) Rp <- bs(p0, knots = attr(MP,"knots"), Boundary.knots = attr(MP,"Boundary.knots"), @@ -186,14 +193,23 @@ C.pt <- unique(P - A) C.pos <- match(C.pt, P - A) MA <- cbind(1, MA) -if (!mode(drtyp) %in% c("character", "numeric")) +if (!mode(dr.extr) %in% c("character", "numeric")) stop("\"dr.extr\" must be of mode \"character\" or \"numeric\".\n") -if (is.character(drtyp)) - wt <- if (toupper(substr(drtyp, 1, 1)) == "W") - D - else rep(1, length(D)) -if (is.numeric(drtyp)) - wt <- drtyp +if (is.character(dr.extr)) + { + wt <- rep(1, length(D) ) + drtyp <- "1-weights" + if( toupper(substr(dr.extr, 1, 1)) %in% c("W","T","D") ) + { wt <- D + drtyp <- "D-weights" } + if( toupper(substr(dr.extr, 1, 1)) %in% c("L","R") ) + { wt <- (Y^2)/D + drtyp <- "Y2/D-weights" } + if( toupper(substr(dr.extr, 1, 1)) %in% c("Y") ) + { wt <- Y + drtyp <- "Y-weights" } + } +if ( is.numeric(dr.extr) ) wt <- dr.extr Rp <- matrix(Rp, nrow = 1) Rc <- matrix(Rc, nrow = 1) xP <- detrend(rbind(Rp, MP), c(p0, P), weight = c(0, wt)) @@ -206,7 +222,7 @@ m.APC <- update(m.0, . ~ . - 1 + MA + I(P - p0) + MPr + MCr) drift <- rbind( ci.exp(m.APC, subset = "I\\(", alpha = alpha), ci.exp(m.Ad , subset = "I\\(", alpha = alpha) ) - rownames(drift) <- c("APC", "A-d") + rownames(drift) <- c(paste("APC (",drtyp,")",sep=""), "A-d") if (parm == "ADCP") m.APC <- update(m.0, . ~ . - 1 + MA + I(P - A - c0) + MPr + MCr) if (parm == "APC") { @@ -301,11 +317,11 @@ if( !ref.p & parm %in% c("APC","ADPC") ) cat( "No reference period given:\n", "Reference period for age-effects is chosen as\n", - "the median date of event: ", p0 ) + "the median date of event: ", p0, ".\n" ) if( !ref.c & parm %in% c("ACP","ADCP") ) cat( "No reference period given:\n", "Reference period for age-effects is chosen as\n", - "the median date of birth for persons with event: ", c0 ) + "the median date of birth for persons with event: ", c0, ".\n" ) class(res) <- "apc" invisible(res) } diff -Nru r-cran-epi-2.0/R/apc.LCa.R r-cran-epi-2.7/R/apc.LCa.R --- r-cran-epi-2.0/R/apc.LCa.R 1970-01-01 00:00:00.000000000 +0000 +++ r-cran-epi-2.7/R/apc.LCa.R 2016-05-25 07:19:03.000000000 +0000 @@ -0,0 +1,53 @@ +################################################################################### +### apc-LCa comparison +apc.LCa <- +function( data, # cohort reference for the interactions + keep.models = FALSE, ... ) +{ +models <- c("APa","ACa","APaC","APCa","APaCa" ) +LCa.list <- list() +length( LCa.list ) <- 5 +names( LCa.list ) <- models +for( mod in models ){ cat( mod, ":\n" ) + LCa.list[[mod]] <- LCa.fit( data = data, + model = mod, ... ) } +APC <- apc.fit( data, npar = list( A=LCa.list[[2]]$knots$a.kn, + P=LCa.list[[1]]$knots$p.kn, + C=LCa.list[[1]]$knots$c.kn ) ) +dev <- c( APC$Anova[c(2,5,3,4),2], + sapply( LCa.list, function(x) x$deviance ) ) +df <- c( APC$Anova[c(2,5,3,4),1], + sapply( LCa.list, function(x) x$df ) ) +names(dev)[1:4] <- names(df)[1:4] <- +gsub( "rift","", gsub("eriod","", gsub("ohort","", gsub("-","", +gsub( "ge", "", rownames(APC$Anova)[c(2,5,3,4)]))))) +if( keep.models ) return( list( dev = cbind( dev, df ), + apc = APC, + LCa = LCa.list ) ) + else return( cbind( dev, df ) ) +} + +show.apc.LCa <- +function( x, dev.scale=TRUE, top="Ad", ... ) +{ +if( is.list(x) ) x <- x[[1]] +TM <- matrix(NA,9,9) +rownames( TM ) <- colnames( TM ) <- +paste( rownames(x), "\n", formatC(x[,1],format="f",digits=1) ) +TM[1,2:3] <- +TM[2,c(4,5)] <- +TM[3,c(4,6)] <- +TM[4,c(7,8)] <- +TM[5,7] <- +TM[6,8] <- +TM[c(7,8),9] <- 1 +TM +bp <- list( x=c(50,30,70,50,10,90,30,70,50), + y=c(90,70,70,50,50,50,30,30,10) ) +if( !dev.scale ) + boxes.matrix( TM, boxpos=bp, hm=1.5, wm=1.5, ... ) +else { + bp$y <- 5+90*(pmin(x[top,1],x[,1])-x[9,1])/(x[top,1]-x[9,1]) + boxes.matrix( TM, boxpos=bp, hm=1.5, wm=1.5, ... ) + } +} diff -Nru r-cran-epi-2.0/R/boxes.MS.R r-cran-epi-2.7/R/boxes.MS.R --- r-cran-epi-2.0/R/boxes.MS.R 2015-05-27 15:26:55.000000000 +0000 +++ r-cran-epi-2.7/R/boxes.MS.R 2016-06-24 06:13:59.000000000 +0000 @@ -154,8 +154,8 @@ if( show.BE ) { # Derive the persons at start and at end of study in each state - Beg <- table( status( obj, at="entry", by.id=TRUE ) ) - End <- table( status( obj, at="exit" , by.id=TRUE ) ) + Beg <- as.vector( table( status( obj, at="entry", by.id=TRUE ) ) ) + End <- as.vector( table( status( obj, at="exit" , by.id=TRUE ) ) ) BE.line <- paste( ifelse( noz & Beg==0, rep(" ",nchar(BE.sep[1])+2), paste(BE.sep[1], diff -Nru r-cran-epi-2.0/R/cal.yr.R r-cran-epi-2.7/R/cal.yr.R --- r-cran-epi-2.0/R/cal.yr.R 2015-05-27 15:26:55.000000000 +0000 +++ r-cran-epi-2.7/R/cal.yr.R 2016-02-26 11:35:35.000000000 +0000 @@ -8,16 +8,12 @@ if( inherits( x, "data.frame" ) & is.null(wh) & missing(format) ) { # Indicator of where a date-type variable is - wh <- sapply( x, inherits, cl.typ ) - # The positions - wh <- (1:length(wh))[wh] + wh <- which( sapply( x, inherits, cl.typ ) ) } if( inherits( x, "data.frame" ) & is.null(wh) & !missing(format) ) { # Indicator of where the character variables are - wh <- sapply( x, is.character ) - # The positions - wh <- (1:length(wh))[wh] + wh <- which( sapply( x, is.character ) ) } if( inherits( x, "data.frame" ) & is.vector(wh) ) { diff -Nru r-cran-epi-2.0/R/clear.R r-cran-epi-2.7/R/clear.R --- r-cran-epi-2.0/R/clear.R 2015-05-27 15:26:55.000000000 +0000 +++ r-cran-epi-2.7/R/clear.R 2016-10-01 08:51:58.000000000 +0000 @@ -1,17 +1,17 @@ -# Clear the workspace of all objects whose names don't start with a full stop +# Clear the workspace of all objects whose names don't start with a +# full stop, and remove objects attached immediately after the global clear <- -function() +function () { -env <- as.environment(1) -to.go <- ls(env, all.names=FALSE) -continue <- TRUE -while (continue) { - nxt <- search()[[2]] - # bit of a botch - if (substr(nxt, 1, 8)!="package:") - detach() - else - continue <- FALSE - } -remove(list=to.go, envir=env) + env <- as.environment(1) + to.go <- ls(env, all.names = FALSE) + continue <- TRUE + while (continue) { + nxt <- search()[[2]] + if (substr(nxt, 1, 8) != "package:") + detach() + else continue <- FALSE + } + remove(list = to.go, envir = env) } + diff -Nru r-cran-epi-2.0/R/Cplot.R r-cran-epi-2.7/R/Cplot.R --- r-cran-epi-2.0/R/Cplot.R 2015-05-27 15:26:55.000000000 +0000 +++ r-cran-epi-2.7/R/Cplot.R 2016-04-20 11:25:57.000000000 +0000 @@ -25,19 +25,29 @@ ... ) { -# First convert the age-period table to an age-cohort table +# Convert the age-period table to an age-cohort table rt <- as.table( rates ) dimnames( rt ) <- list( age = age, per = per ) rtf <- data.frame( rt ) rtf$age <- as.numeric( as.character( rtf$age ) ) rtf$per <- as.numeric( as.character( rtf$per ) ) + +# Check that age and period classe have equal lengths: +wa <- diff( ag <- sort( unique(rtf$age) ) ) +wp <- diff( pg <- sort( unique(rtf$per) ) ) +if( unique(wa) != unique(wp) ) + stop("Age and period intervals must have same width:\n", + "Ages:", ag, "\n", "Periods:", pg, "\n", + "No plot with cohort produced.\n" ) + +# Table by age and cohort ac <- tapply( rtf$Freq, list( rtf$age, rtf$per-rtf$age ), mean ) coh <- as.numeric( colnames( ac ) ) if( is.null( c.lim ) ) c.lim <- range( coh, na.rm=TRUE ) + c(0,diff( range( coh ) )/30) * ann # Plot the frame -if( ann ) c.lim <- c.lim - c(diff( range( coh ) ) * xannx,0) +if( ann ) c.lim <- c.lim - c(diff( range( coh ) ) * xannx,0) matplot( coh, t(ac), type="n", xlim=c.lim, ylim=ylim, xlab=c.lab, ylab=ylab, log=log.ax, las=las, yaxt=if( !is.null( at ) ) "n" else "s" ) diff -Nru r-cran-epi-2.0/R/erl.R r-cran-epi-2.7/R/erl.R --- r-cran-epi-2.0/R/erl.R 1970-01-01 00:00:00.000000000 +0000 +++ r-cran-epi-2.7/R/erl.R 2016-10-01 07:45:53.000000000 +0000 @@ -0,0 +1,193 @@ +erl <- +function( int, + muW, + muD, + lam = NULL, + age.in = 0, + A = NULL, + immune = is.null(lam), + yll = TRUE, + note = TRUE ) +{ +# Computes expected residual life time for Well and Dis states +# respectively in an illness-death model, optionally ignoring +# the well->ill transition + +# Utility to integrate a survival function from the last point where +# it is 1, assuming points are 1 apart +trsum <- +function( x ) +{ +x[c(diff(x)==0,TRUE)] <- NA +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" ) + +# 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 ) +else sW <- surv2( int=int, muW, muD, lam, age.in = age.in, A = A ) + +# Area under the survival functions +erl <- cbind( apply( sW[,-1], 2, trsum ), + apply( sD[,-1], 2, trsum ) ) * int +colnames( erl ) <- c("Well","Dis") +rownames( erl ) <- colnames( sW )[-1] + +# Should we compute years of life lost? +if( yll ) erl <- cbind( erl, YLL = erl[,"Well"] - erl[,"Dis"] ) + +# Cautionary note +if( immune ) + { + attr( erl, "NOTE" ) <- "Calculations assume that Well persons cannot get Ill (quite silly!)." + if( note ) cat("NOTE:", attr( erl, "NOTE" ), "\n" ) + } +return( erl ) +} + +# yll is just a simple wrapper for erl, only selecting the YLL column. +yll <- +function( int, + muW, + muD, + lam = NULL, + age.in = 0, + A = NULL, + immune = is.null(lam), + note = TRUE ) erl( int = int, + muW = muW, + muD = muD, + lam = lam, + age.in = age.in, + A = A, + immune = immune, + yll = TRUE, + note = note )[,"YLL"] + +surv1 <- +function( int, mu, age.in=0, A=NULL ) +{ +# Computes the survival function from age A till the end, assuming +# that mu is a vector of mortalities in intervals of length int. +# int and mu should be in compatible units that is T and T^-1 for +# some unit T (months, years, ...) + +# age-class boundaries +age <- 0:length(mu)*int + age.in + +# cumulative rates and survival at the boundaries +Mu <- c( 0, cumsum( mu )*int ) +Sv <- exp( -Mu ) +surv <- data.frame( age=age, surv=Sv ) + +# if a vector of conditioning ages A is given +if( cond <- !is.null(A) ) + { + j <- 0 + # actual conditioning ages + cage <- NULL + for( ia in A ) + { + j <- j+1 + # Where is the age we condition on + cA <- which( diff(age>ia)==1 ) + surv <- cbind( surv, pmin( 1, surv$surv/(surv$surv[cA]) ) ) + cage[j] <- surv$age[cA] + } + } +names( surv )[-1] <- paste( "A", c( age.in, if( cond ) cage else NULL ), sep="" ) +rownames( surv ) <- NULL +return( surv ) +} + +erl1 <- +function( int, mu, age.in = 0 ) +{ +# Computes expected residual life time at all ages +age <- 0:length(mu)*int + age.in + +# Small utility: cumulative cumulative sum from the end of a vector +musmuc <- function( x ) rev( cumsum( rev(x) ) ) + +# The survival function with a 0 at end, and the integral from the upper end +surv <- surv1( int = int, mu = mu, age.in = age.in )[,2] +cbind( age = age, + surv = surv, + erl = c( musmuc( ( surv[-1]-diff(surv)/2 ) ) / + surv[-length(surv)], 0 ) * int ) +} + +surv2 <- +function( int, muW, muD, lam, age.in=0, A=NULL ) +{ +# check the vectors +if( length(muW) != length(muD) | + length(muD) != length(lam) ) + stop( "Vectors with rates must have same length:\n", + "length(muW)=", length(muW), + ", length(muD)=", length(muD), + ", length(lam)=", length(lam) ) + +# 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 +# state is muD, and that all refer to constant rates intervals of +# length int starting from age.in, conditional on survival to A +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 <- 0:length(muW)*int + age.in + +# cumulative rates at the boundaries, given survival to A +MuW <- cumsum( c( 0, muW ) * ( age > A ) ) * int +MuD <- cumsum( c( 0, muD ) * ( age > A ) ) * int +Lam <- cumsum( c( 0, lam ) * ( age > A ) ) * int + +# probability of being well +pW <- exp( -( Lam + MuW ) ) + +# probability of diagnosis at s --- first term in the integral for +# P(DM at a). Note that we explicitly add a 0 at the start so we get a +# probability of 0 of transition at the first age point +Dis <- c(0,lam) * ( age > A ) * exp( -(Lam+MuW) ) * int + +# for each age (age[ia]) we compute the integral over the range +# [0,age] of the product of the probability of diagnosis and the +# probability of surviving from diagnosis till age ia +pDM <- Dis * 0 +for( ia in 1:length(age) ) + pDM[ia] <- sum( Dis[1:ia] * exp( -(MuD[ia]-MuD[1:ia]) ) ) + # 1st term as function of s (1:ia) + # 2nd term integral over range s:age + # upper integration limit is age (ia) and the lower + # limit is the intermediate age (at DM) (1:ia) +# Finally, we add the probabilities of being in Well resp. DM to get +# the overall survival: +surv <- data.frame( age = age, surv = pDM + pW ) +return( surv ) +} + +# survival from start +surv <- wsurv2( int, muW, muD, lam, age.in=age.in, A=0 ) + +# add columns for conditioning ages +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] ) + } + } +Al <- A +for( i in 1:length(A) ) Al[i] <- max( surv$age[surv$age <= A[i]] ) +colnames( surv )[-1] <- paste( "A", c( age.in, Al ), sep="" ) + +# done! +return( surv ) +} diff -Nru r-cran-epi-2.0/R/LCa.fit.R r-cran-epi-2.7/R/LCa.fit.R --- r-cran-epi-2.0/R/LCa.fit.R 1970-01-01 00:00:00.000000000 +0000 +++ r-cran-epi-2.7/R/LCa.fit.R 2016-05-24 22:18:00.000000000 +0000 @@ -0,0 +1,547 @@ +###################################################################### +### estimation method +LCa.fit <- +function( data, A, P, D, Y, + model = "APa", # or one of "ACa", "APaC", "APCa" or "APaCa" + 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 + c.ref, # cohort reference for the interactions + npar = c(a = 6, # no. knots for main age-effect + p = 6, # no. knots for peroid-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? + alpha = 0.05, # 1 minus confidence level + eps = 1e-6, # convergence criterion + maxit = 100, # max. no iterations + quiet = TRUE ) # cut the crap +{ +# "model" must have values in c(APa/ACa/APaC/APCa/APaCa)? +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 + 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))) +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 +if( !missing(data) ) + { + if (length(match(c("A", "P", "D", "Y"), names(data))) != 4) + stop("Data frame ", deparse(substitute(data)), + " has columns:\n", names(data), + "\nmust have variables:\n", "A (age), P (period), D (cases) and Y (person-time)") + data <- data[,c("A","P","D","Y")] + data <- data[complete.cases(data),] + A <- data$A + P <- data$P + 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) + if (any(nm)) + stop("Variable", if (sum(nm) > 1) + "s", paste(c(" A", " P", " D", " Y")[nm], collapse = ","), + " missing from input") + # and that they have the same length + if( diff(range( lv <- c( length(A), + length(P), + length(D), + length(Y) ) )) != 0 ) + stop( "\nLengths of variables (", paste(paste(names(lv), + lv, sep = ":"), collapse = ", "), ") are not the same." ) + } # end of data acquisition + +# code-simplifier for knot calculation +eqqnt <- function(n) round( (1:n-0.5)/n, 2 ) +# Define knots - we compute also the knots not needed +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 ) ) + 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"] ) ) + } + +# 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(ci.ref) ) ci.ref <- median( rep( A,D) ) + +############################################################################ +# Here starts the actual modelling +commence <- Sys.time() + +# Matrices to extract the age-interaction terms at reference points +Mp <- Ns( rep(pi.ref,length(A)), knots=pi.kn, intercept=TRUE ) +Mc <- Ns( rep(ci.ref,length(A)), knots=ci.kn, intercept=TRUE ) + +# Current age-effects (in the iteration these will be term predictions) +ba <- cbind( rep(1,length(A)), 1 ) +# set to 0 if term is not in model at all +if( !mainP ) ba[,1] <- 0 +if( !mainC ) ba[,2] <- 0 + +# Main effects model with (at least one) age-interactions +mat <- glm( D ~ -1 + Ns( A, knots=a.kn, intercept=TRUE ) + + Ns( P , knots=p.kn, ref=p.ref):ba[,1] + + Ns( P-A, knots=c.kn, ref=c.ref):ba[,2], + offset = log(Y), + family = poisson ) +oldmb <- oldmat <- mat$deviance + +# Terms prediction --- three terms here. +# No need to divide by the ba, it is eiter 1 or 0 +pat <- predict( mat, type="terms" ) + +# iteration counter and continuation indicator +nit <- 0 +one.more <- TRUE + +# For simple formatting of the iteration output +fC <- function(x,d) formatC(x,format="f",digits=d) + +# now, iterate till convergence +while( one.more ) + { +nit <- nit+1 + +# Terms with main effects should be either in interaction or offset, +# so one of these should always be 0 +Pint <- Poff <- pat[,2] +Cint <- Coff <- pat[,3] +if( intP ) Poff <- Poff*0 else Pint <- Pint*0 +if( intC ) Coff <- Coff*0 else Cint <- Cint*0 +# Iteration of the age-interaction +mb <- glm( D ~ -1 + Ns( A, knots=pi.kn, intercept=TRUE ):Pint + + Ns( A, knots=ci.kn, intercept=TRUE ):Cint, + offset = pat[,1] + Poff + Coff + log(Y), + family = poisson ) + +# Get the age-interaction terms only, and if one is not needed set to 0 +ba <- predict( mb, type="terms" ) / cbind(Pint,Cint) / + cbind( ci.lin( mb, subset="pi.kn", ctr.mat=Mp)[,1], + ci.lin( mb, subset="ci.kn", ctr.mat=Mc)[,1] ) +ba[is.na(ba)] <- 0 + +# If no interaction only main should be fitted; if no main effect, set to 0 +if( !intP ) ba[,1] <- rep(1,length(A)) * mainP +if( !intC ) ba[,2] <- rep(1,length(A)) * mainC +# apc model with assumed known interactions with age +mat <- glm( D ~ -1 + Ns( A, knots=a.kn, intercept=TRUE ) + + Ns( P , knots=p.kn, ref=p.ref):ba[,1] + + Ns( P-A, knots=c.kn, ref=c.ref):ba[,2], + offset = log(Y), + family = poisson ) + +# extract age and period terms +pat <- predict( mat, type="terms" ) / cbind( 1, ba ) +pat[is.na(pat)] <- 0 + +# convergence? Check bot that the two models give the same deviance +# and that the chnage in each is small +newmat <- mat$deviance +newmb <- mb$deviance +conv <- ( reldif <- max( (abs(newmat-newmb)/ + (newmat+newmb)*2), + (oldmat-newmat)/newmat, + (oldmb -newmb )/newmb ) ) < eps +one.more <- ( !conv & ( nit < maxit ) ) +oldmat <- newmat +oldmb <- newmb +if( !quiet & nit==1 ) + cat( " Deviances: model(AT) model(A) Rel. diff.\n" ) +if( !quiet ) cat( "Iteration", formatC( nit, width=3, flag=" "), "", + fC(mat$deviance,3), + fC( mb$deviance,3), + fC( reldif, 7 ), "\n" ) + } # end of iteration loop + +# Deviance and d.f - there is a "+1" because the intercept is in both models +# but not explicit, (both models fitted with "-1"), hence the df.null +# is the total no. observations +dev <- mb$deviance +df <- mat$df.null - ( mb$df.null- mb$df.res # no. parms in mb + + mat$df.null-mat$df.res # no. parms in mat + - 1 ) # common intercept +if( conv ) cat( "LCa.fit convergence in ", nit, + " iterations, deviance:", dev, "on", df, "d.f.\n") +if( !conv ) cat( "LCa.fit *not* converged in ", nit, + " iterations:\ndeviance (AT):", mat$deviance, + ", deviance (B):" , mb$deviance, "\n", + if( VC ) "...no variance-covariance computed.\n" ) +fin <- Sys.time() +if( !quiet ) cat("...using", round(difftime(fin,commence,units="secs"),1), "seconds.\n") + +# unique values of A, P and C in the dataset for reporting effects +a.pt <- sort(unique( A)) +p.pt <- sort(unique(P )) +c.pt <- sort(unique(P-A)) + +# extract effects from final models after convergence +ax <- ci.exp( mat, subset= "a.kn", ctr.mat=Ns(a.pt,knots= a.kn,intercept=TRUE ) ) +kp <- ci.exp( mat, subset= "p.kn", ctr.mat=Ns(p.pt,knots= p.kn,ref=p.ref) ) +kc <- ci.exp( mat, subset= "c.kn", ctr.mat=Ns(c.pt,knots= c.kn,ref=c.ref) ) +pi <- ci.exp( mb , subset="pi.kn", ctr.mat=Ns(a.pt,knots=pi.kn,intercept=TRUE), Exp=FALSE ) +ci <- ci.exp( mb , subset="ci.kn", ctr.mat=Ns(a.pt,knots=ci.kn,intercept=TRUE), Exp=FALSE ) + +# Label the estimated effects +rownames( ax ) <- +rownames( pi ) <- +rownames( ci ) <- a.pt +rownames( kp ) <- p.pt +rownames( kc ) <- c.pt + +# do we bother about the correct variance-covariance? +if( VC & conv ) # ...certainly not without convergence +{ +commence <- Sys.time() +if( !quiet ) cat("...computing Hessian by numerical differentiation...\n") + +# the number of parameters for each of the 5 effects +na <- length( grep( "a.kn", names(coef(mat)) ) ) +np <- length( grep( "p.kn", names(coef(mat)) ) ) +nc <- length( grep( "c.kn", names(coef(mat)) ) ) +npi <- length( grep( "pi.kn", names(coef(mb )) ) ) +nci <- length( grep( "ci.kn", names(coef(mb )) ) ) + +# get only the parameters for effects that are non-zero (the others +# are in the models but they are 0) +ml.cf <- c( coef(mat)[c(rep(TRUE,na), + rep(mainP,np), + rep(mainC,nc))], + coef(mb)[c(rep(intP,npi), + rep(intC,nci))] ) +# and some more snappy names for the parameters: first all names +all.nam <- c( paste("ax",1:na,sep=""), + paste("kp",1:np,sep=""), + paste("kc",1:nc,sep=""), + paste("pi",1:npi,sep=""), + paste("ci",1:nci,sep="") ) +# ...then those actually present in the model +names( ml.cf ) <- all.nam[c(rep( TRUE,na), + rep(mainP,np), + rep(mainC,nc), + rep(intP,npi), + rep(intC,nci))] + +# We need the variance-covariance of the estimates as the 2nd +# derivative of the log-likelihood, D*log(lambda) - lambda*Y, +# or for eta=log(lambda), D*eta - exp(eta)*Y, +# assuming the sequence of parameters is ax, kp, kc, pi, ci +# (first A, P, C from model mat, then Pa, Ca from model mb) +# Note that we cannot simplify this calculation because the model is +# non-linear in pi,kp resp. ci,kc + +# Matrices to use in calculation of the terms of the model for each parms +MA <- Ns( A, knots= a.kn, intercept=TRUE ) +Mp <- Ns( P , knots= p.kn, ref=p.ref ) +Mc <- Ns( P-A, knots= c.kn, ref=c.ref ) +Mpi <- Ns( A, knots=pi.kn, intercept=TRUE ) +Mci <- Ns( A, knots=ci.kn, intercept=TRUE ) + +# Computing the log-likelihood for any set of parameters +llik <- +function( parms ) +{ + ax <- MA %*% parms[ 1:na] ; nn <- na +if( mainP ) { kp <- Mp %*% parms[nn+1:np] ; nn <- nn+np } else kp = rep(0,length(ax)) +if( mainC ) { kc <- Mc %*% parms[nn+1:nc] ; nn <- nn+nc } else kc = rep(0,length(ax)) +if( intP ) { pi <- Mpi %*% parms[nn+1:npi] ; nn <- nn+npi} else pi = rep(1,length(ax)) +if( intC ) { ci <- Mci %*% parms[nn+1:nci] } else ci = rep(1,length(ax)) +eta <- ax + pi*kp + ci*kc +sum( D*eta - exp(eta)*Y ) +} + +# Numerical calculation of the Hessian for llik +ivar <- -numDeriv::hessian( llik, ml.cf ) + +# Sometimes not quite positive definite, fix that after inverting the Hessian +vcov <- Matrix::nearPD( solve( ivar ) ) +vcov <- as.matrix( vcov$mat ) +fin <- Sys.time() +if( !quiet ) cat("...done - in", round(difftime(fin,commence,units="secs"),1), "seconds.\n") + +# Since we now have the variances of the parameters for each of the +# effects we can compute corrected c.i.s for the effects +se.eff <- +function( sub, cM, Alpha=alpha ) +{ +wh <- grep( sub, names(ml.cf) ) +res <- cbind( cM %*% ml.cf[wh], + sqrt( diag( cM %*% vcov[wh,wh] %*% t(cM) ) ) ) %*% ci.mat(alpha=Alpha) +colnames(res)[1] <- paste( "Joint", colnames(res)[1] ) +res +} + +# Append the corrected c.i.s to the effect objects + ax <- cbind( ax, exp( se.eff( "ax", Ns(a.pt,knots= a.kn,intercept=TRUE) ) ) ) +if( mainP ) kp <- cbind( kp, exp( se.eff( "kp", Ns(p.pt,knots= p.kn,ref=p.ref) ) ) ) +if( mainC ) kc <- cbind( kc, exp( se.eff( "kc", Ns(c.pt,knots= c.kn,ref=c.ref) ) ) ) +if( intP ) pi <- cbind( pi, se.eff( "pi", Ns(a.pt,knots=pi.kn,intercept=TRUE) ) ) +if( intC ) ci <- cbind( ci, se.eff( "ci", Ns(a.pt,knots=ci.kn,intercept=TRUE) ) ) +} + +# Collect refs and knots in lists +klist <- list( a.kn=a.kn, pi.kn=pi.kn, p.kn=p.kn, ci.kn=ci.kn, c.kn=c.kn ) +rlist <- list( pi.ref=pi.ref, p.ref=p.ref, ci.ref=ci.ref, c.ref=c.ref ) + +# Finally output object +res <- list( model = model, + ax = ax, + pi = if( intP ) pi else NULL, + kp = if( mainP ) kp else NULL, + ci = if( intC ) ci else NULL, + kc = if( mainC ) kc else NULL, + mod.at = mat, + mod.b = mb, + coef = if( VC & conv ) ml.cf else NULL, + vcov = if( VC & conv ) vcov else NULL, + knots = klist, + refs = rlist, + deviance = dev, + df.residual = df, + iter = nit ) +# Remove redundant stuff before returning +res <- res[!sapply(res,is.null)] +class( res ) <- "LCa" +invisible( res ) +} + +###################################################################### +### utility to determine the types of terms in a model +model.terms <- +function( x ) list( intP = as.logical(length(grep("Pa",x$model))), + mainP = as.logical(length(grep("P" ,x$model))), + intC = as.logical(length(grep("Ca",x$model))), + mainC = as.logical(length(grep("C" ,x$model))) ) + +###################################################################### +### print method +print.LCa <- +function( x, ... ) +{ +# terms in the model +mt <- model.terms( x ) + +# no. of parameters in each term +npar <- sapply(x$knots,length)-1 +npar[1] <- npar[1] + 1 +npar <- npar[c(TRUE,mt$intP,mt$mainP,mt$intC,mt$mainC)] +npar <- paste( paste( npar, c(rep(", ",length(npar)-2)," and ",""), sep=""), collapse=rep("") ) +cat( paste( x$model, ": Lee-Carter model with natural splines:\n", + " log(Rate) = ax(Age)", + if( mt$mainP ) " + ", + if( mt$intP ) "pi(Age)", + if( mt$mainP ) "kp(Per)", + if( mt$mainC ) " + ", + if( mt$intC ) "ci(Age)", + if( mt$mainC ) "kc(Coh)", "\nwith ", + npar, " parameters respectively.\n", + "Deviance: ", round(x$deviance,3), " on ", x$df, " d.f.\n", + sep=""), + if( !("vcov" %in% names(x)) ) "(only conditional c.i.s for effects)\n" ) +} + +###################################################################### +### summary method +summary.LCa <- +function( object, show.est=FALSE, ... ) +{ +# terms in the model +mt <- model.terms( object ) + +print( object ) + +# the number of parameters for each of the 5 effects +na <- length( grep( "a.kn", names(coef(object$mod.at)) ) ) +np <- length( grep( "p.kn", names(coef(object$mod.at)) ) ) +nc <- length( grep( "c.kn", names(coef(object$mod.at)) ) ) +npi <- length( grep( "pi.kn", names(coef(object$mod.b )) ) ) +nci <- length( grep( "ci.kn", names(coef(object$mod.b )) ) ) + +# Show knots used +cat( "\nKnots used:\n") +for( nm in names(object$knots[c(TRUE,mt$mainP,mt$intP,mt$mainC,mt$intC)]) ) + { cat( nm,"\n" ) ; print(object$knots[[nm]] ) } +cf <- rbind( ci.lin(object$mod.at)[c(rep(TRUE ,na), + rep(mt$mainP,np), + rep(mt$mainC,nc)),1:2], + ci.lin(object$mod.b )[c(rep(mt$intP ,npi), + rep(mt$intC ,nci)),1:2] ) +if( "vcov" %in% names(object) ) + { + cf <- cbind( cf, sqrt( diag(object$vcov) ) ) + colnames( cf )[-1] <- c("cond.se","joint.se") + } +if( show.est ) print( cf ) +invisible( cf ) +} + +###################################################################### +### plot method +plot.LCa <- +function( x, ... ) +{ +# terms in the model +mt <- model.terms( x ) + +# A small plot utility to exploit the structure of the effects +plc <- +function( x, ... ) matplot( as.numeric( rownames(x) ), x[,ncol(x)-2:0], + type="l", lty=1, lwd=c(3,1,1), ... ) + +plc( x$ax, col="black", xlab="Age", ylab="Age-specific rates", log="y" ) +rug( x$knots$a.kn, lwd=2 ) + +if( mt$intP ){ +plc( x$pi, col="black", xlab="Age", ylab="Relative period log-effect multiplier" ) +abline(h=1,v=x$refs$pi.ref) +rug( x$knots$pi, lwd=2 ) +} + +if( mt$mainP ){ +plc( x$kp, col="black", log="y", xlab="Date of follow-up (period)", ylab="Period effect (RR)" ) +abline(h=1,v=x$refs$p.ref) +rug( x$knots$kp, lwd=2 ) +} + +if( mt$intC ){ +plc( x$ci, col="black", xlab="Age", ylab="Relative cohort log-effect multiplier" ) +abline(h=1,v=x$refs$ci.ref) +rug( x$knots$ci, lwd=2 ) +} + +if( mt$mainC ){ +plc( x$kc, col="black", log="y", xlab="Date of birth (cohort)", ylab="Cohort effect (RR)" ) +abline(h=1,v=x$refs$c.ref) +rug( x$knots$kc, lwd=2 ) +} + +} + +###################################################################### +### predict method +predict.LCa <- +function( object, + newdata, + alpha = 0.05, + level = 1-alpha, + sim = ( "vcov" %in% names(object) ), + ... ) +{ +# What main effects and interactions are in the model +mt <- model.terms( object ) + +# is person-years suppied, otherwise use units as in the model +if( "Y" %in% names(newdata) ) Y <- newdata$Y else + Y <- rep(1,nrow(newdata)) + +# Matrices to extract effects at newdata rows + Ma <- Ns( newdata$A, knots = object$knots$a.kn, intercept = TRUE) + Mp <- Ns( newdata$P , knots = object$knots$p.kn, ref=object$refs$p.ref ) + Mc <- Ns( newdata$P-newdata$A, knots = object$knots$c.kn, ref=object$refs$c.ref ) +Mpi <- Ns( newdata$A, knots = object$knots$pi.kn, intercept = TRUE) +Mci <- Ns( newdata$A, knots = object$knots$ci.kn, intercept = TRUE) + +# Default terms values for models without interactions +kp <- kc <- rep( 0, nrow(newdata) ) +pi <- ci <- rep( 1, nrow(newdata) ) + +# P, C and interaction term(s) if included in the model +if( mt$intP ) { + pi <- ci.lin( object$mod.b , subset="pi.kn", ctr.mat=Mpi )[,1] + kp <- ci.lin( object$mod.at, subset= "p.kn", ctr.mat=Mp )[,1] + } +if( mt$intC ) { + ci <- ci.lin( object$mod.b , subset="ci.kn", ctr.mat=Mci )[,1] + kc <- ci.lin( object$mod.at, subset= "c.kn", ctr.mat=Mc )[,1] + } + +# 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. +pr0 <- ci.exp( object$mod.at, alpha=alpha, ctr.mat=cbind(Ma,Mp*pi,Mc*ci) ) + +# Then fitted values from mod.b +# But mod.b has an offset beyond log(Y), namely all the APC terms +lp.b <- ci.lin( object$mod.b , ctr.mat=cbind(Mpi*kp,Mci*kc) )[,1:2] +lp.b[,1] <- lp.b[,1] + ci.lin( object$mod.at, ctr.mat=cbind(Ma,Mp*(!mt$intP),Mc*(!mt$intC)) )[,1] +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... +if( sim ) # also renders TRUE if sim is numerical (and not 0) + { +if( is.logical(sim) & sim ) sim <- 1000 +# Check that there is a vcov component of the model +if( !( "vcov" %in% names(object) ) ) + warning( + "No variance-covariance in LCa object, only conditional c.i.s available.\n", + "no simulation (prametric bootstrap) is done.\n" ) +else { +# require( MASS ) +# using the parametric bootstrap based on the parameters and the +# (numerically computed) Hessian +eta <- NArray( list( pt = 1:nrow(pr0), + it = 1:sim ) ) +parms <- MASS::mvrnorm( n = sim, + mu = object$coef, + Sigma = object$vcov ) +na <- ncol( Ma ) +npi <- ncol( Mpi ) + np <- ncol( Mp ) +nci <- ncol( Mci ) + nc <- ncol( Mc ) +# Compute the linear predictor in each of the simulated samples +# period and cohort effects if not in the model +kp <- kc <- rep( 0, nrow(newdata) ) +pi <- ci <- rep( 1, nrow(newdata) ) +for( i in 1:sim ){ + ax <- Ma %*% parms[i, 1:na] ; nn <- na +if( mt$mainP ) { kp <- Mp %*% parms[i,nn+1:np] ; nn <- nn+np } +if( mt$mainC ) { kc <- Mc %*% parms[i,nn+1:nc] ; nn <- nn+nc } +if( mt$intP ) { pi <- Mpi %*% parms[i,nn+1:npi] ; nn <- nn+npi} +if( mt$intC ) { ci <- Mci %*% parms[i,nn+1:nci] } +eta[,i] <- ax + kp*pi + kc*ci + } +# predicted rates with bootstrap confidence limits +pr.sim <- exp( t( apply( eta, 1, quantile, + probs=c(0.5,alpha/2,1-alpha/2), + na.rm=TRUE ) ) ) +colnames( pr.sim )[1] <- "Joint est." +return( pr.sim ) + } +} +else return( pr0 ) +} + + diff -Nru r-cran-epi-2.0/R/lexis.R r-cran-epi-2.7/R/lexis.R --- r-cran-epi-2.0/R/lexis.R 2016-01-04 13:51:22.000000000 +0000 +++ r-cran-epi-2.7/R/lexis.R 2016-10-03 19:58:32.000000000 +0000 @@ -519,10 +519,10 @@ ### Generic functions -### Methods for data.frame drop Lexis attributes, so we need a Lexis -### method that adds them again +### Methods for data.frame drop Lexis attributes, so we need Lexis +### methods that retain them -subset.Lexis <- function(x, ...) +subset.Lexis <- function(x, ... ) { y <- subset.data.frame(x, ...) attr(y,"breaks") <- attr(x, "breaks") @@ -532,14 +532,25 @@ } `[.Lexis` <- -function( x, ... ) +function(x, ...) { - structure( NextMethod(), - breaks = attr(x, "breaks"), - time.scales = attr(x, "time.scales"), - time.since = attr(x, "time.since") ) + y <- NextMethod(x) + if (is.data.frame(y)) { + for (a in c("class", "time.scales", "time.since", "breaks")) { + data.table::setattr(y ,a, attr(x, a))} + } + y } +# `[.Lexis` <- +# function( x, ... ) +# { +# structure( NextMethod(), +# breaks = attr(x, "breaks"), +# time.scales = attr(x, "time.scales"), +# time.since = attr(x, "time.since") ) +# } + merge.data.frame <- function(x, y, ...) { if (is.Lexis(x)) @@ -607,7 +618,7 @@ is.nul <- sapply( allargs, is.null ) if( !all(is.lex[!is.nul]) ) stop( "All arguments must be Lexis objects,\n", - "arguments number ", which(!is.lex & !is.nul), " are not." ) + "arguments number ", which(!is.lex & !is.null), " are not." ) # Put them all together allargs <- allargs[!is.nul] res <- plyr::rbind.fill( allargs ) diff -Nru r-cran-epi-2.0/R/Ns.r r-cran-epi-2.7/R/Ns.r --- r-cran-epi-2.0/R/Ns.r 2015-08-12 12:30:51.000000000 +0000 +++ r-cran-epi-2.7/R/Ns.r 2016-04-21 13:34:47.000000000 +0000 @@ -1,83 +1,177 @@ +# ------------------------------------------------------------- +# extension of ns() from splines, to allow for clamping. + +ns.ld <- function(x, df = NULL, knots = NULL, intercept = FALSE, + Boundary.knots = range(x), fixsl=c(FALSE,FALSE) ) +{ + nx <- names(x) + x <- as.vector(x) + nax <- is.na(x) + if(nas <- any(nax)) x <- x[!nax] + + + if(!missing(Boundary.knots)) { + Boundary.knots <- sort(Boundary.knots) + outside <- (ol <- x < Boundary.knots[1L]) | (or <- x > Boundary.knots[2L]) + } + + else outside <- FALSE + + if(!missing(df) && missing(knots)) { + ## df = number(interior knots) + 1 + intercept + nIknots <- df - 1 - intercept + if(nIknots < 0) { + nIknots <- 0 + warning("'df' was too small; have used ", 1 + intercept) + } + + knots <- if(nIknots > 0) { + knots <- seq.int(0, 1, + length.out = nIknots + 2L)[-c(1L, nIknots + 2L)] + stats::quantile(x[!outside], knots) + } ## else NULL + + } else nIknots <- length(knots) + Aknots <- sort(c(rep(Boundary.knots, 4L), knots)) + + if(any(outside)) { + basis <- array(0, c(length(x), nIknots + 4L)) + + if(any(ol)) { + k.pivot <- Boundary.knots[1L] + xl <- cbind(1, x[ol] - k.pivot) + tt <- splines::spline.des(Aknots, rep(k.pivot, 2L), 4, c(0, 1))$design + basis[ol, ] <- xl %*% tt + } + + if(any(or)) { + k.pivot <- Boundary.knots[2L] + xr <- cbind(1, x[or] - k.pivot) + tt <- splines::spline.des(Aknots, rep(k.pivot, 2L), 4, c(0, 1))$design + basis[or, ] <- xr %*% tt + } + + + if(any(inside <- !outside)) + basis[inside, ] <- splines::spline.des(Aknots, x[inside], 4)$design + } + + else basis <- splines::spline.des(Aknots, x, 4)$design + +# Fix the clamping +if( !is.logical(fixsl) ) warning( "fixsl elements must be of mode logical" ) +# Only the 4th parameter affected, should be either 1 or 2 in the two positions +const <- splines::spline.des( Aknots, Boundary.knots, 4, c(2-fixsl[1],2-fixsl[2]) )$design + + if(!intercept) { + const <- const[, -1 , drop = FALSE] + basis <- basis[, -1 , drop = FALSE] + } + + qr.const <- qr(t(const)) + basis <- as.matrix((t(qr.qty(qr.const, t(basis))))[, - (1L:2L), drop = FALSE]) + n.col <- ncol(basis) + + if(nas) { + nmat <- matrix(NA, length(nax), n.col) + nmat[!nax, ] <- basis + basis <- nmat + } + + dimnames(basis) <- list(nx, 1L:n.col) + a <- list(degree = 3, knots = if(is.null(knots)) numeric(0) else knots, + Boundary.knots = Boundary.knots, intercept = intercept) + attributes(basis) <- c(attributes(basis), a) + class(basis) <- c("cns", "basis", "matrix") + basis +} + #------------------------------------------------------------------------------- -# A wrapper for ns that automatically takes the smallest and largest knots -# as boundary knots without further ado, but also allows cenering -# around a reference and detrending by means of projection -Ns <- function( x, ref = NULL, +# An extension of ns that automatically takes the smallest and largest knots +# as boundary knots without further ado, but also allows centering +# around a reference and detrending by means of projection, as well as +# clamping. + +Ns <- function( x, ref = NULL, df = NULL, knots = NULL, intercept = FALSE, Boundary.knots = NULL, - detrend = FALSE ) + fixsl = c(FALSE,FALSE), + detrend = FALSE ) { -## Check sensibility of arguments -if( !is.null(ref) ) - { - if( !is.vector(ref) ) - stop( "Argument 'ref' must be a scalar, but it is a ", class(ref), "." ) - if( is.vector(ref) & length(ref)>1 ) - stop( "Argument 'ref' must be a scalar, but has length ", length(ref), "." ) - if( intercept ) - { - warning( "ref= specified, hence intercept=TRUE is ignored") + ## Check sensibility of arguments + if( !is.null(ref) ) { + if( !is.vector(ref) ) + stop( "Argument 'ref' must be a scalar, but it is a ", class(ref), "." ) + if( is.vector(ref) & length(ref)>1 ) + stop( "Argument 'ref' must be a scalar, but has length ", length(ref), "." ) + if( intercept ) { + warning( "ref= specified, hence intercept=TRUE is ignored") + intercept <- FALSE + } + } + ## Detrending required? + if( any(detrend>0) ) { # covers both logical and vector + if( any(detrend<0) ) + stop( "Some elements of weight are <0, e.g. no", + (ww <- which(detrend<0))[1:min(5,length(ww))], "." ) + if( !(length(detrend) %in% c(1,length(x))) ) { + warning( "Weights in inner product diagonal matrix set to 1") + weight <- rep(1,length(x)) + } + else weight <- if( is.numeric(detrend) ) detrend else rep(1,length(x)) + detrend <- TRUE + } + if( detrend & intercept ) { + warning( "detrend= specified, hence intercept=TRUE is ignored") intercept <- FALSE } - } -## Detrending required? -if( any(detrend>0) ) ## covers both logical and vector - { - if( any(detrend<0) ) - stop( "Some elements of weight are <0, e.g. no", - (ww <- which(detrend<0))[1:min(5,length(ww))], "." ) - if( !(length(detrend) %in% c(1,length(x))) ) - { - warning( "Weights in inner product diagonal matrix set to 1") - weight <- rep(1,length(x)) + if( detrend & any(!is.na(fixsl)) ) { + warning( "detrend= specified, hence fixsl argument is ignored") + fixsl=c(NA,NA) } - else weight <- if( is.numeric(detrend) ) detrend else rep(1,length(x)) - detrend <- TRUE - } -if( detrend & intercept ) - { - warning( "detrend= specified, hence intercept=TRUE is ignored") - intercept <- FALSE - } -## Here is the specification of the spline basis -## df= specified -if( !is.null(df) ) - MM <- ns( x, df = df, intercept = (intercept & is.null(ref)) ) -else -## knots= specified -{ -if( is.null( Boundary.knots ) ) + ## Here is the specification of the spline basis + ## df= specified + if( !is.null(df) ) + MM <- ns.ld( x, df = df, intercept = (intercept & is.null(ref)), fixsl = fixsl ) + else + ## knots= specified { - if( !is.null( knots ) ) + if( is.null( Boundary.knots ) ) { - knots <- sort( unique( knots ) ) - ok <- c(1,length(knots)) - Boundary.knots <- knots[ok] - knots <- knots[-ok] + if( !is.null( knots ) ) + { + knots <- sort( unique( knots ) ) + ok <- c(1,length(knots)) + Boundary.knots <- knots[ok] + knots <- knots[-ok] + } } + MM <- ns.ld( x, knots = knots, + Boundary.knots = Boundary.knots, + intercept = (intercept & is.null(ref)), + fixsl = fixsl ) } -MM <- ns( x, knots = knots, Boundary.knots = Boundary.knots, - intercept = (intercept & is.null(ref)) ) -} -## Reference point specified ? -if( !is.null(ref) ) + ## Reference point specified ? + if( !is.null(ref) ) + { + MM <- MM - ns.ld( rep(ref,length(x)), + knots = attr(MM,"knots"), + Boundary.knots = attr(MM,"Boundary.knots"), + fixsl = fixsl ) + } + + ## Detrending required ? + if( detrend ) { - MM <- MM - ns( rep(ref,length(x)), - knots = attr(MM,"knots"), - Boundary.knots = attr(MM,"Boundary.knots") ) - } -## Detrending required ? -if( detrend ) - { - DD <- detrend( MM, x, weight=weight ) - ## NOTE: detrend does not preserve attributes - for( aa in c("degree","knots","Boundary.knots","intercept","class") ) - attr( DD, aa ) <- attr( MM, aa ) - attr( DD, "detrend" ) <- TRUE - attr( DD, "proj.wt" ) <- weight - MM <- DD + DD <- detrend( MM, x, weight=weight ) + ## NOTE: detrend does not preserve attributes + for( aa in c("degree","knots","Boundary.knots","intercept","class") ) + attr( DD, aa ) <- attr( MM, aa ) + attr( DD, "detrend" ) <- TRUE + attr( DD, "proj.wt" ) <- weight + MM <- DD } -return( MM ) + return( MM ) } diff -Nru r-cran-epi-2.0/R/rm.tr.R r-cran-epi-2.7/R/rm.tr.R --- r-cran-epi-2.0/R/rm.tr.R 1970-01-01 00:00:00.000000000 +0000 +++ r-cran-epi-2.7/R/rm.tr.R 2016-07-20 13:37:36.000000000 +0000 @@ -0,0 +1,45 @@ +rm.tr <- +function( obj, from, to ) + { + # checks + if( from==to) stop( "Don't be silly, 'from' and 'to' are identical." ) + if( !(from %in% levels(obj) ) ) stop( "'from' not a state in the object." ) + if( !( to %in% levels(obj) ) ) stop( " 'to' not a state in the object." ) +### These things do not change over the purging iterations + # Sort the rows of the object and count them + obj <- obj[order(obj$lex.id,obj[,timeScales(obj)[1]]),] + nrw <- nrow( obj ) + # First obs for each person + no1 <- !duplicated( obj$lex.id ) + wh1 <- which( no1 ) +### Utility function doing the work inside the loop +purge.one <- +function( obj, from, to, nrw, wh1 ) + { + # Where are the illegal transitions + chX <- ( paste( obj$lex.Cst, obj$lex.Xst ) == paste( from, to ) ) + whX <- which( chX ) + if( length(whX)>0 ) + { + # Change lex.Xst in this record + obj$lex.Xst[whX] <- obj$lex.Cst[whX] + # and lex.Cst in next record, but only if it is not a new person + whX <- setdiff( whX[whX Follow-up data with the Epi package , and the corresponding R-code . -
  • +
  • Simulation in multistate models with multiple timescales - , and the corresponding R-code . + , and the corresponding R-code . Here is the website for the Epi package. diff -Nru r-cran-epi-2.0/vignettes/simLexis.rnw r-cran-epi-2.7/vignettes/simLexis.rnw --- r-cran-epi-2.0/vignettes/simLexis.rnw 2016-01-04 08:17:14.000000000 +0000 +++ r-cran-epi-2.7/vignettes/simLexis.rnw 2016-10-02 07:41:34.000000000 +0000 @@ -31,7 +31,6 @@ \usepackage[ae,hyper]{Rd} \usepackage[dvipsnames]{xcolor} \usepackage[super]{nth} -\usepackage[retainorgcmds]{IEEEtrantools} \usepackage{makeidx,Sweave,floatflt,amsmath,amsfonts,amsbsy,enumitem,dcolumn,needspace} \usepackage{ifthen,calc,eso-pic,everyshi} \usepackage{booktabs,longtable,rotating,graphicx}