发帖
雷达卡
Age-period-cohort analysis can alternatively be done by the package Epi.
Examples######################### Belgian lung cancer######## 1. Get apc.data.list# This is ready made. For other data construct list using apc.data.listdata.list <- data.Belgian.lung.cancer()objects(data.list)data.list######## 2. Plot data# Plot all data.# Note a warning is produced because the defaults settings# lead to an unbalanced grouping of data.apc.plot.data.all(data.list)# Or make individual plots.# Plot data sums.apc.plot.data.sums(data.list)# Plot sparsity to see where data are thin.# Plots are blank with default settings# ... therefore change sparsity.limits.apc.plot.data.sparsity(data.list)dev.new()apc.plot.data.sparsity(data.list,sparsity.limits=c(5,10))# Plot data using different pairs of the three time scales.# This plot is done for mortality ratios.# All plots appear to have approximately parallel lines.# This indicates that interpretation should be done carefully.apc.plot.data.within(data.list,"m",1)######## 3. Get a deviance table# Need to input distribution.# The table show that the sub-models "AC" and "Ad"# cannot be rejected relative to the unrestricted "APC" modelapc.fit.table(data.list,"poisson.dose.response")######## 4. Estimate selected models# Consider "APC" and "Ad"# Consider also the sub-model "A", which is not supported by# the tests in the deviance tablefit.apc <- apc.fit.model(data.list,"poisson.dose.response","APC")fit.at <- apc.fit.model(data.list,"poisson.dose.response","Ad")fit.a <- apc.fit.model(data.list,"poisson.dose.response","A")# Get coefficients for canonical parameters throughfit.apc$coefficients.canonicalfit.at$coefficients.canonical######## 5. Plot probability transforms of responses given fit# Black circle are used for central part of distribution.# Triangles are used in tails, green/blue/red as responses are further in tail# No sign of mis-specification for "APC" and "Ad": there are many# black circles and only few coloured triangles.# In comparison the model "A" yields more extreme observations.# That model is not supported by the data. # To get numerical values see apc.plot.fit.ptapc.plot.fit.pt(fit.apc)apc.plot.fit.pt(fit.at)apc.plot.fit.pt(fit.a)######## 6. Plot estimated coefficients # Consider "APC" and "Ad"# The first row of plots show double differences of paramters# The second row of plots shows level and slope determining a linear plane# The third row shows double sums of double differences,# all identified to be zero at the begining and at the end.# Thus the plots in third row must be interpreted jointly with those in the# second row. The interpretation of the third row plots# is that they show deviations from linear trends. The third row plots are# not invariant to changes to data arrayapc.plot.fit(fit.apc)dev.new()apc.plot.fit(fit.at)######## 7. Recursive analysis# Cut the first period group and redo analysisdata.list.subset.1 <- apc.data.list.subset(data.list,0,0,1,0,0,0)apc.fit.table(data.list.subset.1,"poisson.dose.response")######## 8. Effect of ad hoc identification# At first a subset is chosen where youngest age and cohort groups# are truncated. This way sparsity is eliminated# and ad hoc identification effects are dominated by estimation# uncertainty. Then consider# Plot 1: parameters estimated from data without first age groups# Plot 2: parameters estimated from all data# Note that estimates for double difference very similar.# Estimates for linear slopes are changed because the indices used# for parametrising these are changed# Estimates for detrended double sums of age and cohort double differences# are changed, because they rely on a particular ad hoc identifications# that have changed. Nonetheless these plots are useful to evaulate# variation in time trends over and above linear trends.data.list <- data.Belgian.lung.cancer()data.list.subset <- apc.data.list.subset(data.list,2,0,0,0,0,0)fit.apc <- apc.fit.model(data.list,"poisson.dose.response","APC")fit.apc.subset <- apc.fit.model(data.list.subset,"poisson.dose.response","APC")apc.plot.fit(fit.apc.subset,main.outer="1. Belgian lung cancer: cut first two age groups")dev.new()apc.plot.fit(fit.apc,main.outer="2. Belgian lung cancer data: all data")
好乱啊
京公网安备 11010802022788号
