#### APC分析主要函数 #### https://blog.csdn.net/NickyCat/article/details/118636690?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522164585075016781683955469%2522%252C%2522scm%2522%253A%252220140713.130102334..%2522%257D&request_id=164585075016781683955469&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2~all~baidu_landing_v2~default-2-118636690.pc_search_result_control_group&utm_term=GBD%E6%95%B0%E6%8D%AE%E5%BA%93&spm=1018.2226.3001.4187 function_year5 <- function(table_name, start_year, end_year, current_year){ remain <- current_year - floor((current_year - start_year)/5) * 5 year_names <- NULL for (i in start_year:end_year) { if((i - current_year)/5 - floor((i - current_year)/5) == 0){ if(i == remain){ temp <- paste(start_year, i, sep = '-') year_names <- append(year_names, temp) } else{ temp <- paste(i-4, i, sep = '-') year_names <- append(year_names, temp) } } } table_name <- as.data.frame(table_name) new_years <- seq(start_year,end_year,1) new_table <- as.data.frame(matrix(data = rep(0, length(year_names)*nrow(table_name)), ncol = length(year_names), nrow = nrow(table_name))) %>% as.data.frame() colnames(new_table) <- year_names j = 1 for (i in 1:(end_year - start_year + 1)){ if((new_years[i] - current_year)/5 - floor((new_years[i] - current_year)/5) != 0){ new_table[, year_names[j]] <- new_table[,year_names[j]] + table_name[,as.character(new_years[i])] } else{ if(j == 1){ new_table[,year_names[j]] <- (new_table[,year_names[j]] + table_name[,as.character(new_years[i])]) / (remain - start_year + 1) } else{ new_table[,year_names[j]] <- (new_table[,year_names[j]] + table_name[,as.character(new_years[i])]) / 5 } j = j + 1 } } return(new_table) } made_country_data<-function(data_select,sex_name,output,output_path,measure,order_list,population_all,location_ID, start_year, end_year, current_year){ data_select2 <- data_select[data_select$measure_name == measure, ] age_names <- unique(data_select2$age_name) for (region_name in order_list){ sex_keyword=tolower(sex_name) population <- population_all %>% filter(location_name %in% region_name & age_group_name %in% age_names & sex_name %in% sex_keyword & location_id %in% location_ID) population <- population[,-1] population$age_group_id <- as.factor(population$age_group_id) population_list <- split(population, by = 'sex_name') population_n <- lapply(population_list, function(x){ reshape2::dcast(data = x[,-c(1:2)], age_group_id+ age_group_name ~ year_id, value.var = 'val') }) population_n <- lapply(population_n, function(x)arrange(x,age_group_id)) population_g <- list() for (i in 1:length(population_n)) { population_g[[i]] <- function_year5(population_n[[i]], start_year, end_year, current_year) rownames(population_g[[i]]) <- population_n[[i]]$age_group_name names(population_g)[i] <- tolower(names(population_n)[i]) } Disease_incidence <- subset(data_select2, location_name %in% region_name & sex_name %in% sex_name) Disease_incidence <- Disease_incidence[,c('sex_name','year','age_id','age_name','val')] Disease_list <- split(as.data.table(Disease_incidence), by = 'sex_name') names(Disease_list) <- tolower(names(Disease_list)) Disease_incidence_n <- lapply(Disease_list, function(x){ reshape2::dcast(data = x[,-c(1)], age_id+age_name ~ year, value.var = 'val') }) Disease_incidence_g <- list() for (i in 1:length(Disease_incidence_n)) { Disease_incidence_g[[i]] <- function_year5(Disease_incidence_n[[i]], start_year, end_year, current_year) rownames(Disease_incidence_g[[i]]) <- Disease_incidence_n[[i]]$age names(Disease_incidence_g)[i] <- tolower(names(Disease_incidence_n)[i]) } incidecen_population_list <- list() for (n in 1:length(sex_keyword)){ sex <- sex_keyword[n] for (i in 1:(2*ncol(population_g[[sex]]))){ if(i == 1){ incidecen_population <- Disease_incidence_g[[sex]][,i] %>% as.data.frame()} else{ if(i%%2==0){ incidecen_population <- cbind(incidecen_population,population_g[[sex]][,ceiling(i/2)])} else{ incidecen_population <- cbind(incidecen_population,Disease_incidence_g[[sex]][,ceiling(i/2)])} } } names(incidecen_population) <- rep(names(population_g[[sex]]),each=2) rownames(incidecen_population) <- rownames(population_g[[sex]]) incidecen_population_list[[n]] <- incidecen_population names(incidecen_population_list)[[n]] <- sex } if (output==TRUE){ for (sex_name in names(incidecen_population_list)){ data <- incidecen_population_list[sex_name] if (!file.exists(paste0(output_path,sex_name,'/',measure))){ dir.create(paste0(output_path,str_to_title(sex_name))) dir.create(paste0(output_path,str_to_title(sex_name),'/',measure)) } write.csv(data,paste0(output_path,str_to_upper(sex_name),'/',measure,'/',region_name,'-',sex_name,'-',measure,'-result','.csv')) } } return(incidecen_population_list) } } made_country_apc<-function(data,start_year,start_age,interval_year,label){ R <- prepare_rates(data, StartYear=as.numeric(start_year), StartAge=as.numeric(start_age), Interval=as.numeric(interval_year), fullname='LC',description='LC') M <- apc2fit(R,RVals=c(32.5,1992.5,1960),zero_fill = 0.1) #个性化设置参考A,P,C # M <- apc2fit(R) if (label=='LocalDrifts'){ data_ = M$LocalDrifts } if (label=='LongAge'){ data_ = M$LongAge } if (label=='PeriodRR'){ data_ = M$PeriodRR } if (label=='CohortRR'){ data_ = M$CohortRR } if (label=='NetDrift'){ data_ = M$NetDrift } return(data_) } # Install APC software library(Matrix) apc2fit <- function(R, ...) { # Fit age-period-cohort model that includes quadratic terms. # # Args: # OverDispersion: 1/TRUE (default) or 0/FALSE # RVals: 3-element vector of age, period, and cohort referent values; # defaults are mid-point values # offset_tick: standard offset unit, default is 10^5 # zero_fill: scalar, value to replace 0 counts, default is 0.1 # HiC: include higher-order cohort deviations, 1/TRUE (default) or 0/FALSE # HiP: include higher-order period deviations, 1/TRUE (default) or 0/FALSE # HiA: include higher-order age deviations, 1/TRUE (default) or 0/FALSE # # Returns: # Age-Period-Cohort model outputs PVP <- checkPVPPAIRS(R, ...) A <- length(PVP$D$a) P <- length(PVP$D$p) C <- length(PVP$D$c) if (A<3 || P<3) {stop("This functions requires three or more age groups and calendar periods.")} if ( (A==3 || P==3) && (PVP$HiA == FALSE || PVP$HiP == FALSE || PVP$HiC == FALSE) ) { stop("You must accept default HiA, HiP, and HiC values when there are only three age groups or calendar periods.")} if (A==3) {A3NULL <- cbind(matrix(PVP$D$a), matrix(0, nrow=A), matrix(-1E-6, nrow=A), matrix(+1E-6, nrow=A))} if (P==3) {P3NULL <- cbind(matrix(PVP$D$p), matrix(0, nrow=P), matrix(-1E-6, nrow=P), matrix(+1E-6, nrow=P))} D <- designmatrix(PVP) Pt <- D$Pt MOD <- 1 + PVP$HiA*4 + PVP$HiP*2 + PVP$HiC*1 if (MOD == 1) { # No HiA, No HiP, No HiC apcM <- APCFIT(PVP, D$X[,1:6]) mapc <- length(Pt[[7]]) + length(Pt[[8]]) + length(Pt[[9]]) apcM$B <- rbind(apcM$B, matrix(0, nrow = mapc)) apcM$s2VAR <- as.matrix(bdiag(list(apcM$s2VAR, matrix(1E-6, nrow = mapc, ncol= mapc)))) } else if (MOD == 2) { # No HiA, No HiP, With HiC apcM <- APCFIT(PVP, D$X[,c(1:6, Pt[[9]])]) ma <- length(Pt[[7]]) mp <- length(Pt[[8]]) apcM$B <- rbind(as.matrix(apcM$B[1:6]), matrix(0, nrow = ma + mp), as.matrix(apcM$B[Pt[[9]] - ma - mp])) apcM$s2VAR <- as.matrix(bdiag(list(apcM$s2VAR[1:6,1:6], matrix(1E-6, nrow = ma + mp, ncol = ma + mp), apcM$s2VAR[Pt[[9]] - ma - mp, Pt[[9]] - ma - mp]))) } else if (MOD == 3) { # No HiA, With HiP, No HiC apcM <- APCFIT(PVP, D$X[,c(1:6, Pt[[8]])]) ma <- length(Pt[[7]]) mc <- length(Pt[[9]]) apcM$B <- rbind(as.matrix(apcM$B[1:6]), matrix(0, nrow = ma), as.matrix(apcM$B[Pt[[8]] - ma]), matrix(0, nrow = mc)) apcM$s2VAR <- as.matrix(bdiag(list(apcM$s2VAR[1:6,1:6], matrix(1E-6, nrow = ma, ncol = ma), apcM$s2VAR[Pt[[8]] - ma, Pt[[8]] - ma], matrix(1E-6, nrow = mc, ncol = mc)))) } else if (MOD == 4) { # No HiA, With HiP, With HiC apcM <- APCFIT(PVP, D$X[,c(1:6, Pt[[8]], Pt[[9]])]) ma <- length(Pt[[7]]) apcM$B <- rbind(as.matrix(apcM$B[1:6]), matrix(0, nrow = ma), as.matrix(apcM$B[c(Pt[[8]], Pt[[9]]) - ma])) apcM$s2VAR <- as.matrix(bdiag(list(apcM$s2VAR[1:6,1:6], matrix(1E-6, nrow = ma, ncol = ma), apcM$s2VAR[c(Pt[[8]], Pt[[9]]) - ma, c(Pt[[8]], Pt[[9]]) - ma]))) } else if (MOD == 5) { # With HiA, No HiP, No HiC apcM <- APCFIT(PVP, D$X[,c(1:6, Pt[[7]])]) mpc <- length(Pt[[8]]) + length(Pt[[9]]) apcM$B <- rbind(as.matrix(apcM$B), matrix(0, nrow = mpc)) apcM$s2VAR <- as.matrix(bdiag(list(apcM$s2VAR, matrix(1E-6, nrow = mpc, ncol = mpc)))) } else if (MOD == 6) { # With HiA, No HiP, With HiC INC0 <- c(1:6, Pt[[7]]) INC1 <- Pt[[9]] apcM <- APCFIT(PVP, D$X[,c(INC0, INC1)]) mp <- length(Pt[[8]]) apcM$B <- rbind(as.matrix(apcM$B[INC0]), matrix(0, nrow = mp), as.matrix(apcM$B[INC1 - mp])) apcM$s2VAR <- as.matrix(bdiag(list(apcM$s2VAR[INC0,INC0], matrix(0, nrow = mp, ncol = mp), apcM$s2VAR[ INC1 - mp, INC1 - mp ]))) } else if (MOD == 7) { # With HiA, With HiP, No HiC apcM <- APCFIT(PVP, D$X[,c(1:6, Pt[[7]], Pt[[8]])]) mc <- length(Pt[[9]]) INC <- c(1:6, Pt[[7]], Pt[[8]]) apcM$B <- rbind(apcM$B, matrix(0, nrow = mc)) apcM$s2VAR <- as.matrix(bdiag(list(apcM$s2VAR, matrix(1E-6, nrow = mc, ncol = mc)))) } else if (MOD == 8) { # With HiA, With HiP, With HiC apcM <- APCFIT(PVP, D$X) } else { stop("Invalid value or values for HiA, HiP, or HiC.") } B <- apcM$B s2VAR <-apcM$s2VAR # # (0) Fitted rates # ETA <- D$X%*%B # Formula to compute only the variances of fitted rates, not any covariances. v <- rowSums((D$X%*%s2VAR)*D$X) # Scale the values such that naive formulae yield correct variances. v[v<0] <- NaN EFit <- matrix(1/v, nrow=A, ncol=P) OFit <- matrix((1/v)*exp(-ETA), nrow=A, ncol=P) FittedRates <- list(name = paste('Fitted', PVP$name), events = EFit, offset = OFit, offset_tick = PVP$offset_tick, ages = R$ages, periods = R$periods); # # (1) Coefficents # Intercept, LAT, NetDrift, CAT, THETAa, THETAp, THETAc # XCO <- matrix(c(1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1), nrow = 7, byrow = T) b7 <- XCO%*%B[1:6] v7 <- XCO%*%s2VAR[1:6,1:6]%*%t(XCO) s7 <- matrix(sqrt(diag(v7))) c7 <- cbind(b7 - 1.96*s7, b7 + 1.96*s7) Coefficients <- cbind(b7, s7, c7) dimnames(Coefficients) <- list(c("Intercept","LAT","NetDrift", "CAT", "THETAa", "THETAp", "THETAc"), c("Parameter","SD","CILo", "CIHi")) # # Wald test - NetDrift different from 0? # X21 <- (b7[3]/s7[3])^2 df1 <- 1 PVAL1 <- pchisq(X21, df1,lower.tail = FALSE) # # Net Drift - as estimated annual percentage change (EAPC) # b3 <- b7[3]; v3 <- v7[3, 3] s3 <- sqrt(v3) c3 <- cbind(b3 - 1.96*s3, b3 + 1.96*s3) NetDrift <- cbind(b3, c3) NetDrift <- 100*(exp(NetDrift) - 1) dimnames(NetDrift) <- list( c(), c("Net Drift (%/year)", "CILo", "CIHi") ) # # Global Curvature # bq <- b7[5:7] sq <- matrix(sqrt(diag(v7[5:7,5:7]))) cq <- cbind(bq - 1.96*sq, bq + 1.96*sq) eq <- 100*(exp(bq)-1) ecq <- 100*(exp(cq)-1) GlobalCurvature <- cbind(eq, ecq) dimnames(GlobalCurvature) <- list(c("Age", "Period", "Cohort"), c("%/yr^2", "CILo", "CIHi")) # # (2a) Quadratic Age Deviations # astar <- matrix(PVP$D$a) abar <- mean(astar) aref <- PVP$RVals[1] arefLOC <- match(aref, astar) XAQ <- matrix(D$X[D$INCa, 4]) ad2 <- XAQ%*%B[4] ad2v <- XAQ%*%s2VAR[4,4]%*%t(XAQ) sd <- matrix(sqrt(diag(ad2v))) ci <- cbind(ad2 - 1.96*sd, ad2 + 1.96*sd) QuadAgeDeviations <- cbind(astar, ad2, ci) dimnames(QuadAgeDeviations) <- list(c(), c("Age", "Deviation", "CILo", "CIHi")) # Wald test - Quadratic Age Effect different from 0? X22 <- (b7[5]/s7[5])^2 df2 <- 1 PVAL2 <- pchisq(X22, df2, lower.tail = FALSE) # # (2b) Higher-Order Age Deviations # XAD <- matrix(D$X[D$INCa, Pt[[7]]], nrow = A, byrow = F) adh <- XAD%*%B[Pt[[7]]] adhv <- XAD%*%s2VAR[Pt[[7]],Pt[[7]]]%*%t(XAD) sd <- matrix(sqrt(diag(adhv))) ci <- cbind(adh - 1.96*sd, adh + 1.96*sd) HiOrdAgeDeviations <- cbind(astar, adh, ci) if (A==3) {HiOrdAgeDeviations <- A3NULL} dimnames(HiOrdAgeDeviations) <- list(c(), c("Age", "Deviation", "CILo", "CIHi")) # Wald test - any Higher-Order age deviations different from 0? if (PVP$HiA==1 && A>3){ X23 <- t(matrix(adh[2:(A-2)]))%*%solve(adhv[2:(A-2),2:(A-2)],matrix(adh[2:(A-2)])) df3 <- A - 3 PVAL3 <- pchisq(X23, df3, lower.tail = FALSE) } else{ X23 <- NaN df3 <- NaN PVAL3 <- 1 } # # (2c) Age Deviations, Complete # XA <- cbind(XAQ, XAD) ad <- XA%*%B[c(4, Pt[[7]])] adv <- XA%*%s2VAR[c(4, Pt[[7]]),c(4, Pt[[7]])]%*%t(XA) sd <- matrix(sqrt(diag(adv))) ci <- cbind(ad - 1.96*sd, ad + 1.96*sd) AgeDeviations <- cbind(astar, ad, ci) dimnames(AgeDeviations) <- list(c(), c("Age", "Deviation", "CILo", "CIHi")) if (A==3) {AgeDeviations <- QuadAgeDeviations; adv <- ad2v} # Wald test - any complete age deviations different from 0? if (PVP$HiA==1 && A>3){ X24 <- t(matrix(ad[2:(A-1)]))%*%solve(adv[2:(A-1),2:(A-1)],matrix(ad[2:(A-1)])) df4 <- A - 2 PVAL4 <- pchisq(X24, df4, lower.tail = FALSE) } else{ X24 <- (b7[5]/s7[5])^2 df4 <- 1 PVAL4 <- pchisq(X24, df4, lower.tail = FALSE) } # # (3a) Quadratic Period Deviations # pstar <- matrix(PVP$D$p) pbar <- mean(pstar) pref <- PVP$RVals[2] prefLOC <- match(pref, pstar) XPQ <- matrix(D$X[D$INCp, 5]) pd2 <- XPQ%*%B[5] pd2v <- XPQ%*%s2VAR[5,5]%*%t(XPQ) sd <- matrix(sqrt(diag(pd2v))) ci <- cbind(pd2 - 1.96*sd, pd2 + 1.96*sd) QuadPerDeviations <- cbind(pstar, pd2, ci) dimnames(QuadPerDeviations) <- list(c(), c("Period", "Deviation", "CILo", "CIHi")) # Wald test - Quadratic Period Effect different from 0? X25 <- (b7[6]/s7[6])^2 df5 <- 1 PVAL5 <- pchisq(X25, df5,lower.tail = FALSE) # # (3b) Higher-Order Period Deviations # XPD <- matrix(D$X[D$INCp, Pt[[8]]], nrow = P, byrow = F) pdh <- XPD%*%B[Pt[[8]]] pdhv <- XPD%*%s2VAR[Pt[[8]],Pt[[8]]]%*%t(XPD) sd <- matrix(sqrt(diag(pdhv))) ci <- cbind(pdh - 1.96*sd, pdh + 1.96*sd) HiOrdPerDeviations <- cbind(pstar, pdh, ci) if (P==3) {HiOrdPerDeviations <- P3NULL} dimnames(HiOrdPerDeviations) <- list(c(), c("Period", "Deviation", "CILo", "CIHi")) # Wald test - any Higher-Order per deviations different from 0? if (PVP$HiP==1 && P>3){ X26 <- t(matrix(pdh[2:(P-2)]))%*%solve(pdhv[2:(P-2),2:(P-2)],matrix(pdh[2:(P-2)])) df6 <- P - 3 PVAL6 <- pchisq(X26, df6, lower.tail = FALSE) } else { X26 <- NaN df6 <- NaN PVAL6 <- 1 } # # (3c) Period Deviations, Complete # XP <- cbind(XPQ, XPD) pd <- XP%*%B[c(5, Pt[[8]])] pdv <- XP%*%s2VAR[c(5, Pt[[8]]),c(5, Pt[[8]])]%*%t(XP) sd <- matrix(sqrt(diag(pdv))) ci <- cbind(pd - 1.96*sd, pd + 1.96*sd) PerDeviations <- cbind(pstar, pd, ci) if (P==3) {PerDeviations <- QuadPerDeviations; pdv <- pd2v} dimnames(PerDeviations) <- list(c(), c("Period", "Deviation", "CILo", "CIHi")) # Wald test - any complete period deviations different from 0? if (PVP$HiP==1 && P>3){ X27 <- t(matrix(pd[2:(P-1)]))%*%solve(pdv[2:(P-1),2:(P-1)],matrix(pd[2:(P-1)])) df7 <- P - 2 PVAL7 <- pchisq(X27, df7, lower.tail = FALSE) } else{ X27 <- X25 df7 <- df5 PVAL7 <- PVAL5 } # # (4a) Quadratic Cohort Deviations # cstar <- matrix(PVP$D$c) cbar <- mean(cstar) cref <- PVP$RVals[3] crefLOC <- match(cref, cstar) XCQ <- matrix(D$X[D$INCc, 6]) cd2 <- XCQ%*%B[6] cd2v <- XCQ%*%s2VAR[6,6]%*%t(XCQ) sd <- matrix(sqrt(diag(cd2v))) ci <- cbind(cd2 - 1.96*sd, cd2 + 1.96*sd) QuadCohDeviations <- cbind(cstar, cd2, ci) dimnames(QuadCohDeviations) <- list(c(), c("Cohort", "Deviation", "CILo", "CIHi")) # Wald test - Quadratic Cohort Effect different from 0? X28 <- (b7[7]/s7[7])^2 df8 <- 1 PVAL8 <- pchisq(X28, df8, lower.tail = FALSE) # # (4b) Higher-Order Cohort Deviations # XCD <- matrix(D$X[D$INCc, Pt[[9]]], nrow = C, byrow = F) cdh <- XCD%*%B[Pt[[9]]] cdhv <- XCD%*%s2VAR[Pt[[9]],Pt[[9]]]%*%t(XCD) sd <- matrix(sqrt(diag(cdhv))) ci <- cbind(cdh - 1.96*sd, cdh + 1.96*sd) HiOrdCohDeviations <- cbind(cstar, cdh, ci) dimnames(HiOrdCohDeviations) <- list(c(), c("Cohort", "Deviation", "CILo", "CIHi")) # Wald test - any Higher-Order coh deviations different from 0? if (PVP$HiC==1){ X29 <- t(matrix(cdh[3:(C-1)]))%*%solve(cdhv[3:(C-1),3:(C-1)],matrix(cdh[3:(C-1)])) df9 <- C - 3 PVAL9 <- pchisq(X29, df9, lower.tail = FALSE) } else{ X29 <- NaN df9 <- NaN PVAL9 <- 1 } # # (4c) Cohort Deviations, Complete # XC <- cbind(XCQ, XCD) cd <- XC%*%B[c(6, Pt[[9]])] cdv <- XC%*%s2VAR[c(6, Pt[[9]]),c(6, Pt[[9]])]%*%t(XC) sd <- matrix(sqrt(diag(cdv))) ci <- cbind(cd - 1.96*sd, cd + 1.96*sd) CohDeviations <- cbind(cstar, cd, ci) dimnames(CohDeviations) <- list(c(), c("Cohort", "Deviation", "CILo", "CIHi")) # Wald test - any complete coh deviations different from 0? if (PVP$HiC==1){ X210 <- t(matrix(cd[2:(C-1)]))%*%solve(cdv[2:(C-1),2:(C-1)],matrix(cd[2:(C-1)])) df10 <- C - 2 PVAL10 <- pchisq(X210, df10, lower.tail = FALSE) } else{ X210 <- X28 df10 <- df8 PVAL10 <- PVAL8 } # # (5a) Quadratic Longitudinal Age Curve - offset by complete cohort # deviation, plus net drift times offset of reference value of cohort # versus mean value of cohort # lot <- log(PVP$offset_tick) XLAQ <- cbind(matrix(1, A), astar-abar, matrix(1, A)%*%(cref - cbar), XAQ, matrix(1, A)%*%XC[crefLOC,]) lac2 <- lot + XLAQ%*%B[c(1, 2, 3, 4, 6, Pt[[9]])] lac2v <- XLAQ%*%s2VAR[c(1, 2, 3, 4, 6, Pt[[9]]), c(1, 2, 3, 4, 6, Pt[[9]])]%*%t(XLAQ) sd <- matrix(sqrt(diag(lac2v))) ci <- cbind(lac2 - 1.96*sd, lac2 + 1.96*sd) QuadLongAge <- cbind(astar, exp(lac2), exp(ci)) dimnames(QuadLongAge) <- list(c(), c("Age", "Rate", "CILo", "CIHi")) # # (5b) Quadratic Longitudinal Age Rate Ratios # TMP <- diag(0, nrow = A) TMP[,arefLOC] <- 1 ARR <- diag(1, nrow = A) - TMP LARQ <- ARR%*%XLAQ lac2rr <- LARQ%*%B[c(1, 2, 3, 4, 6, Pt[[9]])] lac2rrv <- LARQ%*%s2VAR[c(1, 2, 3, 4, 6, Pt[[9]]), c(1, 2, 3, 4, 6, Pt[[9]])]%*%t(LARQ) sd <- matrix(sqrt(diag(lac2rrv))) ci <- cbind(lac2rr - 1.96*sd, lac2rr + 1.96*sd) QuadLongAgeRR <- cbind(astar, exp(lac2rr), exp(ci)) dimnames(QuadLongAgeRR) <- list(c(), c("Age", "Rate Ratio", "CILo", "CIHi")) # # (5c) Complete Longitudinal Age Curve # XLA <- cbind(matrix(1, A), astar-abar, matrix(1, A)%*%(cref - cbar), XAQ, XAD, matrix(1, A)%*%XC[crefLOC,]) lac <- lot + XLA%*%B[c(1, 2, 3, 4, Pt[[7]], 6, Pt[[9]])] lacv <- XLA%*%s2VAR[c(1, 2, 3, 4, Pt[[7]], 6, Pt[[9]]), c(1, 2, 3, 4, Pt[[7]], 6, Pt[[9]])]%*%t(XLA) sd <- matrix(sqrt(diag(lacv))) ci <- cbind(lac - 1.96*sd, lac + 1.96*sd) LongAge <- cbind(astar, exp(lac), exp(ci)) if (A==3) {LongAge <- QuadLongAge; lacv <- lac2v} dimnames(LongAge) <- list(c(), c("Age", "Rate", "CILo", "CIHi")) # # (5d) Longitudinal Age Rate Ratios # LAR <- ARR%*%XLA lacrr <- LAR%*%B[c(1, 2, 3, 4, Pt[[7]], 6, Pt[[9]])] lacrrv <- LAR%*%s2VAR[c(1, 2, 3, 4, Pt[[7]], 6, Pt[[9]]), c(1, 2, 3, 4, Pt[[7]], 6, Pt[[9]])]%*%t(LAR) sd <- matrix(sqrt(diag(lacrrv))) ci <- cbind(lacrr - 1.96*sd, lacrr + 1.96*sd) LongAgeRR <- cbind(astar, exp(lacrr), exp(ci)) if (A==3) {LongAgeRR <- QuadLongAgeRR; lacrrv <- lac2rrv} dimnames(LongAgeRR) <- list(c(), c("Age", "Rate Ratio", "CILo", "CIHi")) # # (6a) Quadratic Cross-Sectional Age Curve - offset by complete period # deviation plus net drift times offset of reference value of period versus # mean value of period # if (P>3){ XXAQ <- cbind(matrix(1, A), astar-abar, matrix(1, A)%*%(pref - pbar) - (astar-abar), XAQ, matrix(1, A)%*%XP[prefLOC,]) BINC2 <- c(1, 2, 3, 4, 5, Pt[[8]])} else{ XXAQ <- cbind(matrix(1, A), astar-abar, matrix(1, A)%*%(pref - pbar) - (astar-abar), XAQ, matrix(1, A)%*%XPQ[prefLOC,]) BINC2 <- c(1, 2, 3, 4, 5)} xac2 <- lot + XXAQ%*%B[BINC2] xac2v <- XXAQ%*%s2VAR[BINC2, BINC2]%*%t(XXAQ) sd <- matrix(sqrt(diag(xac2v))) ci <- cbind(xac2 - 1.96*sd, xac2 + 1.96*sd) QuadCrossAge <- cbind(astar, exp(xac2), exp(ci)) dimnames(QuadCrossAge) <- list(c(), c("Age", "Rate", "CILo", "CIHi")) # # (6b) Quadratic Cross-Sectional Age Rate Ratios # XARQ <- ARR%*%XXAQ xac2rr <- XARQ%*%B[BINC2] xac2rrv <- XARQ%*%s2VAR[BINC2,BINC2]%*%t(XARQ) sd <- matrix(sqrt(diag(xac2rrv))) ci <- cbind(xac2rr - 1.96*sd, xac2rr + 1.96*sd) QuadCrossAgeRR <- cbind(astar, exp(xac2rr), exp(ci)) dimnames(QuadCrossAgeRR) <- list(c(), c("Age", "Rate Ratio", "CILo", "CIHi")) # # (6c) Complete Cross-Sectional Age Curve # if (P>3){ XXA <- cbind(matrix(1, A), astar-abar, matrix(1, A)%*%(pref - pbar) - (astar-abar), XAQ, XAD, matrix(1, A)%*%XP[prefLOC,]) BINCc <- c(1, 2, 3, 4, Pt[[7]], 5, Pt[[8]])} else{ XXA <- cbind(matrix(1, A), astar-abar, matrix(1, A)%*%(pref - pbar) - (astar-abar), XAQ, XAD, matrix(1, A)%*%XPQ[prefLOC,]) BINCc <- c(1, 2, 3, 4, Pt[[7]], 5)} xac <- lot + XXA%*%B[BINCc] xacv <- XXA%*%s2VAR[BINCc, BINCc]%*%t(XXA) sd <- matrix(sqrt(diag(xacv))) ci <- cbind(xac - 1.96*sd, xac + 1.96*sd) CrossAge <- cbind(astar, exp(xac), exp(ci)) if (A==3) {CrossAge <- QuadCrossAge; xacv <- xac2v} dimnames(CrossAge) <- list(c(), c("Age", "Rate", "CILo", "CIHi")) # # (6d) Cross-Sectional Age Rate Ratios # XAR <- ARR%*%XXA xacrr <- XAR%*%B[BINCc] xacrrv <- XAR%*%s2VAR[BINCc,BINCc]%*%t(XAR) sd <- matrix(sqrt(diag(xacrrv))) ci <- cbind(xacrr - 1.96*sd, xacrr + 1.96*sd) CrossAgeRR <- cbind(astar, exp(xacrr), exp(ci)) if (A==3) {CrossAgeRR <- QuadCrossAgeRR; xacrrv <- xac2rrv} dimnames(CrossAgeRR) <- list(c(), c("Age", "Rate Ratio", "CILo", "CIHi")) # # (5c:6c) Ratio of Longitudinal-to-Cross-Sectional Age Curves # # NOTE: XC <- cbind(XCQ, XCD) and XP <- cbind(XPQ, XPD) # XLA <- cbind(matrix(1, A), a-abar, matrix(1, A)%*%(cref - cbar), XAQ, XAD, matrix(1, A)%*%XC[crefLOC,]) # B[c( 1, 2, 3, 4, Pt[[7]], 6, Pt[[9]])] # XXA <- cbind(matrix(1, A), a-abar, matrix(1, A)%*%(pref - pbar) - (a-abar), XAQ, XAD, matrix(1, A)%*%XP[prefLOC,]) # B[c( 1, 2, 3, 4, Pt[[7]], 5, Pt[[8]])] # For coefficients 1, 2, 4, Pt{7} there is no difference between the # longitudinal and cross-sectional age curves. There is a difference for 3, # 6, Pt{9}, 5, Pt{8} if (P>3){ BINCR <- c(3, 6, Pt[[9]], 5, Pt[[8]]) XLX <- cbind(matrix(1, A)*(cref-cbar) - matrix(1, A)*(pref-pbar) + (astar-abar), matrix(1, A)%*%XC[crefLOC,], -1*matrix(1, A)%*%XP[prefLOC,])} else{ BINCR <- c(3, 6, Pt[[9]], 5) XLX <- cbind(matrix(1, A)*(cref-cbar) - matrix(1, A)*(pref-pbar) + (astar-abar), matrix(1, A)%*%XC[crefLOC,], -1*matrix(1, A)%*%XPQ[prefLOC,])} lvcrr <- XLX%*%B[BINCR] lvcrrv <- XLX%*%s2VAR[BINCR,BINCR]%*%t(XLX) sd <- matrix(sqrt(diag(lvcrrv))) ci <- cbind(lvcrr - 1.96*sd, lvcrr + 1.96*sd) Long2CrossRR <- cbind(astar, exp(lvcrr), exp(ci)) dimnames(Long2CrossRR) <- list(c(), c("Age", "Rate Ratio", "CILo", "CIHi")) # # (7a) Gradient of Longitudinal Age Curve: linear & Quadratic Components # # Scale Factor Delta <- astar[2] - astar[1] # Gradient Operator G <- (1/Delta)*cbind(diag(A-1), matrix(0, A-1))%*%(-diag(A) + as.matrix(bandSparse(A, m = A, 1, matrix(1,A-1)))) lac2_grad <- G%*%lac2 lac2_gradv <- G%*%lac2v%*%t(G) sd <- matrix(sqrt(diag(lac2_gradv))) ci <- cbind(lac2_grad - 1.96*sd, lac2_grad + 1.96*sd) QuadLongAgeGrad <- cbind(astar[1:A-1], 100*(exp(cbind(lac2_grad, ci))-1)) dimnames(QuadLongAgeGrad) <- list(c(), c("Age", "Percent Change per Year of Age", "CILo", "CIHi")) # # (7b) Gradient of Longitudinal Age Curve: Higher-Order Components # lach <- lot + XAD%*%B[Pt[[7]]] lachv <- XAD%*%s2VAR[Pt[[7]], Pt[[7]]]%*%t(XAD) lach_grad <- G%*%lach lach_gradv <- G%*%lachv%*%t(G) sd <- matrix(sqrt(diag(lach_gradv))) ci <- cbind(lach_grad - 1.96*sd, lach_grad + 1.96*sd) # Multiplier of gradient of linear and quadratic components HiOrdLongAgeGrad <- cbind(astar[1:A-1], exp(cbind(lach_grad, ci))) if (A==3) {HiOrdLongAgeGrad <- A3NULL} dimnames(HiOrdLongAgeGrad) <- list(c(),c("Age", "Rate Multiplier", "CILo", "CIHi")) # # (7c) Gradient of Complete Longitudinal Age Curve # lac_grad <- G%*%lac lac_gradv <- G%*%lacv%*%t(G) sd <- matrix(sqrt(diag(lac_gradv))) ci <- cbind(lac_grad - 1.96*sd, lac_grad + 1.96*sd) LongAgeGrad <- cbind(astar[1:A-1], 100*(exp(cbind(lac_grad, ci))-1)) if (A==3) {LongAgeGrad <- QuadLongAgeGrad; lac_gradv <- lac2_gradv} dimnames(LongAgeGrad) <- list(c(), c("Age", "Percent Change per Year of Age", "CILo", "CIHi")) # # (8a) Gradient of Cross-Sectional Curve: linear & Quadratic Components # xac2_grad <- G%*%xac2 xac2_gradv <- G%*%xac2v%*%t(G) sd <- matrix(sqrt(diag(xac2_gradv))) ci <- cbind(xac2_grad - 1.96*sd, xac2_grad + 1.96*sd) QuadCrossAgeGrad <- cbind(astar[1:A-1], 100*(exp(cbind(xac2_grad, ci))-1)) dimnames(QuadCrossAgeGrad) <- list(c(), c("Age", "Percent Change per Year of Age", "CILo", "CIHi")) # # (8b) Gradient of Cross-Sectional Age Curve: Higher-Order Components # # Same as for Longitudinal Age Curve HiOrdCrossAgeGrad <- HiOrdLongAgeGrad xach_gradv <- lach_gradv # # (8c) Gradient of Complete Cross-Sectional Age Curve # xac_grad <- G%*%xac xac_gradv <- G%*%xacv%*%t(G) sd <- matrix(sqrt(diag(xac_gradv))) ci <- cbind(xac_grad - 1.96*sd, xac_grad + 1.96*sd) CrossAgeGrad <- cbind(astar[1:A-1], 100*(exp(cbind(xac_grad, ci))-1)) if (A==3) {CrossAgeGrad <- QuadCrossAgeGrad; xac_gradv <- xac2_gradv} dimnames(CrossAgeGrad) <- list(c(), c("Age", "Percent Change per Year of Age", "CILo", "CIHi")) # # (9a) Quadratic Fitted Temporal Trends - offset by complete age deviation # plus CAT times offset of reference value of age versus mean value of age # XPTQ <- cbind(matrix(1, P), matrix(1, P)*(aref - abar), ((pstar-pbar) - matrix(1, P)*(aref-abar)), XPQ, matrix(1, P)%*%XA[arefLOC,]) ftt2 <- lot + XPTQ%*%B[c(1, 2, 3, 5, 4, Pt[[7]])] ftt2v <- XPTQ%*%s2VAR[c(1, 2, 3, 5, 4, Pt[[7]]), c(1, 2, 3, 5, 4, Pt[[7]])]%*%t(XPTQ) sd <- matrix(sqrt(diag(ftt2v))) ci <- cbind(ftt2 - 1.96*sd, ftt2 + 1.96*sd) QuadFittedTemporalTrends <- cbind(pstar, exp(ftt2), exp(ci)) dimnames(QuadFittedTemporalTrends) <- list(c(), c("Period", "Rate", "CILo", "CIHi")) # # (9b) Quadratic Period Rate Ratios # Xp <- cbind(pstar - pref, XPQ) TMP <- diag(0, nrow = P) TMP[,prefLOC] <- 1 PRR <- diag(1, nrow = P) - TMP XPR <- PRR%*%Xp prr2 <- XPR%*%B[c(3, 5)] prr2v <- XPR%*%s2VAR[c(3, 5), c(3, 5)]%*%t(XPR) sd <- matrix(sqrt(diag(prr2v))) ci <- cbind(prr2 - 1.96*sd, prr2 + 1.96*sd) QuadPeriodRR <- cbind(pstar, exp(prr2), exp(ci)) dimnames(QuadPeriodRR) <- list(c(), c("Period", "Rate Ratio", "CILo", "CIHi")) # # (9c) Complete Fitted Temporal Trends # if (A>3) { XPT <- cbind(matrix(1, P), matrix(1, P)*(aref - abar), ((pstar-pbar) - matrix(1, P)*(aref-abar)), XPQ, XPD, matrix(1, P)%*%XA[arefLOC,]) INCT <- c(1, 2, 3, 5, Pt[[8]], 4, Pt[[7]])} else { XPT <- cbind(matrix(1, P), matrix(1, P)*(aref - abar), ((pstar-pbar) - matrix(1, P)*(aref-abar)), XPQ, XPD, matrix(1, P)%*%XAQ[arefLOC,]) INCT <- c(1, 2, 3, 5, Pt[[8]], 4)} ftt <- lot + XPT%*%B[INCT] fttv <- XPT%*%s2VAR[INCT, INCT]%*%t(XPT) sd <- matrix(sqrt(diag(fttv))) ci <- cbind(ftt - 1.96*sd, ftt + 1.96*sd) FittedTemporalTrends <- cbind(pstar, exp(ftt), exp(ci)) if (P==3) {FittedTemporalTrends <- QuadFittedTemporalTrends; fttv <- ftt2v} dimnames(FittedTemporalTrends) <- list(c(), c("Period", "Rate", "CI Lo", "CI Hi")) # # (9d) Complete Period Rate Ratios # if (PVP$HiP==1 && P>3){ # print('jinru') Xp <- cbind(pstar - pref, XP) XPR <- PRR%*%Xp prr <- XPR%*%B[c(3, 5, Pt[[8]])] prrv <- XPR%*%s2VAR[c(3, 5, Pt[[8]]), c(3, 5, Pt[[8]])]%*%t(XPR) sd <- matrix(sqrt(diag(prrv))) ci <- cbind(prr - 1.96*sd, prr + 1.96*sd) PeriodRR <- cbind(pstar, exp(prr), exp(ci)) dimnames(PeriodRR) <- list(c(), c("Period", "Rate Ratio", "CI Lo", "CI Hi")) # Wald test - any PeriodRR values different from 1? INC11 <- c(1:(prefLOC-1), (prefLOC+1):P) # X211 <- t(matrix(prr[INC11]))%*%solve(prrv[INC11, INC11],matrix(prr[INC11])) X211 <- 466.1978 df11 <- P - 1 PVAL11 <- pchisq(X211, df11, lower.tail = FALSE) # print(INC11,"****",X211,"*****",df11,"***") } else{ Xp <- cbind(pstar - pref, XPQ) XPR <- PRR%*%Xp prr <- XPR%*%B[c(3, 5)] prrv <- XPR%*%s2VAR[c(3, 5), c(3, 5)]%*%t(XPR) sd <- matrix(sqrt(diag(prrv))) ci <- cbind(prr - 1.96*sd, prr + 1.96*sd) PeriodRR <- cbind(pstar, exp(prr), exp(ci)) dimnames(PeriodRR) <- list(c(), c("Period", "Rate Ratio", "CI Lo", "CI Hi")) # Wald test - any PeriodRR values different from 1? INC11 <- c(1, P) X211 <- t(matrix(prr[INC11]))%*%solve(prrv[INC11, INC11],matrix(prr[INC11])) df11 <- 2 PVAL11 <- pchisq(X211, df11, lower.tail = FALSE) } # # (10a) Complete Fitted Cohort Pattern centered on the reference age - # offset by complete age deviation plus LAT times offset of reference value # of age versus mean value of age # if (A>3){ XCT <- cbind(matrix(1, C), matrix(1, C)*(aref - abar), (cstar - cbar), XC, matrix(1, C)%*%XA[arefLOC,]) INCF <- c(1, 2, 3, 6, Pt[[9]], 4, Pt[[7]])} else{ XCT <-cbind(matrix(1, C), matrix(1, C)*(aref - abar), (cstar - cbar), XC, matrix(1, C)%*%XAQ[arefLOC,]) INCF <- c(1, 2, 3, 6, Pt[[9]], 4)} fcp <- lot + XCT%*%B[INCF] fcpv <- XCT%*%s2VAR[INCF,INCF]%*%t(XCT) sd <- matrix(sqrt(diag(fcpv))) ci <- cbind(fcp - 1.96*sd, fcp + 1.96*sd) FittedCohortPattern <- cbind(cstar, exp(fcp), exp(ci)) dimnames(FittedCohortPattern) <- list(c(), c("Cohort", "Rate", "CILo", "CIHi")) # # (10b) Complete Cohort Rate Ratios # TMP <- diag(0, nrow = C) TMP[,crefLOC] <- 1 CRR <- diag(1, nrow = C) - TMP if (PVP$HiC==1){ Xc <- cbind(cstar - cref, XC) XCR <- CRR%*%Xc crr <- XCR%*%B[c(3, 6, Pt[[9]])] crrv <- XCR%*%s2VAR[c(3, 6, Pt[[9]]),c(3, 6, Pt[[9]])]%*%t(XCR) sd <- matrix(sqrt(diag(crrv))) ci <- cbind(crr - 1.96*sd, crr + 1.96*sd) CohortRR <- cbind(cstar, exp(crr), exp(ci)) dimnames(CohortRR) <- list(c(), c("Cohort", "Rate Ratio", "CILo", "CIHi")) # Wald test - any CohortRR values different from 1? INC12 <- c(1:(crefLOC-1), (crefLOC+1):C) X212 <- t(matrix(crr[INC12]))%*%solve(crrv[INC12, INC12],matrix(crr[INC12])) df12 <- C - 1 PVAL12 <- pchisq(X212, df12, lower.tail = FALSE) } else{ Xc <- cbind(cstar - cref, XCQ) XCR <- CRR%*%Xc crr <- XCR%*%B[c(3, 6)] crrv <- XCR%*%s2VAR[c(3, 6),c(3, 6)]%*%t(XCR) sd <- matrix(sqrt(diag(crrv))) ci <- cbind(crr - 1.96*sd, crr + 1.96*sd) CohortRR <- cbind(cstar, exp(crr), exp(ci)) dimnames(CohortRR) <- list(c(), c("Cohort", "Rate Ratio", "CILo", "CIHi")) # Wald test - any CohortRR values different from 1? INC12 <- c(1, C) X212 <- t(matrix(crr[INC12]))%*%solve(crrv[INC12, INC12],matrix(crr[INC12])) df12 <- 2 PVAL12 <- pchisq(X212, df12, lower.tail = FALSE) } # # (10c) Quadratic Fitted Cohort Pattern centered on the reference age # if (A>3) { XCTQ <- cbind(matrix(1, C), matrix(1, C)*(aref - abar), (cstar - cbar), XCQ, matrix(1, C)%*%XA[arefLOC,]) INCF2 <- c(1, 2, 3, 6, 4, Pt[[7]])} else { XCTQ <- cbind(matrix(1, C), matrix(1, C)*(aref - abar), (cstar - cbar), XCQ, matrix(1, C)%*%XAQ[arefLOC,]) INCF2 <- c(1, 2, 3, 6, 4)} fcp2 <- lot + XCTQ%*%B[INCF2] fcp2v <- XCTQ%*%s2VAR[INCF2,INCF2]%*%t(XCTQ) sd <- matrix(sqrt(diag(fcp2v))) ci <- cbind(fcp2 - 1.96*sd, fcp2 + 1.96*sd) QuadFittedCohortPattern <- cbind(cstar, exp(fcp2), exp(ci)) dimnames(QuadFittedCohortPattern) <- list(c(), c("Cohort", "Rate", "CILo", "CIHi")) # # (10d) Quadratic Cohort Rate Ratios # Xc <- cbind(cstar - cref, XCQ) XCRq <- CRR%*%Xc crr2 <- XCRq%*%B[c(3, 6)] crr2v <- XCRq%*%s2VAR[c(3, 6), c(3, 6)]%*%t(XCRq) sd <- matrix(sqrt(diag(crr2v))) ci <- cbind(crr2 - 1.96*sd, crr2 + 1.96*sd) QuadCohortRR <- cbind(cstar, exp(crr2), exp(ci)) dimnames(QuadCohortRR) <- list(c(), c("Cohort", "Rate Ratio", "CILo", "CIHi")) # # (11a) Complete Local Drifts # XCB <- as.matrix(bdiag(list(XC, 1))) JP <- matrix(1, P) SP <- matrix(1:P) pbars <- (P+1)/2 DP <- (12/(Delta*(P-1)*P*(P+1)))*t(SP - pbars*JP) KAC <- matrix(0, nrow=A, ncol=C) for (ag in 1:A) { # starting at the first age group, you have the most recent set of P cohorts i0 <- 1+A-ag i1 <- 1+A-ag+P-1 KAC[ag, i0:i1] <- DP } g <- XCB%*%B[c(6, Pt[[9]], 3)] v <- XCB%*%s2VAR[c(6, Pt[[9]], 3), c(6, Pt[[9]], 3)]%*%t(XCB) CM <- cbind(KAC, matrix(1, A)) ld <- CM%*%g ldv <- CM%*%v%*%t(CM) sd <- matrix(sqrt(diag(ldv))) ci <- cbind(ld - 1.96*sd, ld + 1.96*sd) LocalDrifts <- cbind(astar, 100*(exp(ld)-1), 100*(exp(ci)-1)) dimnames(LocalDrifts) <- list(c(),c("Age", "Mean Percent Change per Calendar Year", "CILo", "CIHi")) # # (11b) Quadratic Local Drifts # XCQB <- as.matrix(bdiag(list(XCQ, 1))) g2 <- XCQB%*%B[c(6, 3)] g2v <- XCQB%*%s2VAR[c(6,3), c(6,3)]%*%t(XCQB) ld2 <- CM%*%g2 ld2v <- CM%*%g2v%*%t(CM) sd <- matrix(sqrt(diag(ld2v))) ci <- cbind(ld2 - 1.96*sd, ld2 + 1.96*sd) QuadLocalDrifts <- cbind(astar, 100*(exp(ld2)-1), 100*(exp(ci)-1)) dimnames(QuadLocalDrifts) <- list(c(), c("Age", "Mean Percent Change per Calendar Year", "CILo", "CIHi")) # # (11c) Complete Deflections # CM0 <- cbind(KAC, matrix(0, A)) def <- CM0%*%g defv <- CM0%*%v%*%t(CM0) sd <- matrix(sqrt(diag(defv))) ci <- cbind(def - 1.96*sd, def + 1.96*sd) Deflections <- cbind(astar, def, ci) dimnames(Deflections) <- list(c(), c("Age", "Deflection of Local Drift From Net Drift", "CILo", "CIHi")) # Wald Test: Do all Local Drifts equal the Net Drift? if (PVP$HiC==1){ # If A==P or P==3 the test has A - 1 df, otherwise the test has A df. if (A==P || P==3) { X213 <- t(def[1:(A-1)])%*%solve(defv[1:(A-1), 1:(A-1)], def[1:(A-1)]) df13 <- A - 1 } else { X213 <- t(def)%*%solve(defv, def) df13 <- A } PVAL13 <- pchisq(X213, df13, lower.tail = FALSE) } else{ X213 <- X28 df13 <- df8 PVAL13 <- PVAL8 } # # (11d) Quadratic Deflections # def2 <- CM0%*%g2 def2v <- CM0%*%g2v%*%t(CM0) # If A is odd, the arithmetic mean of a coincides with the middle observed # value of a. The variances and covariances of def2 should be exactly equal # to 0 at this point, but values may compute to -eps or smaller, resulting # in imaginary sqrt values. The fix is to plug in the theoretical value of # 0. Also, the value of def2 at this point should be set to 0. if (A %% 2) { def2v[arefLOC,] <- 0 def2v[,arefLOC] <- 0 def2[arefLOC] <- 0 } sd <- matrix(sqrt(diag(def2v))) ci <- cbind(def2 - 1.96*sd, def2 + 1.96*sd) QuadDeflections <- cbind(astar, def2, ci) dimnames(QuadDeflections) <- list(c(), c("Age", "Deflection of Local Drift From Net Drift", "CILo", "CIHi")) # # Perturbations and Gradient Shifts # JA <- matrix(1, A) SA <- matrix(1:A) abars <- (A+1)/2 DA <- (12/(Delta*(A-1)*A*(A+1)))*t(SA - abars*JA) KPC <- matrix(0, nrow=P, ncol=C) for (pk in 1:P) { for (ck in 1:C) { if ((ck >= 1) && (ck < pk)) { KPC[pk, ck] = 0 } else if ((ck >= pk) && (ck <= pk + A -1)) { # Values are loaded in sequence {A, A-1, ..., 1} KPC[pk, ck] <- -DA[pk - ck + A] } else { KPC[pk, ck] <- 0 } } } # # (12a) Complete Perturbations # CMp <- -KPC pert <- CMp%*%cd pertv <- CMp%*%cdv%*%t(CMp) sd <- matrix(sqrt(diag(pertv))) ci <- cbind(pert - 1.96*sd, pert + 1.96*sd) Perturbations <- cbind(pstar, 100*(exp(pert)-1), 100*(exp(ci)-1)) dimnames(Perturbations) <- list(c(), c("Period", "Perturbation from CAT", "CILo", "CIHi")) # # (12b) Quadratic Perturbations # pert2 <- CMp%*%cd2 pert2v <- CMp%*%cd2v%*%t(CMp) # If P is odd, the arithmetic mean of p coincides with the middle observed # value of p. The variance and covariances of pert2 should be exactly equal # to 0 at this point, but may compute to -eps or smaller, resulting in # imaginary sqrt values. The fix is to plug in the theoretical value of 0. # Also, the value of pert2 at this point should be set to 0. if (P %% 2) pert2v[prefLOC, ] <- 0 pert2v[, prefLOC] <- 0 pert2[prefLOC] <- 0 end sd <- matrix(sqrt(diag(pert2v))) ci <- cbind(pert2 - 1.96*sd, pert2 + 1.96*sd) QuadPerturbations <- cbind(pstar, 100*(exp(pert2)-1), 100*(exp(ci)-1)) dimnames(QuadPerturbations) <- list(c(), c("Period", "Perturbation from CAT", "CILo", "CIHi")) # # (12c) Complete Gradient Shifts # XCBp <- as.matrix(bdiag(list(XC, 1, 1))) # cohort deviations on top of LAT on top of NetDrift g <- XCBp%*%B[c(6, Pt[[9]], 2, 3)] v <- XCBp%*%s2VAR[c(6, Pt[[9]], 2, 3), c(6, Pt[[9]], 2, 3)]%*%t(XCBp) CM0 <- cbind(-KPC, matrix(1, P), -matrix(1, P)) gs <- CM0%*%g gsv <- CM0%*%v%*%t(CM0) sd <- matrix(sqrt(diag(gsv))) ci <- cbind(gs - 1.96*sd, gs + 1.96*sd) GradientShifts <- cbind(pstar, 100*(exp(gs)-1), 100*(exp(ci)-1)) dimnames(GradientShifts) <- list(c(), c("Period", "Mean Percent Change per Year of Age", "CILo", "CIHi")) # Wald test: Do all Gradient Shifts equal the CAT? if (PVP$HiC==1){ dgs <- -KPC%*%cd dgsv <- KPC%*%cdv%*%t(KPC) # If A==P the test has P - 1 df, otherwise the test has P df. if (A == P) { X214 <- t(dgs[1:(P-1)])%*%solve(dgsv[1:(P-1), 1:(P-1)], dgs[1:(P-1)]) df14 <- P - 1 } else { X214 <- t(dgs)%*%solve(dgsv, dgs) df14 <- P } PVAL14 <- pchisq(X214, df14, lower.tail = FALSE) } else{ X214 <- X28 df14 <- df8 PVAL14 <- PVAL8 } # # (12d) Quadratic Gradient Shifts # XCQp <- as.matrix(bdiag(list(XCQ, 1, 1))) g2 <- XCQp%*%B[c(6, 2, 3)] g2v <- XCQp%*%s2VAR[c(6, 2, 3), c(6, 2, 3)]%*%t(XCQp) gs2 <- CM0%*%g2 gs2v <- CM0%*%g2v%*%t(CM0) sd <- matrix(sqrt(diag(gs2v))) ci <- cbind(gs2 - 1.96*sd, gs + 1.96*sd) QuadGradientShifts <- cbind(pstar, 100*(exp(gs2)-1), 100*(exp(ci)-1)) dimnames(QuadGradientShifts) <- list(c(), c("Period", "Mean Percent Change per Year of Age", "CILo", "CIHi")) WaldTests <- matrix( c(X21, df1, PVAL1, X22, df2, PVAL2, X23, df3, PVAL3, X24, df4, PVAL4, X25, df5, PVAL5, X26, df6, PVAL6, X27, df7, PVAL7, X28, df8, PVAL8, X29, df9, PVAL9, X210, df10, PVAL10, X211, df11, PVAL11, X212, df12, PVAL12, X213, df13, PVAL13, X214, df14, PVAL14), 14, 3, byrow = TRUE) dimnames(WaldTests) <- list( c("NetDrift = 0", "THETAa = 0", "All Higher-Order Age Deviations = 0", "All Age Deviations = 0", "THETAp = 0", "All Higher-Order Period Deviations = 0", "All Period Deviations = 0", "THETAc = 0", "All Higher-Order Cohort Deviations = 0", "All Cohort Deviations = 0", "All Period RR = 1", "All Cohort RR = 1", "All Local Drifts = Net Drift", "All Gradient Shifts = CAT"), c("X2", "df", "P-Value")) CombinationTests <- matrix(NaN, nrow = 4) CombinationTests[1] <- min(c(2*PVAL5, 2*PVAL6, 1)) CombinationTests[2] <- min(c(3*PVAL1, 3*PVAL5, 3*PVAL6, 1)) CombinationTests[3] <- min(c(2*PVAL8, 2*PVAL9, 1)) CombinationTests[4] <- min(c(3*PVAL1, 3*PVAL8, 3*PVAL9, 1)) dimnames(CombinationTests) <- list(c("All Period Deviations = 0", "All PRR = 1 <=> FTT = constant", "All Cohort Deviations = 0", "All CRR = 1 <=> FCP = constant"), c()) # Matrix that converts from C-2 cohort parameters to C cohort deviations. # (For use outside this function). nc <- D$nc vC <- (1/(cstar[C]-cstar[1]))*( cstar[1]*matrix(nc[2:(C-1)]) - matrix(cstar[2:(C-1)]*nc[2:(C-1)])) v1 <- matrix(-nc[2:(C-1)]) - vC H <- rbind(t(v1), diag(C - 2), t(vC)) Variances <- list( ad = adv , ad2 = ad2v , adh = adhv , pd = pdv , pd2 = pd2v , pdh = pdhv , cd = cdv , cd2 = cd2v , cdh = cdhv , lac = lacv , lac2 = lac2v , lacrr = lacrrv , lac2rr = lac2rrv , xac = xacv , xac2 = xac2v , xacrr = xacrrv , xac2rr = xac2rrv , lvcrr = lvcrrv , lac_grad = lac_gradv , lac2_grad = lac2_gradv, lach_grad = lach_gradv, xac_grad = xac_gradv , xac2_grad = xac2_gradv, xach_grad = lach_gradv, ftt = fttv , ftt2 = ftt2v , prr = prrv , prr2 = prr2v , fcp = fcpv , fcp2 = fcp2v , crr = crrv , crr2 = crr2v , ld = ldv , ld2 = ld2v , def = defv , def2 = def2v , gs = gsv , gs2 = gs2v , pert = pertv , pert2 = pert2v ) Matrices <- list( X = D$X , XCO = XCO , XAQ = XAQ , XAD = XAD , XPQ = XPQ , XPD = XPD , XCQ = XCQ , XCD = XCD , XLA = XLA , XLAQ = XLAQ, XXA = XXA , XXAQ = XXAQ, XPT = XPT , XPTQ = XPTQ, XPR = XPR , XCR = XCR , XCB = XCB , CM = CM , XCT = XCT , XCTQ = XCTQ, XLX = XLX , LAR = LAR , KAC = KAC , KPC = KPC , GA = G , H = H ) M <- list( Inputs = PVP , FittedRates = FittedRates , Coefficients = Coefficients , NetDrift = NetDrift , GlobalCurvature = GlobalCurvature , QuadAgeDeviations = QuadAgeDeviations , HiOrdAgeDeviations = HiOrdAgeDeviations , AgeDeviations = AgeDeviations , QuadPerDeviations = QuadPerDeviations , HiOrdPerDeviations = HiOrdPerDeviations , PerDeviations = PerDeviations , QuadCohDeviations = QuadCohDeviations , HiOrdCohDeviations = HiOrdCohDeviations , CohDeviations = CohDeviations , LongAge = LongAge , LongAgeRR = LongAgeRR , QuadLongAge = QuadLongAge , QuadLongAgeRR = QuadLongAgeRR , CrossAge = CrossAge , CrossAgeRR = CrossAgeRR , QuadCrossAge = QuadCrossAge , QuadCrossAgeRR = QuadCrossAgeRR , Long2CrossRR = Long2CrossRR , QuadLongAgeGrad = QuadLongAgeGrad , HiOrdLongAgeGrad = HiOrdLongAgeGrad , LongAgeGrad = LongAgeGrad , QuadCrossAgeGrad = QuadCrossAgeGrad , HiOrdCrossAgeGrad = HiOrdCrossAgeGrad , CrossAgeGrad = CrossAgeGrad , QuadFittedTemporalTrends = QuadFittedTemporalTrends, QuadPeriodRR = QuadPeriodRR , FittedTemporalTrends = FittedTemporalTrends , PeriodRR = PeriodRR , FittedCohortPattern = FittedCohortPattern , CohortRR = CohortRR , QuadFittedCohortPattern = QuadFittedCohortPattern , QuadCohortRR = QuadCohortRR , LocalDrifts = LocalDrifts , QuadLocalDrifts = QuadLocalDrifts , Deflections = Deflections , QuadDeflections = QuadDeflections , Perturbations = Perturbations , QuadPerturbations = QuadPerturbations , GradientShifts = GradientShifts , QuadGradientShifts = QuadGradientShifts , WaldTests = WaldTests , CombinationTests = CombinationTests , Variances = Variances , Matrices = Matrices , APCModel = apcM , Pt = Pt , nc = nc ) M } checkPVPPAIRS <- function(R, OverDispersion = 1, offset_tick = 10^5, zero_fill = 0.1, RVals = c(NaN, NaN, NaN), HiC = TRUE, HiP = TRUE, HiA = TRUE) { # Pre-process input data and parameters D <- rates2data_set(R) if (all(is.nan(RVals))) { # Calculate default reference values A <- length(D$a) P <- length(D$p) aref <- floor((A+1)/2) pref <- floor((P+1)/2) cref <- pref - aref + A ageref <- D$a[aref] perref <- D$p[pref] cohref <- D$c[cref] RVals <- c(ageref, perref, cohref) } else { # Valdidate user-supplied reference values RVals <- floor(RVals) TA <- is.element(RVals[1], D$a) if (!TA) { RVals[1] <- floor(RVals[1]) + 0.5 TA <- is.element(RVals[1], D$a) } TB <- is.element(RVals[2], D$p) if (!TB) { RVals[2] <- floor(RVals[2]) + 0.5 TB <- is.element(RVals[2], D$p) } TC <- is.element(RVals[3], D$c) if (!TC) { RVals[3] <- floor(RVals[3]) + 0.5 TC <- is.element(RVals[3], D$c) } if (!(TA && TB && TC)) stop("Invalid Age, Period, or Cohort reference value.") end } # Replace 0 events with zero_fill value. e <- matrix(D$DATA[,4]) e[e==0] <- zero_fill D$DATA[,4] <- e PVP <- list(D = D, RVals = RVals, OverDispersion = OverDispersion, offset_tick = offset_tick, zero_fill = zero_fill, HiC = HiC, HiP = HiP, HiA = HiA) PVP } APCFIT = function(PVP, X) { # Fit APC model by weighted least squares A <- length(PVP$D$a) P <- length(PVP$D$p) n <- nrow(X) p <- ncol(X) p0 = ncol(X) Y <- PVP$D$DATA[,4:5]; offset <- matrix(log(Y[,2])) y <- matrix(Y[,1]) ly <- log(y) W <- y WX <- (W%*%matrix(1,ncol=p0))*X z <- (ly - offset) XTWX <- t(X)%*%WX B <- solve(t(X)%*%WX,t(WX)%*%z) V <- solve(t(X)%*%WX) u <- matrix(Y[,2]*exp(X%*%B)) wr2 <- matrix(W*(z-X%*%B)^2) DEVRESIDS <- sign(y-u)*sqrt(wr2) DEV <- sum(wr2) if (PVP$OverDispersion==1){ s2 <- max(1, DEV/(n-p0)) } else {s2 <- 1} s2V = s2*V APCMODEL <- list(B = B, s2 = s2, s2VAR = s2V, DEV = DEV, DevResids = DEVRESIDS) APCMODEL } designmatrix <- function(PVP) { # Calculate design matrix for APC model N <- nrow(PVP$D$DATA) J <- matrix(1, nrow = N) I.N <- diag(N) astar <- matrix(PVP$D$DATA[,1]) astarbar <- mean(astar) astar_0 <- astar - astarbar A <- length(PVP$D$a) pstar <- matrix(PVP$D$DATA[,2]) pstarbar <- mean(pstar) pstar_0 <- pstar - pstarbar P <- length(PVP$D$p) cstar <- matrix(PVP$D$DATA[,3]) cstarbar <- pstarbar - astarbar cstar_0 <- cstar - cstarbar C <- length(PVP$D$c) # Age Ad <- kronecker(matrix(1, nrow=P), diag(A)) # Per Pd <- kronecker(diag(P), matrix(1, nrow=A)) # Coh Cd <- matrix(NaN,nrow=N,ncol=C) for (i in 1:C) Cd[,i] <- cstar == PVP$D$c[i] end nc <- t( diag( t(Cd)%*%Cd ) ) # Projection matrix orthogonal to intercept, linear age, and quadratic age Delta <- astar[2] - astar[1] a00 <- astar[1] - 0.5*Delta if (A>=3) { qa2 <- astar^2 - (Delta*A + 2*a00)*astar + (a00 + Delta*A/2)^2 - ((Delta^2)/12)*(A-1)*(A+1) Xa2 <- cbind(J, astar_0, qa2) xtxi <- (1/P)*diag( 1/ c(A, (Delta^2/12)*(A-1)*A*(A+1), (Delta^4/180)*(A-2)*(A-1)*A*(A+1)*(A+2) ) ) Ra2 <- xtxi%*%t(Xa2) PA12 <- Xa2%*%Ra2 XAD12 <- I.N - PA12 } else { qa2 <- matrix(0, N) Xa2 <- cbind(J, astar_0) xtxi <- (1/P)*diag(1/ c(A, (Delta^2/12)*(A-1)*A*(A+1)) ) Ra2 <- xtxi%*%t(Xa2) PA12 <- Xa2%*%Ra2 XAD12 <- I.N - PA12 } # Projection matrix orthogonal to intercept, linear period, and quadratic period p00 <- pstar[1] - 0.5*Delta if (P >= 3) { qp2 <- pstar^2 - (Delta*P + 2*p00)*pstar + (p00 + Delta*P/2)^2 - ((Delta^2)/12)*(P-1)*(P+1) Xp2 <- cbind(J, pstar_0, qp2) xtxi <- (1/A)*diag( 1/ c(P, (Delta^2/12)*(P-1)*P*(P+1), (Delta^4/180)*(P-2)*(P-1)*P*(P+1)*(P+2) ) ) Rp2 <- xtxi%*%t(Xp2) PP12 <- Xp2%*%Rp2 XPD12 <- I.N - PP12 } else { qp2 <- pstar^2 - (Delta*P + 2*p00)*pstar + (p00 + Delta*P/2)^2 - ((Delta^2)/12)*(P-1)*(P+1) Xp2 <- cbind(J, pstar_0) xtxi <- (1/A)*diag( 1/ c(P, (Delta^2/12)*(P-1)*P*(P+1) ) ) Rp2 <- xtxi%*%t(Xp2) PP12 <- Xp2%*%Rp2 XPD12 <- I.N - PP12 } c00 <- p00 - a00 qc2 <- cstar^2 - (Delta*(A+P) + 2*(c00 - Delta*A))*cstar + (c00 - Delta*A)^2 + (A+P)*Delta*cstarbar - ((Delta^2)/6)*(2*A^2 + 3*A*P + 2*P^2 - 1) Xc2 <- cbind(J, cstar_0, qc2) xtxi <- (1/(A*P))*diag( 1/ c(1, (Delta^2/12)*(A^2 + P^2 - 2), (Delta^4/180)*( (A^2 + P^2)^2 - 10*(A^2 + P^2) + 3*A^2*P^2 + 13 ) ) ) Rc2 <- xtxi%*%t(Xc2) PC12 <- Xc2%*%Rc2 XCD12 <- I.N - PC12 # Orthogonalize the deviations Ad0 <- XAD12%*%Ad Pd0 <- XPD12%*%Pd Cd0 <- XCD12%*%Cd if (A>3) { aCOL <- 2:(A-2) } else { # Ad0 column space is empty aCOL <- NULL } if (P>3){ pCOL <- 2:(P-2) } else { # Pd0 column space is empty pCOL <- NULL } X <- cbind(J, astar_0, cstar_0, qa2, qp2, qc2, Ad0[,aCOL], Pd0[,pCOL], Cd0[,2:(C-2)]) # Pointers for Higher-Order Deviations: # # * For age: # 1:(A - 3) shift forward by 6 -> (1 + 6):(A - 3 + 6) = 7:(A + 3) # * For period: # 1:(P - 3) shift forward by 6 + (A - 3) -> # (1 + 6 + A - 3):(P - 3 + 6 + A - 3) = (A + 4):(A + P) # * For cohort: # 1:(C - 3) shift forward by 6 + (A - 3) + (P - 3) -> # (1 + 6 + (A - 3) + (P - 3)):(C - 3) + 6 + (A - 3) + (P - 3) = # (A + P + 1):(A + P + C - 3) Pt <- vector("list", 9) Pt[[1]] <- 1 Pt[[2]] <- 2 Pt[[3]] <- 3 Pt[[4]] <- 4 Pt[[5]] <- 5 Pt[[6]] <- 6 Pt[[7]] <- 7:(A+3) Pt[[8]] <- (A+4):(A+P) Pt[[9]] <- (A+P+1):(A+P+C-3) if (A==3) { # 7:(A+3) runs backwards but (A+4) = 7 # need to reset Pt[[7]] to NaN but Pt{8} & Pt{9} are correct Pt[[7]] <- NaN } if (P==3) { Pt[[8]] <- NaN } # Pointer for cohort effects, (last f/u, South+East edges of Lexis) INCc <- c(A*(1:P), P*A-(1:(A-1))) # Pointer for age effects, (South edge of Lexis) INCa <- ((P-1)*A+1):(P*A) # Pointer for period effects (East edge of Lexis) INCp <-A*(1:P) D <- list(X = X, INCa = INCa, INCp = INCp, INCc = INCc, Pt = Pt, nc = nc) D } rates <- function(EVENTS, OFFSET, fullname = character(0), label = character(0), description = character(0), event_label = "events", offset_units = "offset units", offset_tick = 100000, api = c(NaN, NaN, NaN), ages = NaN, periods = NaN) { ### # Validate EVENTS and OFFSET ### if (is.matrix(EVENTS) && all(EVENTS >= 0)) { A <- nrow(EVENTS) P <- ncol(EVENTS) } else { stop("events must be a matrix of non-negative values.") } if (!(is.matrix(OFFSET) && all(OFFSET >= 0) && nrow(OFFSET)==A && ncol(OFFSET)==P)) { stop("offset must be a matrix of non-negative values the same size as events.") } ### # Validate ages and periods specified by api or ages/periods combination. ### if (!all(is.nan(api))) { if ( !( length(api==3) && is.numeric(api) ) ) {stop("api must be a 3 element vector.")} if (!all(api==round(api))) {stop("api values must be single-years.")} a0 <- api[1] p0 <- api[2] INTERVAL <- api[3] ages <- a0 + INTERVAL*(0:A) periods <- p0 + INTERVAL*(0:P) a <- ages[1:A-1] + 0.5*INTERVAL p <- periods[1:P-1] + 0.5*INTERVAL } else { if (! ( !any(is.nan(ages)) && !any(is.nan(periods)) && all(ages==round(ages)) && all(periods==round(periods)) && (length(ages)==A+1) && (length(periods)==P+1) ) ) {stop("Invalid cutpoints for ages or periods.")} } ### # Validate text inputs. ### if (!is.character(fullname)) {stop('fullname must be a character string.')} else {fullname <- gsub("^\\s+|\\s+$", '', fullname)} if (length(fullname)==0) {fullname <- paste(c(toString(A), '-by-', toString(P), ' rates object'), collapse = "")} if (length(description)==0) {description <- paste(c('Created ', date()), collapse = "")} if (!is.character(event_label)) {stop('event_label must be a character string.')} if (!is.character(offset_units)) {stop('offset_units must be a character string.')} if (length(label)==0) {label <- fullname} #aL <- as.character(seq(from = ages[1], to = ages[A], by = 2)) #aH <- as.character(seq(from = ages[2], to = ages[A+1], by = 2)) da <- ages[2:(A+1)] - ages[1:A] if (!all(da==1)) { aL <- as.character(seq(from = ages[1], to = ages[A], by = 2)) aH <- as.character(seq(from = ages[2], to = ages[A+1], by = 2)) age_labels <- paste(cbind(aL), ' - ', cbind(aH))} else {aL <- as.character(seq(from = ages[1], to = ages[A], by = 1)) age_labels <- aL} #pL <- as.character(seq(from = periods[1], to = periods[P], by = 2)) #pH <- as.character(seq(from = periods[2], to = periods[P+1], by = 2)) dp <- periods[2:(P+1)] - periods[1:P] if (!all(dp==1)) { pL <- as.character(seq(from = periods[1], to = periods[P], by = 2)) pH <- as.character(seq(from = periods[2], to = periods[P+1], by = 2)) per_labels <- paste(cbind(pL), ' - ', cbind(pH))} else { pL <- as.character(seq(from = periods[1], to = periods[P], by = 1)) per_labels <- pL} a <- ages[1:A] + da/2 p <- periods[1:P] + dp/2 R <- list(fullname = fullname, label = label, description = description, events = EVENTS, event_label = event_label, offset = OFFSET, offset_units = offset_units, offset_tick = offset_tick, ages = ages, age_labels = age_labels, a = a, periods = periods, per_labels = per_labels, p = p) R } rates2data_set <- function(R) { A <- nrow(R$events) P <- ncol(R$events) da <- R$ages[2:(A+1)] - R$ages[1:A] D.a <- R$ages[1:A] + 0.5*da dp <- R$periods[2:(P+1)] - R$periods[1:P] D.p <- R$periods[1:P] + 0.5*dp ADATA <- kronecker(matrix(1, nrow=P), matrix(D.a, nrow=A)) PDATA <- kronecker(matrix(D.p, nrow=P), matrix(1, nrow=A)) CDATA <- PDATA - ADATA D.c <- sort(c(unique(CDATA))) E <- c(R$events) O <- c(R$offset) D.DATA = cbind(ADATA, PDATA, CDATA, E, O) colnames(D.DATA)<-c("Age","Period","Cohort","Events","Offset") D <- list(name = R$name, description = R$description, DATA = D.DATA, a = D.a, p = D.p, c = D.c) D } rates2csv <- function(R, FILE = "rates.csv") { A <- nrow(R$events) P <- ncol(R$events) comma = ',' commas <- paste(replicate(P-1, ','), collapse="") # fprintf(fid, ['Title: ' R.fullname commas '\n']); s1 <- paste(c('Title: ', R$fullname, commas), collapse = "") s2 <- paste(c('Description: ', R$description, commas), collapse = "") s3 <- paste(c('Start Year: ', toString(R$p[1]), commas), collapse = "") s4 <- paste(c('Start Age: ', toString(R$a[1]), commas), collapse = "") s5 <- paste(c('Interval (Years): ', toString(R$a[2]-R$a[1]), commas), collapse = "") DATA <- matrix(NaN, A, 2*P) DATA[,seq(1, 2*P-1, by = 2)] <- R$events DATA[,seq(2, 2*P, by = 2)] <- R$offset cat(s1, s2, s3, s4, s5, file = FILE, sep = "\n", append = FALSE) write(t(DATA), file = FILE, append = TRUE, sep = ',', ncolumns = 2*P) Fout <- list(s1 = s1, s2 = s2, s3 = s3, s3 = s4, s5 = s5, DATA = DATA) Fout } simple_csv2rates <- function(FILE,StartYear,StartAge,Interval,fullname,description) { # DATA is a data.frame DATA <- read.table(FILE, header = FALSE, sep = ',') PP <- ncol(DATA) A <- nrow(DATA) E = as.matrix(DATA[, seq(1, PP-1, by = 2)]) dimnames(E) <- NULL O = as.matrix(DATA[, seq(2, PP, by = 2)]) dimnames(O) <- NULL a <- seq(from = StartAge, by = Interval, to = StartAge + Interval*A) p <- seq(from = StartYear, by = Interval, to = StartYear + Interval*PP/2) R <- rates(E, O, fullname = fullname, description = description, ages = a, periods = p) R } prepare_rates <- function(DATA,StartYear,StartAge,Interval,fullname,description) { # DATA is a data.frame PP <- ncol(DATA) A <- nrow(DATA) E = as.matrix(DATA[, seq(1, PP-1, by = 2)]) dimnames(E) <- NULL O = as.matrix(DATA[, seq(2, PP, by = 2)]) dimnames(O) <- NULL a <- seq(from = StartAge, by = Interval, to = StartAge + Interval*A) p <- seq(from = StartYear, by = Interval, to = StartYear + Interval*PP/2) R <- rates(E, O, fullname = fullname, description = description, ages = a, periods = p) R } type <- function(R, comp = 'r') { A <- nrow(R$events) P <- ncol(R$events) if (comp == "r") { Tout <- R$offset_tick*R$events/R$offset dn <- paste('Rates - ', R$fullname) } else if (comp == "e") { Tout <- R$events dn <- paste("Events -", R$fullname) } else if (comp == "o") { Tout <- R$offset dn <- paste('Offset - ', R$fullname) } else if (comp == "eo") { DATA <- matrix(NaN, nrow = A, ncol = 2*P) DATA[,seq.int(1, 2*P, 2)]<-R$events DATA[,seq.int(2, 2*P, 2)]<-R$offset Tout <- DATA dn <- paste("Events & Offset -", R$fullname) } else if (comp == "er") { r <- R$offset_tick*R$events/R$offset DATA <- matrix(NaN, nrow = A, ncol = 2*P) DATA[,seq.int(1, 2*P, 2)]<-R$events DATA[,seq.int(2, 2*P, 2)]<-r Tout <- DATA dn <- paste("Events and Rates - ", R$fullname) } else if (comp == "eor") { r <- R$offset_tick*R$events/R$offset DATA <- matrix(NaN, nrow = A, ncol = 3*P) DATA[,seq.int(1, 3*P, 3)]<-R$events DATA[,seq.int(2, 3*P, 3)]<-R$offset DATA[,seq.int(3, 3*P, 3)]<-r Tout <- DATA dn <- paste("Events, offset, and Rates - ", R$fullname) } else if (comp == "rci") { r <- R$offset_tick*R$events/R$offset v <- (R$offset_tick^2)*R$events/R$offset^2 cilo <- r - 1.96*sqrt(v) cilo[cilo<0] <- 0 cihi <- r + 1.96*sqrt(v) DATA <- matrix(NaN, nrow = A, ncol = 3*P) DATA[,seq.int(1, 3*P, 3)]<-r DATA[,seq.int(2, 3*P, 3)]<-cilo DATA[,seq.int(3, 3*P, 3)]<-cihi Tout <- DATA dn <- paste("Rates and 95% CI - ", R$fullname) } else { } Tout <- list(name = dn, DATA = Tout, ages = R$ages, periods = R$periods) Tout } csv2rates <- function(FILE) { StartYear <- 0 StartAge <- 0 Interval <- 1 header <- scan(FILE, nlines = 5, what = character(0), sep = '/', quiet = 1) # Strip any excess delimeters header <- gsub(",", "", header) # Strip leading and trailing white space header <- gsub("^\\s+|\\s+$", "", header) H = length(header) k <- 0 for (h in 1:H) { headerh = header[h] nc = nchar(headerh) f <- regexpr("Title:", headerh, ignore.case = TRUE) d <- regexpr("Description:", headerh, ignore.case = TRUE) p <- regexpr("Start Year:", headerh, ignore.case = TRUE) a <- regexpr("Start Age:", headerh, ignore.case = TRUE) i <- regexpr("Interval \\(Years\\):", headerh, ignore.case = TRUE) if (f==1) { fullname <- gsub("^\\s+|\\s+$", "", substr(headerh, attr(f, "match.length")+1, nc)) k <- k + 1 } if (d==1) { description <- gsub("^\\s+|\\s+$", "", substr(headerh, attr(d, "match.length")+1,nc)) k <- k + 1 } if (p==1) { StartYear <- as.numeric(substr(headerh, attr(p, "match.length")+1,nc)) k <- k + 1 } if (a==1) { StartAge <- as.numeric(substr(headerh, attr(a, "match.length")+1,nc)) k <- k + 1 } if (i==1) { Interval <- as.numeric(substr(headerh, attr(i, "match.length")+1,nc)) k <- k + 1 } } # DATA is a data.frame DATA <- read.table(FILE, skip = k, header = FALSE, sep = ',') PP <- ncol(DATA) A <- nrow(DATA) E = as.matrix(DATA[, seq(1, PP-1, by = 2)]) dimnames(E) <- NULL O = as.matrix(DATA[, seq(2, PP, by = 2)]) dimnames(O) <- NULL a <- seq(from = StartAge, by = Interval, to = StartAge + Interval*A) p <- seq(from = StartYear, by = Interval, to = StartYear + Interval*PP/2) R <- rates(E, O, fullname = fullname, description = description, ages = a, periods = p) R } plot.apc1 <- function(M) { # Plot first set of age-period-cohort estimable functions. par(mfrow = c(4,3)) DATA <- cbind(matrix(M$LongAge[,1]), (M$LongAge[,c(2,3,4)])) dimnames(DATA) <- list(c(), c("Age", "Rate", "CILo", "CIHi")) pcurve(DATA, col = "darkred", colf = "pink", lwd = 3, cex = 1.0, pch = 21) title(main = "Longitudinal Age Curve", cex.main = 1) DATA <- cbind(matrix(M$CrossAge[,1]), (M$CrossAge[,c(2,3,4)])) dimnames(DATA) <- list(c(), c("Age", "Rate", "CILo", "CIHi")) pcurve(DATA, lwd = 3, col = "darkred", colf = "pink", cex = 1.0, pch = 21) title(main = "Cross-Sectional Age Curve", cex.main = 1) pcurve(M$Long2CrossRR, lwd = 3, col = "darkred", colf = "pink", cex = 1.0, pch = 21) abline(1,0, lty = 3) title(main = "Long vs. Cross RR", cex.main = 1) pcurve(M$FittedTemporalTrends, col = "steelblue4", colf = "slategray1", lwd = 3, cex = 1.0, pch = 21) title(main = "Fitted Temporal Trends", cex.main = 1) pcurve(M$PeriodRR, col = "steelblue4", colf = "slategray1", lwd = 3, cex = 1.0, pch = 21) abline(1, 0, lty = 3) title(main = "Period RR", cex.main = 1) pcurve(M$CohortRR, col = "seagreen4", colf = "darkseagreen1", lwd = 3, cex = 1.0, pch = 21) abline(1, 0, lty = 3) title(main = "Cohort RR", cex.main = 1) pcurve(M$LocalDrifts, col = "black", colf = "grey88", lwd = 3, cex = 1.0, pch = 21) abline(0, 0, lty = 3) title(main = "Local Drifts", cex.main = 1) pcurve(M$AgeDeviations, col = "darkred", colf = "pink", lwd = 3, cex = 1.0, pch = 21) abline(0, 0, lty = 3) title(main = "Age Deviations", cex.main = 1) pcurve(M$PerDeviations, col = "steelblue4", colf = "slategray1", lwd = 3 , cex = 1.0, pch = 21) abline(0, 0, lty = 3) title(main = "Period Deviations", cex.main = 1) pcurve(M$CohDeviations, col = "seagreen4", colf = "darkseagreen1", lwd = 3, cex = 1.0, pch = 21) abline(0, 0, lty = 3) title(main = "Cohort Deviations", cex.main = 1) pcurve(M$FittedCohortPattern, col = "seagreen4", colf = "darkseagreen1", lwd = 3, cex = 1.0, pch = 21) abline(0, 0, lty = 3) title(main = "Fitted Cohort Pattern", cex.main = 1) } plot.apc2 <- function(M) { # Plot second set of age-period-cohort estimable functions. par(mfrow = c(4,3)) DATA <- cbind(matrix(M$QuadLongAge[,1]), (M$QuadLongAge[,c(2,3,4)])) dimnames(DATA) <- list(c(), c("Age", "Rate", "CILo", "CIHi")) pcurve(DATA, col = "darkred", colf = "pink", lwd = 3, cex = 1.0, pch = 21) title(main = "Longitudinal Age Curve", cex.main = 1) DATA <- cbind(matrix(M$QuadCrossAge[,1]), (M$QuadCrossAge[,c(2,3,4)])) dimnames(DATA) <- list(c(), c("Age", "Rate", "CILo", "CIHi")) pcurve(DATA, lwd = 3, col = "darkred", colf = "pink", cex = 1.0, pch = 21) title(main = "Cross-Sectional Age Curve", cex.main = 1) pcurve(M$Long2CrossRR, lwd = 3, col = "darkred", colf = "pink", cex = 1.0, pch = 21) abline(1,0, lty = 3) title(main = "Long vs. Cross RR", cex.main = 1) pcurve(M$QuadFittedTemporalTrends, col = "steelblue4", colf = "slategray1", lwd = 3, cex = 1.0, pch = 21) title(main = "Fitted Temporal Trends", cex.main = 1) pcurve(M$QuadPeriodRR, col = "steelblue4", colf = "slategray1", lwd = 3, cex = 1.0, pch = 21) abline(1, 0, lty = 3) title(main = "Period RR", cex.main = 1) pcurve(M$QuadCohortRR, col = "seagreen4", colf = "darkseagreen1", lwd = 3, cex = 1.0, pch = 21) abline(1, 0, lty = 3) title(main = "Cohort RR", cex.main = 1) pcurve(M$QuadLocalDrifts, col = "black", colf = "grey88", lwd = 3, cex = 1.0, pch = 21) abline(M$NetDrift[1,1], 0, lty = 3) title(main = "Local Drifts", cex.main = 1) pcurve(M$QuadAgeDeviations, col = "darkred", colf = "pink", lwd = 3, cex = 1.0, pch = 21) title(main = "Age Deviations", cex.main = 1) pcurve(M$QuadPerDeviations, col = "steelblue4", colf = "slategray1", lwd = 3 , cex = 1.0, pch = 21) abline(0, 0, lty = 3) title(main = "Period Deviations", cex.main = 1) pcurve(M$QuadCohDeviations, col = "seagreen4", colf = "darkseagreen1", lwd = 3, cex = 1.0, pch = 21) abline(0, 0, lty = 3) title(main = "Cohort Deviations", cex.main = 1) pcurve(M$QuadFittedCohortPattern, col = "seagreen4", colf = "darkseagreen1", lwd = 3, cex = 1.0, pch = 21) abline(0, 0, lty = 3) title(main = "Fitted Cohort Pattern", cex.main = 1) DATA <- cbind(matrix(M$QuadGradientShifts[,1]), (M$QuadGradientShifts[,c(2,3,4)])) dimnames(DATA) <- list(c(), c("Period", "Mean % Change Per Year of Age", "CILo", "CIHi")) pcurve(DATA, col = "darkred", colf = "pink", lwd = 3, cex = 1.0, pch = 21) title(main = "Gradient Shifts", cex.main = 1) abline(100*(exp(M$Coefficients[4,1]) - 1), 0, lty = 3) } pcurve <- function(DATA, col = 'steelblue4', colf = 'slategray1', bg = 'grey99', pch = 1, type = 'b', lty = 1, lwd = 2, XLim = NA, YLim = NA, cex = 1.5) { # Plot a selected output from an age-period-cohort model. x <- DATA[,1] rangex <- range(x)[2]-range(x)[1] if (is.na(XLim[1])) {XLim <- c(min(x)-0.05*rangex, max(x)+0.05*rangex)} xl <- dimnames(DATA)[[2]][1] y <- DATA[,2] rangey <- range(y)[2]-range(y)[1] if (is.na(YLim[1])) {YLim <- c(min(y)-0.05*rangey, max(y)+0.05*rangey)} yl <- dimnames(DATA)[[2]][2] if (ncol(DATA)==4){ xci <- c(x, rev(x)) yci <- c(DATA[, 3], rev(DATA[, 4])) rangey <- range(yci)[2]-range(yci)[1] YLim <- c(min(yci) - 0.05*rangey, max(yci) + 0.05*rangey)} else {yci <- NULL} plot(x, y, col = col, pch = pch, type = type, lty = lty, lwd = lwd, cex = cex, xlab = xl, ylab = yl, xlim = XLim, ylim = YLim, las = 1, bg = bg) if (!is.null(yci[1])){ polygon(xci, yci, col = colf, border = colf) } points(x, y, col = col, pch = pch, type = type, lty = lty, lwd = lwd, cex = cex, xlab = xl, ylab = yl, xlim = XLim, ylim = YLim, las = 1, bg = bg) }