[R-sig-Epi] splitLexis
BXC (Bendix Carstensen)
bxc at steno.dk
Mon Feb 25 21:13:00 CET 2013
Hi Dimitis,
You have a typo in the definition of the cohR Lexis object:
It should be "age = doz-dob" , that is doz instead of doe.
No need to call the timescale "age1" either, it just produces trouble when putting the tables together for analysis.
I assume that doz is the date of screening.
A more serious concern is however the setup with 5-year time bands,
You should use smaller time bands now you have the possibility. In studies where you (at least pre-1970) were limited by population statistics to 5-year bands you were restricted to the assumption that breast cancer incidence was the same for women aged 45.5 and 49.5 years. And dramatically different for those aged 50.5.
Use instead 1-year time bands and assign each cell in the table the average age / date of follow-up, and analyze rates with some smooth function of these. When you actually have individual data you should exploit it.
You can of course replace "1 year" in the above by "3 months" (or even "17 days"); even if you have ages 30-95 and a 35 year period you are still only looking at 36400 cells, quite handy for analysis.
So use breaks=seq(35,95,1/4) instead.
Here is a full account of such analyses underlying a paper in Diabetologia on cancer among diabetes patients:
http://bendixcarstensen.com/DMCa/Diabetologia/Analyses.pdf
(Data here are actually in 1-year Lexis triangles).
p. 70 is an example of the type of model used.
The paper and other stuff is in the same folder.
b.r.
Bendix
______________________________________________
Bendix Carstensen
Senior Statistician
Epidemiology
Steno Diabetes Center A/S
Niels Steensens Vej 2-4
DK-2820 Gentofte
Denmark
+45 44 43 87 38 (direct)
+45 30 75 87 38 (mobile)
bxc at steno.dk
http://BendixCarstensen.com
www.steno.dk
> -----Original Message-----
> From: r-sig-epi-bounces at r-project.org [mailto:r-sig-epi-bounces at r-
> project.org] On Behalf Of dimitris m
> Sent: 25. februar 2013 16:52
> To: r-sig-epi at r-project.org
> Subject: [R-sig-Epi] (no subject)
>
> Hello to all,
>
> I have a problem using the Epi package and I would like your help.
>
> my data set looks like this:
>
> id date of born date of entry date of first
> screen date of exit
>
> the date of entry is fixed for all obs to be a specific date. Now, I want
> to calculate person years of no screening, thus date first screen-date of
> entry and person years of screening therefore, date of exit-date of first
> screen.
>
> I have created two lexis objects, one for the non-screening period
> calculating the age at entry (date of entry-date of birth) , and years of
> no screening and one for the screening period calculating again the same
> age and person years of screening.
>
> Then I want to create an age/time matrix with person years, so I split the
> lexis object first on age and then on calendar year with splitLexis
> function.
> Finally, I tabulate in a matrix and all seems fine for the person years
> for the non exposed period.
>
> However, for the screening period I got a warning message saying that I
> have multiple time bands when I am trying to tabulate into the final
> matrix.
>
> this is how my algorithm looks like
>
>
>
> ########################################
> ####calculation of no-screening years###
> ########################################
>
>
> library(foreign)
>
> data1<-read.dta("T:\\person_no_screening_years.dta")
>
> attach(data1)
> names(data1)
>
> bbb<-data.frame(data1)
>
> # Convert the character dates into numerical variables (fractional years)#
> library(Epi)
> dobt <- cal.yr( bbb$date_born, format="%d/%m/%Y" ) doet <- cal.yr(
> bbb$date_entry, format="%d/%m/%Y" ) doxt <- cal.yr( bbb$min_scr_date ,
> format="%d/%m/%Y" )
>
>
>
> # Set up the dataset as a Lexis object # # with age and calendar time as
> time-scales #
>
> cohh <- Lexis( entry = list( per=doet,
> age=doet-dobt ),
> exit = list( per=doxt ),
> exit.status = cases,
> data = bbb )
>
> # Split time along one time-axis #
> La1 <-splitLexis( cohh, "aget",
> breaks=c(30,35,40,45,50,55,60,65,70,75,80,85,90,95))
> # Split time along the calendar time-axes #
> Lap1 <-splitLexis( La1 , "per",
> breaks=c(1975,1980,1985,1990,1995,2000,2005,2010))
>
> # Tabulate events and person-years #
> PYtab2 <-
> with( Lap1, xtabs( cbind( C=lex.Xst, Y=lex.dur ) ~
> timeBand(Lap1,"aget","left") +
> timeBand(Lap1,"per","left") ) )
> names(dimnames(PYtab2)) <- c("","Period","Aget") T<-round(ftable( PYtab2,
> row.vars=c(3,1) )) T
>
>
>
> write.table(T,"T:\\Years_no_screening.txt")
> #################################################################
> #################################################################
>
>
> ######################################
> ####calculation of screening years####
> ######################################
>
>
> library(foreign)
>
> data2<-read.dta("T:\\person_screening_years_new.dta")
>
> attach(data2)
> names(data2)
>
> vvvv<-data.frame(data2)
>
> # Convert the character dates into numerical variables (fractional years)#
> library(Epi)
> dob <- cal.yr( vvvv$date_born, format="%d/%m/%Y" ) doe <- cal.yr(
> vvvv$date_entry, format="%d/%m/%Y" ) dox <- cal.yr( vvvv$date_exit ,
> format="%d/%m/%Y" ) doz <- cal.yr( vvvv$min_scr_date, format="%d/%m/%Y" )
>
>
> # Set up the dataset as a Lexis object # # with age and calendar time as
> time-scales #
>
> cohR <- Lexis( entry = list( per=doz,
> age1=doe-dob ),
> exit = list( per=dox ),
> exit.status = cases,
> data = vvvv )
>
> # Split time along one time-axis #
> LaR <-splitLexis( cohR, "age1",
> breaks=c(30,35,40,45,50,55,60,65,70,75,80,85,90,95) )
>
>
> # Split time along the calendar time-axes # LapR <-splitLexis( LaR ,
> "per", breaks=c(1975,1980,1985,1990,1995,2000,2005,2010))
>
>
>
> # Tabulate events and person-years #
> PYtab3R <-
> with( LapR, xtabs( cbind( C=lex.Xst, Y=lex.dur ) ~
> timeBand(LapR,"age1","left") +
> timeBand(LapR,"per","left") ) )
> names(dimnames(PYtab3R)) <- c("","Period","Age") FR<-round(ftable(
> PYtab3R, row.vars=c(3,1) )) FR
>
>
>
> write.table(FR,"T:\\Years_screening.txt")
>
> tapply( status(LapR,"exit")==1, list( timeBand(LapR,"age1","left"),
> timeBand(LapR,"per","left") ), sum )
> tapply( dur(LapR), list( timeBand(LapR,"age1","left"),
> timeBand(LapR,"per","left") ), sum )
>
>
>
> Could you please help me with that???
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Epi at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-epi
More information about the R-sig-Epi
mailing list