################################################################################
# APC_workshop_prep.R
#
# Created:  2010-04-07 by Caroline D
# Updated:  2010-04-19 by Caroline D
#           2010-04-20 by BxC
#           2010-04-21 by Caroline D
#
# Uses: breast_i_6206.txt
#
# Package: Epi
#
# Software: R
#
# Aim: Exercise as preparation for the APC workshop
#
################################################################################


# Set working directory - change to your own!
# setwd("Z:/APC_workshop_2010/Data/Breast")

# Install and load the Epi-package, if not installed on your local computer
# install.packages("Epi", dependencies = TRUE)
library(Epi)


################################
# Breast cancer incidence
# for women in the Nordic
# countries, ages 25+
################################

# Read data and have a look at it
NC <- read.table( "http://staff.pubhealth.ku.dk/~bxc/APC/data/breast_i_6206.txt",
                  header = TRUE )  # Read in data
NC$P <- NC$P+0.5                # Mean date of follow-up
str( NC )                       # Display structure of data set
summary( NC )                   # Summary of the variables
head( NC )                      # Returns first part of vector

with( NC, table( A, P ) )                     # Display table (Age*Period)
with( NC, tapply( D, list(A,cou), sum ) )     # List no of cases per age and country

# Make an empty list to hold the results  (see ?list for help)
apc.list <- list()

# Loop over the countries (see Epi manual for help on apc.fit)
for( ic in levels(NC$cou) )
{
apc.list[[ic]] <-                 # Put the result into the i'th element of the list
apc.fit( subset( NC, cou==ic & A>=25 ), # Data
         dr.extr="Holford",             # How drift is extracted using Holford method
         parm="ACP",                    # What estimates to present [rates, RR, ...]
			                                  #  also which effect to put to zero slope
         npar=c(5,6,10),                # Number of parameters (age,period,cohort)
         scale=10^5,                    # Scale to display rates
         ref.c=1920 )                   # Reference cohort
}

# Create APC frame for the graph, type ?apc.frame for help
par( las=1, mar=c(4,4,1,4) )
res <- apc.frame(
    a.lab=seq(20,90,10),                     # Num vector of labels for age axis
    cp.lab=c(seq(1880,2010,10)),             # Num vector of labels for cohort-period axis
    r.lab=c(c(2,5),                          # Num vector of labels for rate axis
            c(1,2,5)*10,
            c(1,2,5)*100),
    a.tic=seq(20,90,5),                      # Location of add tick marks on age scale
    r.tic=c(3:10,2:9*10,2:5*100),            # Location of add tick marks on rate scale
    gap=13,                                  # Gap between age scale and cohort-period scale
    a.txt="Age",                             # Text for age axis
    r.txt="Rate",                            # Text for rate axis
    cp.txt="Date")                           # Text for cohort-period axis

# Create graph (see Epi manual for help on apc.lines)
clr <- c("black","red","blue","green")
for( ic in 1:length(apc.list) )
{
apc.lines(  apc.list[[ic]],
            frame.par=res,
            col=clr[ic],
            lwd=3)
text( 2005-options()[["apc.frame.par"]][1],
      500*(0.8)^ic,
      levels(NC$cou)[ic],
      col=clr[ic],
      font=2,
      adj=1,
      cex=1.2)
}