# ----------------------------------------------------------------------------- # DPhil Chapter 2 # combined.functions.R # Zoe Fannon # Date started: 05/07/18 # Last edited: 24/08/18 # ----------------------------------------------------------------------------- # consolidate into three functions which make calls to internal functions # ----------------------------------------------------------------------------- #### write apc.indiv.est.model #### apc.indiv.est.model <- function(data, unit=1, n.coh.excl.start=0, n.coh.excl.end=0, n.per.excl.start=0, n.per.excl.end=0, n.age.excl.start=0, n.age.excl.end=0, model.design="APC", dep.var=NULL, covariates=NULL, model.family=NULL, NR.controls=NULL, existing.collinear = NULL, existing.design = NULL){ # Check that design and family are correctly specified design.list <- c("APC", "AP", "AC", "PC", "Ad", "Pd", "Cd", "A", "P", "C", "t", "tA", "tP", "tC", "1") if (!isTRUE(model.design %in% c(design.list, "TS"))) stop("model.design not recognised") family.list <- c("gaussian", "binomial") if (!isTRUE(model.family %in% c(family.list))) stop("model.family not recognised") # Estimate any of the APC model or submodels if (model.design %in% design.list){ # First get collinear design matrix collinear <- existing.collinear if (is.null(collinear)){ collinear <- apc.indiv.design.collinear(data=data, unit=unit, n.coh.excl.start = n.coh.excl.start, n.coh.excl.end = n.coh.excl.end, n.per.excl.start = n.per.excl.start, n.per.excl.end = n.per.excl.end, n.age.excl.start = n.age.excl.start, n.age.excl.end = n.age.excl.end) } # Then get model-specific design matrix design <- existing.design if (is.null(design)){ design <- apc.indiv.design.model(collinear, model.design = model.design, dep.var = dep.var, covariates = covariates) } # Finally estimate the model model <- apc.indiv.fit.model(design, model.family = model.family) } else if (model.design=="TS"){ # Estimate the time-saturated model if (model.family == "gaussian"){ # The gaussian model is estimated analytically model <- apc.indiv.estimate.TS(data=data, dep.var = dep.var, covariates = covariates) } else if (model.family == "binomial"){ # The binomial model requires Newton-Rhapson approximation # Set defaults if no user-specified Newton-Rhapson parameters if (is.null(NR.controls$maxit.loop)) NR.controls$maxit.loop <- 10 if (is.null(NR.controls$maxit.linesearch)) NR.controls$maxit.linesearch <- 30 if (is.null(NR.controls$tolerance)) NR.controls$tolerance <- .002 if (is.null(NR.controls$init)) NR.controls$init <- "ols" if (is.null(NR.controls$d1.tol)) NR.controls$d1.tol <- .002 # Estimate the TS model by Newton-Rhapson model <- apc.indiv.logit.TS(data = data, dep.var = dep.var, covariates = covariates, maxit.loop = NR.controls$maxit.loop, maxit.linesearch = NR.controls$maxit.linesearch, tolerance = NR.controls$tolerance, init = NR.controls$init, inv.tol = NR.controls$inv.tol, d1.tol = NR.controls$d1.tol, custom.kappa = NR.controls$custom.kappa, custom.zeta = NR.controls$custom.zeta) } } return(model) } #### write apc.indiv.model.table #### apc.indiv.model.table <- function(data, dep.var, covariates = NULL, unit=1, n.coh.excl.start = 0, n.coh.excl.end = 0, n.age.excl.start = 0, n.age.excl.end = 0, n.per.excl.start = 0, n.per.excl.end = 0, model.family=NULL, NR.controls = NULL){ # Check that family is correctly specified family.list <- c("gaussian", "binomial") if (!isTRUE(model.family %in% c(family.list))) stop("model.family not recognised") # For the gaussian family, use an Ftable if (model.family == "gaussian") table <- apc.indiv.ftable(data = data, dep.var = dep.var, unit = unit, covariates = covariates, n.coh.excl.start = n.coh.excl.start, n.coh.excl.end = n.coh.excl.end, n.age.excl.start = n.age.excl.start, n.age.excl.end = n.age.excl.end, n.per.excl.start = n.per.excl.start, n.per.excl.end = n.per.excl.end) # For the binomial family, use an LR table if (model.family == "binomial"){ # Specify defaults to control Newton-Rhapson of time-saturated model if (is.null(NR.controls$maxit.loop)) NR.controls$maxit.loop <- 10 if (is.null(NR.controls$maxit.linesearch)) NR.controls$maxit.linesearch <- 30 if (is.null(NR.controls$init)) NR.controls$init <- "ols" table <- apc.indiv.LRtable.TS(data = data, dep.var = dep.var, unit = unit, covariates = covariates, n.coh.excl.start = n.coh.excl.start, n.coh.excl.end = n.coh.excl.end, n.age.excl.start = n.age.excl.start, n.age.excl.end = n.age.excl.end, n.per.excl.start = n.per.excl.start, n.per.excl.end = n.per.excl.end, maxit.loop.TS = NR.controls$maxit.loop, maxit.linesearch.TS = NR.controls$maxit.linesearch, init.TS = NR.controls$init) } return(table) } #### write apc.indiv.compare.direct #### apc.indiv.compare.direct <- function(data, big.model, small.model, unit=1, dep.var, covariates=NULL, model.family, n.coh.excl.start=0, n.coh.excl.end=0, n.age.excl.start=0, n.age.excl.end=0, n.per.excl.start=0, n.per.excl.end=0, NR.controls=NULL){ # Check that family is correctly specified family.list <- c("gaussian", "binomial") if (!isTRUE(model.family %in% c(family.list))) stop("model.family not recognised") # Check that all models are correctly specified design.list <- c("APC", "AP", "AC", "PC", "Ad", "Pd", "Cd", "A", "P", "C", "t", "tA", "tP", "tC", "1") if (!isTRUE(big.model %in% c(design.list, "TS"))) stop("big.model not recognised") if (!isTRUE(small.model %in% design.list)) stop("small.model not recognised; TS cannot be small.model") # For the gaussian family, use an Ftest if (model.family == "gaussian"){ if (big.model %in% design.list){ test <- apc.indiv.ftest.fullapc(data = data, big.model = big.model, small.model = small.model, dep.var = dep.var, covariates = covariates, unit = unit, n.coh.excl.start = n.coh.excl.start, n.coh.excl.end = n.coh.excl.end, n.age.excl.start = n.age.excl.start, n.age.excl.end = n.age.excl.end, n.per.excl.start = n.per.excl.start, n.per.excl.end = n.per.excl.end) } else if (big.model == "TS"){ test <- apc.indiv.ftest.TS(data = data, small.model = small.model, dep.var = dep.var, covariates = covariates, unit = unit, n.coh.excl.start = n.coh.excl.start, n.coh.excl.end = n.coh.excl.end, n.age.excl.start = n.age.excl.start, n.age.excl.end = n.age.excl.end, n.per.excl.start = n.per.excl.start, n.per.excl.end = n.per.excl.end) } } # For the binomial family, use likelihood ratio tests if (model.family == "binomial"){ if (big.model %in% design.list){ test <- apc.indiv.LRtest.fullapc(data = data, big.model = big.model, small.model = small.model, dep.var = dep.var, covariates = covariates, unit = unit, n.coh.excl.start = n.coh.excl.start, n.coh.excl.end = n.coh.excl.end, n.age.excl.start = n.age.excl.start, n.age.excl.end = n.age.excl.end, n.per.excl.start = n.per.excl.start, n.per.excl.end = n.per.excl.end) } else if (big.model == "TS"){ # Specify defaults to control Newton-Rhapson of time-saturated model if (is.null(NR.controls$maxit.loop)) NR.controls$maxit.loop <- 10 if (is.null(NR.controls$maxit.linesearch)) NR.controls$maxit.linesearch <- 30 if (is.null(NR.controls$init)) NR.controls$init <- "ols" test <- apc.indiv.LRtest.TS(data = data, small.model = small.model, dep.var = dep.var, covariates = covariates, unit = unit, n.coh.excl.start = n.coh.excl.start, n.coh.excl.end = n.coh.excl.end, n.age.excl.start = n.age.excl.start, n.age.excl.end = n.age.excl.end, n.per.excl.start = n.per.excl.start, n.per.excl.end = n.per.excl.end, maxit.loop.TS = NR.controls$maxit.loop, maxit.linesearch.TS = NR.controls$maxit.linesearch, init.TS = NR.controls$init) } } return(test) } #### write apc.indiv.design.collinear ##### apc.indiv.design.collinear <- function(data, unit=1, n.coh.excl.start=0, n.coh.excl.end=0, n.per.excl.start=0, n.per.excl.end=0, n.age.excl.start=0, n.age.excl.end=0){ # --------------------------------------------------------------------------- # This function generates a collinear design matrix which encompasses all # possible APC submodels and covariate combinations. # --------------------------------------------------------------------------- ## Step 1: Checks #### # check data has appropriate variable names if(isTRUE(("age" %in% colnames(data)) && ("period" %in% colnames(data)) && ("cohort" %in% colnames(data))) == FALSE) stop("apc.error: data missing one of age, period, cohort") # add cell.name if (!"cell.name" %in% colnames(data)) { cell.name <- paste("ik", as.character(data$age), as.character(data$cohort), sep="_") data <- cbind(cell.name, data) } else { warning("Variable cell.name must be unique identifier of age-cohort cells; if preexisting variable cell.name has alternate meaning please rename") } ## Step 2: Get data dimensions and mu.index #### # get data minimum values # first two used in construction of i.value, k.value age1 <- min(data$age ) coh1 <- min(data$cohort) # third used only in naming later per1 <- min(data$period) # use minimum values to construct i, j, and k i.value <- (data$age - age1)/unit + 1 + n.age.excl.start k.value <- (data$cohort - coh1)/unit + 1 + n.coh.excl.start j.value <- i.value + k.value - 1 + n.per.excl.start # get unique combinations of i, j, and k; i.e. list of ik cells data.wijk <- cbind(data, i.value, j.value, k.value) mu.index <- plyr::count(data.wijk, c("i.value", "j.value", "k.value")) mu.index <- mu.index[mu.index$freq>0, ] mu.index <- mu.index[, !names(mu.index)=="freq"] # define L and U: anchor points of plane L <- min(j.value) - 1 # correct bc this is the minimum observed j U <- as.integer((L+3)/2) # stop if something has gone wrong here stopifnot(exists("L")==TRUE) stopifnot(exists("U")==TRUE) # get dimensionality of data I <- as.numeric(length(unique(data.wijk$i.value))) J <- as.numeric(length(unique(data.wijk$j.value))) K <- as.numeric(length(unique(data.wijk$k.value))) ## Step 3: create and fill design matrix - one row per ik cell #### #create space of design matrix n.rows <- nrow(mu.index) n.param.design <- I+J+K-2 # -2 rather than -3 because we allow for a period # slope; in the submodels with P and tP, it's easier to estimate a period # slope than an equality constraint design <- matrix(data=0, nrow=n.rows, ncol=n.param.design+3) # the extra three are for index: i, j, k #begin by defining L.odd L.odd <- !L %% 2 == 0 #fill row-by-row for (row in 1:n.rows) { i <- mu.index[row, 1] j <- mu.index[row, 2] k <- mu.index[row, 3] design[row, 1] <- 1 # the intercept design[row, 2] <- i - U # the first (age) slope design[row, 4] <- k - U # the second (cohort) slope design[row, 3] <- design[row, 2] + design[row, 4] # period slope #age DDs # note 4 plane pieces => start at 4+1 if (i < U) design[row, (4 + i - n.age.excl.start): (4 + U - 1 - n.age.excl.start)] <- seq(1, U - i) # backward cumulation from DD3 to DD(U-1) if (i > U + 1) design[row, (4 + U - n.age.excl.start): (4 + i - 2 - n.age.excl.start)] <- seq(i - U - 1, 1) # forward cumulation from DDU to DDI # note DDI is only the (I-2)nd column of DDs b/c start at DD3 #period DDs # note 4 plane pieces + I-2 age DDs => start at 2+I+1 # CASE 1: L odd if (L.odd && j == (L + 1)) design[row, (2 + I + 1)] <- 1 # single backward cumulation of DD(L+3) if (L.odd && j > (L + 3)) design[row, (2 + I + 1 + L.odd):(2 + I + L.odd + j - (L + 3))] <- seq(j - (L + 3), 1) # forward cumulation from DD(L+4) to DD(L+J) # note DD(L+J) is the (J-3)rd column of DDs here b/c start at DD(L+4); # recall DD(L+3) already covered with backward cumulation # CASE 2: L even if ((!L.odd) && j > (L + 2)) design[row, (2 + I + 1):(2 + I + j - (L + 2))] <- seq(j - (L + 2), 1) # forward cumulation from DD(L+3) to DD(L+J) # DD(L+J) is the (J-2)nd column of DDs here b/c start at DD(L+3) #cohort DDs # note 4 plane pieces + I-2 age + J-2 period DDs => start at I+J+1 # also note symmetry to age DDs if (k < U) design[row, (I + J + k - n.coh.excl.start): (I + J + U - 1 - n.coh.excl.start)] <- seq(1, U - k) # backward cumulation from DD3 to DD(U-1) if (k > U + 1) design[row, (I + J + U - n.coh.excl.start): (I + J + k - 2 - n.coh.excl.start)] <- seq(k - U - 1, 1) # forward cumulation from DDU to DDI # note DDI is only the (I-2)nd column of DDs b/c start at DD3 #index design[row, n.param.design+1] <- i design[row, n.param.design+2] <- j design[row, n.param.design+3] <- k } # convert matrix -> data frame designdf <- as.data.frame(design) # assign names to variables in design dataframe firstfour <- c("level", "age slope", "period slope", "cohort slope") age.names <- paste("DD_age", seq(age1+2*unit, age1+2*unit+(I-2)*unit-1, unit), sep="_") per.names <- paste("DD_period", seq(per1+2*unit, per1+2*unit+(J-2)*unit-1, unit), sep="_") coh.names <- paste("DD_cohort", seq(coh1+2*unit, coh1+2*unit+(K-2)*unit-1, unit), sep="_") index <- c("i.value", "j.value", "k.value") colnames(designdf) <- c(firstfour, age.names, per.names, coh.names, index) ## Step 4: join design matrix with full data using ijk index values #### # associates APC design matrix with each individual datapoint as needed data.with.design <- join(data.wijk, designdf, by=c("i.value", "j.value", "k.value")) # drop information that is redundant in future commands exclude.all <- c("X", "indiv.ID", "age", "period", "cohort", "i.value", "j.value", "k.value") first.exclusion <- names(data.with.design) %in% c(exclude.all) full.design.collinear <- data.with.design[!first.exclusion] # check join was successful if (anyNA(full.design.collinear)) stop("NA in design matrix, check data or n.excl age/per/coh") # --------------------------------------------------------------------------- # Return valuables valuables <- list(structure.design.collinear = designdf, full.design.collinear = full.design.collinear, unit = unit, age1 = age1, per1 = per1, coh1 = coh1, age.max = I, per.max = J, coh.max = K, per.zero = L, per.odd = L.odd, U = U) return (valuables) } #### write apc.indiv.design.model ##### apc.indiv.design.model <- function(apc.indiv.design.collinear, model.design = "APC", dep.var = NULL, covariates = NULL){ ################## # This function reduces the collinear design matrix to the design # matrix for the desired model ################## # Step 1: set-up # check model design design.list <- c("APC", "AP", "AC", "PC", "Ad", "Pd", "Cd", "A", "P", "C", "t", "tA", "tP", "tC", "1") if (!isTRUE(model.design %in% design.list)) stop("model.design not recognised") # Extract elements from apc.indiv.design.collinear structure.design.collinear <- apc.indiv.design.collinear$structure.design.collinear full.design.collinear <- apc.indiv.design.collinear$full.design.collinear # check covariates missing.covariates <- covariates[!covariates %in% names(full.design.collinear)] if(!length(missing.covariates)==0) stop("some listed covariates are not in dataset") # get double differences set.coh.DDs <- colnames(structure.design.collinear)[grep( "^DD_cohort_", colnames(structure.design.collinear))] set.age.DDs <- colnames(structure.design.collinear)[grep( "^DD_age_", colnames(structure.design.collinear))] set.per.DDs <- colnames(structure.design.collinear)[grep( "^DD_period_", colnames(structure.design.collinear))] ################## # Step 2: Select model inclusions if(model.design=="APC") { slopes <- c(1,0,1); difdif <- c(1,1,1); } if(model.design=="AP" ) { slopes <- c(1,0,1); difdif <- c(1,1,0); } if(model.design=="AC" ) { slopes <- c(1,0,1); difdif <- c(1,0,1); } if(model.design=="PC" ) { slopes <- c(1,0,1); difdif <- c(0,1,1); } if(model.design=="Ad" ) { slopes <- c(1,0,1); difdif <- c(1,0,0); } if(model.design=="Pd" ) { slopes <- c(1,0,1); difdif <- c(0,1,0); } if(model.design=="Cd" ) { slopes <- c(1,0,1); difdif <- c(0,0,1); } if(model.design=="A" ) { slopes <- c(1,0,0); difdif <- c(1,0,0); } if(model.design=="P" ) { slopes <- c(0,1,0); difdif <- c(0,1,0); } if(model.design=="C" ) { slopes <- c(0,0,1); difdif <- c(0,0,1); } if(model.design=="t" ) { slopes <- c(1,0,1); difdif <- c(0,0,0); } if(model.design=="tA" ) { slopes <- c(1,0,0); difdif <- c(0,0,0); } if(model.design=="tP" ) { slopes <- c(0,1,0); difdif <- c(0,0,0); } if(model.design=="tC" ) { slopes <- c(0,0,1); difdif <- c(0,0,0); } if(model.design=="1" ) { slopes <- c(0,0,0); difdif <- c(0,0,0); } incl <- vector(mode="character") incl <- c(incl, "level") if (slopes[1]) {incl <- c(incl, "age slope")} if (slopes[2]) {incl <- c(incl, "period slope")} if (slopes[3]) {incl <- c(incl, "cohort slope")} if (difdif[1]) {incl <- c(incl, set.age.DDs)} if (difdif[2]) {incl <- c(incl, set.per.DDs)} if (difdif[3]) {incl <- c(incl, set.coh.DDs)} full.inclusion <- names(full.design.collinear) %in% c(incl, covariates) final.design <- full.design.collinear[full.inclusion] # warn if no DV specification. No stop because can introduce DV in fit stage if (is.null(dep.var)){ warning("Dependent variable unspecified") DV <- dep.var # ie assign NULL } else DV <- full.design.collinear[dep.var] ################## # Return valuables valuables <- list(full.design = final.design, DV = DV, slopes = slopes, difdif = difdif, xi.dim = length(incl), model.design = model.design, unit = apc.indiv.design.collinear$unit, age1 = apc.indiv.design.collinear$age1, per1 = apc.indiv.design.collinear$per1, coh1 = apc.indiv.design.collinear$coh1, age.max = apc.indiv.design.collinear$age.max, per.max = apc.indiv.design.collinear$per.max, coh.max = apc.indiv.design.collinear$coh.max, per.zero = apc.indiv.design.collinear$per.zero, per.odd = apc.indiv.design.collinear$per.odd, U = apc.indiv.design.collinear$U) return(valuables) } #### write apc.indiv.fit.model ##### apc.indiv.fit.model <- function (apc.indiv.design.model, model.family=NULL, DV=NULL){ ################## # This function estimates the model developed in apc.indiv.design.model ################## # Step 1: Set up # Check model.design family.list <- c("gaussian", "binomial") if (!isTRUE(model.family %in% family.list)) stop("model.family must be either 'gaussian' or 'binomial'") # Extract elements from apc.indiv.design.model design <- apc.indiv.design.model$full.design if(is.null(DV)) { DV <- apc.indiv.design.model$DV} xi.dim1 <- apc.indiv.design.model$xi.dim model.design <- apc.indiv.design.model$model.design difdif <- apc.indiv.design.model$difdif slopes <- apc.indiv.design.model$slopes age.max <- apc.indiv.design.model$age.max per.max <- apc.indiv.design.model$per.max coh.max <- apc.indiv.design.model$coh.max age1 <- apc.indiv.design.model$age1 per1 <- apc.indiv.design.model$per1 coh1 <- apc.indiv.design.model$coh1 unit <- apc.indiv.design.model$unit if (is.null(DV)) stop("Dependent variable unspecified") ################## # Step 2: Estimation # Regression using glm.fit if (model.family == "binomial") fit <- glm.fit(design, DV[, 1], family=binomial()) if (model.family == "gaussian") fit <- glm.fit(design, DV[, 1], family=gaussian()) # Extract coefficients and covariance coefficients <- summary.glm(fit)$coefficients covariance <- summary.glm(fit)$cov.scaled n.coeff <- length(coefficients[, 1]) n.coeff.canonical <- xi.dim1 # Manipulate coefficients and covariance based on model design # First case: model design is 1 if (model.design == "1"){ coefficients.canonical <- t(coefficients[(n.coeff- n.coeff.canonical+1):n.coeff, ]) covariance.canonical <- summary.glm(fit)$cov.scaled[(n.coeff- n.coeff.canonical+1):n.coeff, (n.coeff- n.coeff.canonical+1):n.coeff] if(n.coeff - n.coeff.canonical > 0) coefficients.covariates <- t(coefficients[1:(n.coeff-n.coeff.canonical), ]) else coefficients.covariates <- NULL # get standard errors coefficients.canonical[,2] <- sqrt(covariance.canonical) # get t-statistics coefficients.canonical[,3] <- coefficients.canonical[,1] / coefficients.canonical[,2] # get p-values coefficients.canonical[,4] <- 2*pnorm(abs(coefficients.canonical[ ,3]), lower.tail=FALSE) } # Second case: model design is not 1 if (!model.design == "1"){ coefficients.canonical <- summary.glm(fit)$coefficients[(n.coeff- n.coeff.canonical+1):n.coeff, ] covariance.canonical <- summary.glm(fit)$cov.scaled[(n.coeff- n.coeff.canonical+1):n.coeff, (n.coeff-n.coeff.canonical+1):n.coeff] if(n.coeff - n.coeff.canonical > 0) coefficients.covariates <- t(coefficients[1:(n.coeff-n.coeff.canonical), ]) else coefficients.covariates <- NULL # get standard errors coefficients.canonical[,2] <- sqrt(diag(covariance.canonical)) # get t-statistics coefficients.canonical[,3] <- coefficients.canonical[,1] / coefficients.canonical[,2] # get p-values coefficients.canonical[,4] <- 2*pnorm(abs(coefficients.canonical[ ,3]), lower.tail=FALSE) } ################## # Step 3: Use output thus far to get dates for use in plotting later index.age <- NULL index.per <- NULL index.coh <- NULL # changes to index based on max start <- 1+sum(slopes) if(difdif[1]) { index.age <- start+seq(1,age.max-2);start <- start+age.max-2 } if(difdif[2]) { index.per <- start+seq(1,per.max-2);start <- start+per.max-2 } if(difdif[3]) { index.coh <- start+seq(1,coh.max-2);start <- start+coh.max-2 } xi.dim2 <- start # use index, max, unit, xi.dim2, difdif to get dates dates <- matrix(data=NA,nrow=xi.dim2,ncol=1) if(difdif[1]) dates[index.age,1] <- age1+seq(2,age.max-1)*unit if(difdif[2]) dates[index.per,1] <- per1+seq(2,per.max-1)*unit if(difdif[3]) dates[index.coh,1] <- coh1+seq(2,coh.max-1)*unit ################## # Construct likelihood RSS <- sum(fit$residuals^2) N <- length(fit$residuals) lik <- -(N/2)*log(2*pi*RSS/N) - (N/2) # Return valuables valuables <- c(fit, list(coefficients.canonical = coefficients.canonical, covariance.canonical = covariance.canonical, dates = dates, index.age = index.age, index.coh = index.coh, index.per = index.per, coefficients.covariates = coefficients.covariates, model.family = model.family, difdif = difdif, slopes = slopes, age1 = age1, per1 = per1, coh1 = coh1, unit = unit, age.max = age.max, per.max = per.max, coh.max = coh.max, model.design = apc.indiv.design.model$model.design, per.zero = apc.indiv.design.model$per.zero, per.odd = apc.indiv.design.model$per.odd, U = apc.indiv.design.model$U, likelihood = lik)) return(valuables) } #### write apc.indiv.ftest.fullapc ##### apc.indiv.ftest.fullapc <- function(data, big.model="APC", small.model, dep.var, covariates=NULL, model.family="gaussian", unit=1, n.coh.excl.start=0, n.coh.excl.end=0, n.age.excl.start=0, n.age.excl.end=0, n.per.excl.start=0, n.per.excl.end=0, existing.big.model.fit=NULL, existing.small.model.fit=NULL, existing.collinear=NULL){ ###################### # This functions runs an F-test comparing the fit of a single # APC submodel to the full APC model. ###################### # Step 1: Fit reduced APC model (the submodel) if(is.null(existing.small.model.fit)){ if(is.null(existing.collinear)){ collinear.R <- apc.indiv.design.collinear(data, unit = unit, n.coh.excl.start = n.coh.excl.start, n.coh.excl.end = n.coh.excl.end, n.age.excl.start = n.age.excl.start, n.age.excl.end = n.age.excl.end, n.per.excl.start = n.per.excl.start, n.per.excl.end = n.per.excl.end) } else collinear.R <- existing.collinear design.R <- apc.indiv.design.model(apc.indiv.design.collinear =collinear.R, model.design = small.model, dep.var=dep.var, covariates=covariates) fitted.R <- apc.indiv.fit.model(design.R, model.family=model.family) } else fitted.R <- existing.small.model.fit # Get RSS, number of parameters RSS.APC.R <- sum(fitted.R$residuals^2) nparam.APC.R <- length(fitted.R$coefficients) ###################### # Step 2: Fit full APC model if(is.null(existing.big.model.fit)){ # note: can use same collinear design matrix for full model as for submodel design.U <- apc.indiv.design.model(apc.indiv.design.collinear =collinear.R, model.design = big.model, dep.var=dep.var, covariates=covariates) fitted.U <- apc.indiv.fit.model(design.U, model.family=model.family) } else fitted.U <- existing.big.model.fit # Get RSS, number of parameters, number of observations RSS.APC.U <- sum(fitted.U$residuals^2) nparam.APC.U <- length(fitted.U$coefficients) sample.N.APC.U <- length(fitted.U$y) ###################### # Step 3: Run F-test # Get elements of F-test big.model.RSS <- RSS.APC.U small.model.RSS <- RSS.APC.R big.model.nparams <- nparam.APC.U small.model.nparams <- nparam.APC.R sample.N <- sample.N.APC.U df.num <- big.model.nparams - small.model.nparams df.denom <- sample.N - big.model.nparams df <- paste("(", paste(as.character(df.num), as.character(df.denom), sep=", "), ")", sep="") # Construct F-statistic fstat.num <- (small.model.RSS - big.model.RSS)/df.num fstat.denom <- (big.model.RSS - 0)/df.denom fstat <- fstat.num/fstat.denom # Get p-value by comparing F-statistic to F-distribution p.value <- pf(fstat, df1 = df.num, df2 = df.denom, lower.tail = FALSE) ###################### # valuables valuables <- list(fstat = fstat, df = df, df.num = df.num, df.denom = df.denom, p.value = p.value, aic.big = fitted.U$aic, aic.small = fitted.R$aic, lik.big = fitted.U$likelihood, lik.small = fitted.R$likelihood) return(valuables) } #### write var.apc.plot.fit #### var.apc.plot.fit <- function(apc.fit.model,scale=FALSE,sdv.at.zero=TRUE, type="detrend",sub.plot=NULL,main.outer=NULL, main.sub=NULL,cex=NULL,cex.axis=NULL, cex.main=2, mgp=c(2, 1, 0), theight=1, mar=c(4, 3, 2, 0), oma = c(0, 0, 5, 1)) { ################## # This fuction is the same as apc.plot.fit, with a few alterations to visuals ################## # change type if(type=="ss.dd") type<-"sum.sum" ################# # get values from fit coefficients.canonical <- apc.fit.model$coefficients.canonical slopes <- apc.fit.model$slopes difdif <- apc.fit.model$difdif index.age <- apc.fit.model$index.age index.per <- apc.fit.model$index.per index.coh <- apc.fit.model$index.coh dates <- apc.fit.model$dates model.design <- apc.fit.model$model.design model.family <- apc.fit.model$model.family age1 <- apc.fit.model$age1 per1 <- apc.fit.model$per1 coh1 <- apc.fit.model$coh1 unit <- apc.fit.model$unit age.max <- apc.fit.model$age.max per.max <- apc.fit.model$per.max coh.max <- apc.fit.model$coh.max per.odd <- apc.fit.model$per.odd U <- apc.fit.model$U ################# # identify fit apc.id <- apc.identify(apc.fit.model) index.age.max <- apc.id$index.age.max index.per.max <- apc.id$index.per.max index.coh.max <- apc.id$index.coh.max dates.max <- apc.id$dates.max index.age.sub <- apc.id$index.age.sub index.per.sub <- apc.id$index.per.sub index.coh.sub <- apc.id$index.coh.sub dates.sub <- apc.id$dates.sub index.age.dif <- apc.id$index.age.dif index.per.dif <- apc.id$index.per.dif index.coh.dif <- apc.id$index.coh.dif dates.dif <- apc.id$dates.dif coefficients.ssdd <- apc.id$coefficients.ssdd coefficients.detrend <- apc.id$coefficients.detrend coefficients.demean <- apc.id$coefficients.demean coefficients.dif <- apc.id$coefficients.dif ############################## # check model design if(isTRUE(type %in% c("dif"))) { if(model.design=="APC") return(cat("apc.error: differences not identified when model.design is APC. Type cannot be demean or dif \n")) if(model.design %in% c("Ad","Pd","Cd","A","P","C","t","tA","tP","tC","1")) return(cat("apc.error: types demean and dif not implemented for model designs At, Pt, Ct, A, P, C, t, tA, tP, tC, 1 \n")) } ########################################### # construct ingredients to plot depending on type mixed <- FALSE if(model.family=="poisson.response") mixed <- TRUE ########################################### # declare variables v.do.plot <- vector(length=9) l.dates <- list(a=1,b=1,c=1,d=1,e=1,f=1,g=1,h=1,i=1) l.coefficients <- list(a=1,b=1,c=1,d=1,e=1,f=1,g=1,h=1,i=1) v.main.sub <- vector(length=9) v.xlab <- vector(length=9) v.intercept <- vector(length=9) v.tau <- vector(length=9) ########################################### # type is "detrend" or "sum.sum" if(type %in% c("detrend","sum.sum")) { if(type=="detrend") main <- paste("APC canonical parameters & detrended representation","\n", "model.design=",model.design, "(1/2 std blue/red)") if(type=="sum.sum") main <- paste("APC canonical parameters & standard representation","\n", "model.design=",model.design, "(1/2 std blue/red)") ################## # do plot? v.do.plot[1:3] <- difdif v.do.plot[4] <- isTRUE(model.design %in% c("APC","AP","AC","PC","Ad","Pd", "Cd","A","P","t","tA","tP")) v.do.plot[5] <- TRUE v.do.plot[6] <- isTRUE(model.design %in% c("APC","AP","AC","PC","Ad","Pd", "Cd","C","t","tC")) v.do.plot[7:9] <- difdif ################## # sub.main v.main.sub[1] <- expression(paste("(a) ",Delta^2,alpha)) v.main.sub[2] <- expression(paste("(b) ",Delta^2,beta)) v.main.sub[3] <- expression(paste("(c) ",Delta^2,gamma)) if(model.design %in% c("APC","AP","AC","PC","Ad","Pd","Cd","t","A","tA")) v.main.sub[4] <- "(d) first linear trend" if(model.design %in% c("A","tA")) v.main.sub[4] <-"(d) age linear trend" if(model.design %in% c("P","tP")) v.main.sub[4] <- "(d) period linear trend" if(!mixed) v.main.sub[5] <- "(e) level" if(mixed) v.main.sub[5] <- "(e) aggregate mean" if(model.design %in% c("APC","AP","AC","PC","Ad","Pd","Cd","t")) v.main.sub[6] <- "(f) second linear trend" if(model.design %in% c("C","tC")) v.main.sub[6] <- "(f) cohort linear trend" if(type=="detrend") { v.main.sub[7] <- expression(paste("(g) detrended ",Sigma^2,Delta^2,alpha)) v.main.sub[8] <- expression(paste("(h) detrended ",Sigma^2,Delta^2,beta)) v.main.sub[9] <- expression(paste("(i) detrended ",Sigma^2,Delta^2,gamma)) } if(type=="sum.sum") { v.main.sub[7] <- expression(paste("(g) ",Sigma^2,Delta^2,alpha)) v.main.sub[8] <- expression(paste("(h) ",Sigma^2,Delta^2,beta)) v.main.sub[9] <- expression(paste("(i) ",Sigma^2,Delta^2,gamma)) } ################## # dates l.dates[[1]] <- dates[index.age,1] l.dates[[2]] <- dates[index.per,1] l.dates[[3]] <- dates[index.coh,1] if(model.design %in% c("APC","AP","AC","PC","Ad","Pd","Cd","t","A","tA")) l.dates[[4]] <- age1+seq(0,age.max-1)*unit if(model.design %in% c("P","tP")) l.dates[[4]] <- per1+seq(0,per.max-1)*unit l.dates[[5]] <- c(0,1) # matrix(data=c(0,1) ,nrow=2 ,ncol=1) l.dates[[6]] <- coh1+seq(0,coh.max-1)*unit l.dates[[7]] <- dates.max[index.age.max,1] l.dates[[8]] <- dates.max[index.per.max,1] l.dates[[9]] <- dates.max[index.coh.max,1] ################## # coefficients l.coefficients[[1]] <- coefficients.canonical[index.age,] l.coefficients[[2]] <- coefficients.canonical[index.per,] l.coefficients[[3]] <- coefficients.canonical[index.coh,] if(type=="detrend") { coefficients.sum.sum <- coefficients.detrend UU <- 1 } if(type=="sum.sum") { coefficients.sum.sum <- coefficients.ssdd UU <- U } if(model.design %in% c("APC","AP","AC","PC","Ad","Pd","Cd","t","A","tA")) l.coefficients[[4]] <- matrix(data=seq(1,age.max)-UU , nrow=age.max,ncol=1) %*% coefficients.sum.sum[2,] if(model.design %in% c("P","tP")) l.coefficients[[4]] <- matrix(data=seq(1,per.max)-per.odd-1 , nrow=per.max,ncol=1) %*% coefficients.sum.sum[2,] l.coefficients[[5]] <- matrix(data=c(1,1) , nrow=2 ,ncol=1) %*% coefficients.sum.sum[1,] if(model.design %in% c("APC","AP","AC","PC","Ad","Pd","Cd","t")) l.coefficients[[6]] <- matrix(data=seq(1,coh.max)-UU , nrow=coh.max,ncol=1) %*% coefficients.sum.sum[3,] if(model.design %in% c("C","tC")) l.coefficients[[6]] <- matrix(data=seq(1,coh.max)-UU , nrow=coh.max,ncol=1) %*% coefficients.sum.sum[2,] l.coefficients[[7]] <- coefficients.sum.sum[index.age.max,] l.coefficients[[8]] <- coefficients.sum.sum[index.per.max,] l.coefficients[[9]] <- coefficients.sum.sum[index.coh.max,] #################### # xlab v.xlab[1:3] <- c("age","period","cohort") if(model.design %in% c("APC","AP","AC","PC","Ad","Pd","Cd","t","A","tA")) v.xlab[4] <- "age" if(model.design %in% c("P","tP")) v.xlab[4] <- "period" if(!mixed) v.xlab[5] <- "age, period, cohort" if(mixed) v.xlab[5] <- "" v.xlab[6] <- "cohort" v.xlab[7:9] <- c("age","period","cohort") #################### # intercept v.intercept[5] <- TRUE #################### # tau if(mixed) v.tau[5] <- TRUE } ########################################### # type is "dif" if(type %in% c("dif")) { main <- paste("Difference parameters & demeaned representation", "\n","model.design=",model.design, "(1/2 std blue/red)") ################## # do plot? v.do.plot[1:3] <- difdif v.do.plot[5] <- TRUE v.do.plot[7:9] <- difdif ################## # sub.main v.main.sub[1] <- expression(paste("(a) ",Delta,alpha)) v.main.sub[2] <- expression(paste("(b) ",Delta,beta)) v.main.sub[3] <- expression(paste("(c) ",Delta,gamma)) if(!mixed) v.main.sub[5] <- "(e) level" if(mixed) v.main.sub[5] <- "(e) aggregate mean" v.main.sub[7] <- expression(paste("(g) demeaned ",Sigma,Delta,alpha)) v.main.sub[8] <- expression(paste("(h) demeaned ",Sigma,Delta,beta)) v.main.sub[9] <- expression(paste("(i) demeaned ",Sigma,Delta,gamma)) ################## # dates l.dates[[1]] <- dates.dif[index.age.dif,1] l.dates[[2]] <- dates.dif[index.per.dif,1] l.dates[[3]] <- dates.dif[index.coh.dif,1] l.dates[[5]] <- c(0,1) # matrix(data=c(0,1) ,nrow=2 ,ncol=1) l.dates[[7]] <- dates.sub[index.age.sub,1] l.dates[[8]] <- dates.sub[index.per.sub,1] l.dates[[9]] <- dates.sub[index.coh.sub,1] ################## # coefficients l.coefficients[[1]] <- coefficients.dif[index.age.dif,] l.coefficients[[2]] <- coefficients.dif[index.per.dif,] l.coefficients[[3]] <- coefficients.dif[index.coh.dif,] l.coefficients[[5]] <- matrix(data=c(1,1) ,nrow=2 , ncol=1) %*% coefficients.demean[1,] l.coefficients[[7]] <- coefficients.demean[index.age.sub,] l.coefficients[[8]] <- coefficients.demean[index.per.sub,] l.coefficients[[9]] <- coefficients.demean[index.coh.sub,] #################### # xlab v.xlab[1:3] <- c("age","period","cohort") if(!mixed) v.xlab[5] <- "age, period, cohort" if(mixed) v.xlab[5] <- "" v.xlab[7:9] <- c("age","period","cohort") #################### # intercept v.intercept[5] <- TRUE #################### # tau if(mixed) v.tau[5] <- TRUE } ############################################# # use arguments of function if(is.null(main.outer)==FALSE) main <- main.outer if(is.null(main.sub)==FALSE) v.main.sub <- main.sub if(scale==1 && model.family=="binomial.dose.response") scale <- 2 ####################################################### # Internal function to plot estimates with sdv ####################################################### function.plot.est.sdv <- function(dates,coefficients,xlab="",main="", scale=0,sdv.at.zero=FALSE,intercept=FALSE, tau=FALSE,cex=NULL,cex.axis=NULL) # BN 2 Dec 2013 { # function.plot.est.sdv # define function that can move to exponential scale function.scale <- function(x,scale=0) { if(scale==0) x.scale <- x if(scale==1) x.scale <- exp(x) if(scale==2) x.scale <- exp(x)/(1+exp(x)) return(x.scale) } ################ # IF MORE THAN ONE OBSERVATION if(length(dates)>1) { ################ # get estimates and sdv dat <- as.vector(dates) est <- as.vector(coefficients[,1]) sdv <- as.vector(coefficients[,2]) # get center for sdv sdv0 <- (1-sdv.at.zero)*est ################ # set ylim y.lower <- min(0,min(function.scale(est,scale)),max(2*min(function.scale( est,scale)),max(function.scale(sdv0-sdv,scale),na.rm=TRUE))) y.upper <- max(0,max(function.scale(est,scale)),min(2*max(function.scale( est,scale)),min(function.scale(sdv0+sdv,scale),na.rm=TRUE))) if(max(est)-min(est)|z|)") rownames(coefficients.TS) <- rownames(Sigma.hat)[1:nrow(kappahat)] coefficients.covariates <- matrix(nrow=nrow(as.matrix(zetahat)), ncol=4) coefficients.covariates[,1] <- zetahat coefficients.covariates[,2] <- Std.error[(nrow(kappahat)+1):nrow(Sigma.hat)] coefficients.covariates[,3] <- coefficients.covariates[,1]/coefficients.covariates[,2] coefficients.covariates[,4] <- 2*pnorm(abs(coefficients.covariates[,3]), lower.tail=FALSE) colnames(coefficients.covariates) <- c("Estimate", "Std. Error", "z value", "Pr(>|z|)") rownames(coefficients.covariates) <- rownames(Sigma.hat)[(nrow(kappahat)+1):nrow(Sigma.hat)] } ###################### ##### Case without covariates else{ # Effective regression of Y on T rhohat.withname <- ddply(data, .variables=c("cell.name"), function(subdata, colnames){colMeans(subdata[colnames])}, dep.var) colnames(rhohat.withname)[2] <- "rhohat" rhohat <- rhohat.withname[, 2] data.all.coeffs <- join(data, rhohat.withname, by="cell.name") # Construct residuals from regression of Y on T y.residuals <- data.all.coeffs[dep.var] - data.all.coeffs["rhohat"] colnames(y.residuals) <- "y.residuals" final.data <- cbind(data.all.coeffs, y.residuals) zetahat <- NA kappahat <- rhohat # Step 6: Calculate covariance matrix # Homoskedasticity assumed so can use var(beta) = (X'X)^{-1} sigma^2 # get the estimate of the variance of the residuals small.sigma.hat <- sum(final.data["y.residuals"]^2)/nrow(final.data) # in case there are more levels than unique values observed (would be e.g. # if using data where cell.name already defined and it was a subset) final.data$cell.name <- droplevels(data$cell.name) # get parts of design matrix Tdiag.withname <- ddply(data, .variables=c("cell.name"), function(subdata){(nrow(subdata))}) colnames(Tdiag.withname)[2] <- "cell.counts" Tdiag <- Tdiag.withname[,2] # Invert top left (T'T) Tinv.diag <- 1/Tdiag invD <- diag(Tinv.diag) Sigma.hat <- invD*small.sigma.hat Std.error <- sqrt(diag(Sigma.hat)) coefficients.TS <- matrix(nrow=nrow(as.matrix(kappahat)), ncol=4) coefficients.TS[,1] <- kappahat coefficients.TS[,2] <- Std.error[1:nrow(as.matrix(kappahat))] coefficients.TS[,3] <- coefficients.TS[,1]/coefficients.TS[,2] coefficients.TS[,4] <- 2*pnorm(abs(coefficients.TS[,3]), lower.tail=FALSE) colnames(coefficients.TS) <- c("Estimate", "Std. Error", "z value", "Pr(>|z|)") rownames(coefficients.TS) <- rownames(Sigma.hat)[1:nrow(as.matrix(kappahat))] coefficients.covariates <- NULL } ###################### # Calculate likelihood and AIC RSS <- sum(final.data["y.residuals"]^2) n.coeff <- length(rhohat) + length(zetahat) - as.numeric(is.na(zetahat[1])) #because length(NA)=1 N <- nrow(final.data) lik <- -(N/2)*log(2*pi*RSS/N) - (N/2) aic <- -2*lik + 2*(n.coeff + 1) # Tidy names zetahat old.names.zetahat <- names(zetahat) new.names.zetahat <- gsub("^resid.", "", old.names.zetahat) names(zetahat) <- new.names.zetahat ###################### # Return valuables valuables <- list(n.cells = length(unique(final.data$cell.name)), TSS = sum(final.data[dep.var]^2), RSS = RSS, zetahat = zetahat, rhohat = rhohat, kappahat = kappahat, coefficients.covariates = coefficients.covariates, coefficients.TS = coefficients.TS, Sigmahat = Sigma.hat, likelihood = lik, aic = aic) } #### write apc.indiv.ftest.TS ##### apc.indiv.ftest.TS <- function(data, small.model="APC", dep.var, covariates=NULL, model.family="gaussian", unit=1, n.coh.excl.start=0, n.coh.excl.end = 0, n.age.excl.start=0, n.age.excl.end = 0, n.per.excl.start=0, n.per.excl.end = 0, existing.small.model.fit=NULL, existing.big.model.fit=NULL, existing.collinear=NULL){ ###################### # This function runs an F-test comparing the fit of a single # APC model or submodel to the temporally-saturated model. ###################### # Step 1: Fit APC model if(is.null(existing.small.model.fit)){ if(is.null(existing.collinear)){ collinear <- apc.indiv.design.collinear(data, unit = unit, n.coh.excl.start = n.coh.excl.start, n.coh.excl.end = n.coh.excl.end, n.age.excl.start = n.age.excl.start, n.age.excl.end = n.age.excl.end, n.per.excl.start = n.per.excl.start, n.per.excl.end = n.per.excl.end) } else collinear <- existing.collinear design <- apc.indiv.design.model(collinear, model.design = small.model, dep.var=dep.var, covariates=covariates) fitted <- apc.indiv.fit.model(design, model.family=model.family) } else fitted <- existing.small.model.fit # get RSS, number of parameters, number of observations for APC model. RSS.APC <- sum(fitted$residuals^2) nparam.APC <- length(fitted$coefficients) sample.N.APC <- length(fitted$y) ###################### # Step 2: Fit TS model if(is.null(existing.big.model.fit)){ dummymodel <- apc.indiv.estimate.TS(data, dep.var = dep.var, covariates = covariates) } else dummymodel <- existing.big.model.fit # get RSS, number of parameters from TS model. RSS.dummy <- dummymodel$RSS if (!is.na(dummymodel$zetahat[1])){ nparam.dummy <- dummymodel$n.cells+length(dummymodel$zetahat) } else{ nparam.dummy <- dummymodel$n.cells } ###################### # step 3: Run F-test # Elements of F-test big.model.RSS <- RSS.dummy small.model.RSS <- RSS.APC big.model.nparams <- nparam.dummy small.model.nparams <- nparam.APC sample.N <- sample.N.APC df.num <- big.model.nparams - small.model.nparams df.denom <- sample.N - big.model.nparams df <- paste("(", paste(as.character(df.num), as.character(df.denom), sep=", "), ")", sep="") # Combine elements in F-statistic fstat.num <- (small.model.RSS - big.model.RSS)/df.num fstat.denom <- (big.model.RSS - 0)/df.denom fstat <- fstat.num/fstat.denom # Get p-value by comparing F-statistic to F-distribution p.value <- pf(fstat, df1 = df.num, df2 = df.denom, lower.tail = FALSE) ###################### # valuables valuables <- list(fstat = fstat, df = df, df.num = df.num, df.denom = df.denom, p.value = p.value, aic.small = fitted$aic, aic.big = dummymodel$aic, lik.small = fitted$likelihood, lik.big = dummymodel$likelihood) return(valuables) } #### write apc.indiv.ftable ##### apc.indiv.ftable <- function(data, dep.var, covariates=NULL, unit=1, n.coh.excl.start=0, n.coh.excl.end=0, n.age.excl.start=0, n.age.excl.end=0, n.per.excl.start=0, n.per.excl.end=0){ ###################### # This function generates a table comparing all APC submodels to the # TS model and the full APC model in terms of an F-test of fit. ###################### # Step 1: write function that uses output from an F-test against the #TS and an F-test against the full # APC model to generate a line suitable for the final table. get.table.line <- function(ftest1, ftest2, first=FALSE, second=FALSE){ if(isTRUE(first)) # first line is TS model so aic only line <- c(NA, NA, NA, NA, NA, NA, round(ftest1$aic.big, 3), round(ftest1$lik.big, 3)) else if(isTRUE(second)) # second line is full APC model so not compared to itself line <- c(round(ftest1$fstat, 3), ftest1$df.num, round(ftest1$p.value,3), NA, NA, NA, round(ftest1$aic.small,3), round(ftest1$lik.small, 3)) else line <- c(round(ftest1$fstat, 3), ftest1$df.num, round(ftest1$p.value,3), round(ftest2$fstat, 3), ftest2$df.num, round(ftest2$p.value, 3), round(ftest2$aic.small,3), round(ftest2$lik.small, 3)) return(line) } ###################### # Step 2: Other setup # List of submodels to be investigated model.design.list <- c("TS", "APC","AP","AC","PC","Ad","Pd","Cd","A","P", "C","t","tA","tP","tC", "1") # Empty table fit.tab <- matrix(nrow=length(model.design.list),ncol=8,data=NA) # Full APC model for comparison fullAPCcollinear <- apc.indiv.design.collinear(data, unit=unit, n.coh.excl.start=n.coh.excl.start, n.coh.excl.end=n.coh.excl.end, n.age.excl.start = n.age.excl.start, n.age.excl.end = n.age.excl.end, n.per.excl.start = n.per.excl.start, n.per.excl.end = n.per.excl.end) fullAPCdesign <- apc.indiv.design.model(fullAPCcollinear, dep.var = dep.var, covariates = covariates) fullAPCfit <- apc.indiv.fit.model(fullAPCdesign, model.family="gaussian") # TS model for comparison TSfit <- apc.indiv.estimate.TS(data, dep.var=dep.var, covariates=covariates) ###################### # Step 3: Generate table # First line: TS model ftest.full.vs.TS <- apc.indiv.ftest.TS(existing.small.model.fit = fullAPCfit, existing.big.model.fit = TSfit) first.line <- get.table.line(ftest1 = ftest.full.vs.TS, first=TRUE) fit.tab[1, 1:8] <- first.line # Second line: full APC model second.line <- get.table.line(ftest1 = ftest.full.vs.TS, second=TRUE) fit.tab[2, 1:8] <- second.line # Remaining submodels for (i in 3:length(model.design.list)){ againstTS <- apc.indiv.ftest.TS(data, small.model = model.design.list[i], dep.var = dep.var,covariates = covariates, unit = unit, n.coh.excl.start = n.coh.excl.start, n.coh.excl.end = n.coh.excl.end, n.age.excl.start = n.age.excl.start, n.age.excl.end = n.age.excl.end, n.per.excl.start = n.per.excl.start, n.per.excl.end = n.per.excl.end, existing.big.model.fit = TSfit, existing.collinear = fullAPCcollinear) againstAPC <- apc.indiv.ftest.fullapc(data, small.model = model.design.list[i], dep.var=dep.var, covariates=covariates, unit = unit, n.coh.excl.start = n.coh.excl.start, n.coh.excl.end = n.coh.excl.end, n.age.excl.start = n.age.excl.start, n.age.excl.end = n.age.excl.end, n.per.excl.start = n.per.excl.start, n.per.excl.end = n.per.excl.end, existing.big.model.fit = fullAPCfit, existing.collinear = fullAPCcollinear) fit.tab[i, 1:8] <- get.table.line(againstTS, againstAPC) } # Change table column and row names df.againstTS <- paste("DF ( * , ", as.character(againstTS$df.denom), ")", sep="") df.againstAPC <- paste("DF ( * , ", as.character(againstAPC$df.denom), ")", sep="") table.colnames <- c("F-test vs TS", df.againstTS, "p-value", "F-test vs APC", df.againstAPC, "p-value", "AIC", "lik") colnames(fit.tab) <- table.colnames rownames(fit.tab) <- model.design.list ###################### valuables <- list(table = fit.tab, NR.report = NULL) return(valuables) } #### write apc.indiv.logit.TS #### apc.indiv.logit.TS <- function(data, dep.var, covariates=NULL, maxit.loop=10, maxit.linesearch=30, tolerance=.002, init="ols", inv.tol=NULL, d1.tol=.002, custom.kappa=NULL, custom.zeta=NULL){ ################## # This function breaks down the logit estimator step-by-step to exploit the # identity structure of the TS model. Y is the dependent variable, mu the # linear predictor of form mu = kappa*T + zeta*Z where T are the time # dummies and Z are the covariates. ################## if(is.null(inv.tol)){ inv.tol <- .Machine$double.eps } if(!is.null(covariates)){ #### Calculate things that will be used repeatedly #### ## get Y and Z Y <- as.matrix(data[dep.var]) Z <- as.matrix(data[names(data) %in% covariates]) ## get cell name # add cell.name to data: age-cohort cell identifier (if not already there) if (!"cell.name" %in% colnames(data)) { cell.name <- paste("ik", as.character(data$age), as.character(data$cohort), sep="_") data <- cbind(cell.name, data) } else { warning("Variable cell.name must be unique identifier of age-cohort cells; if preexisting variable cell.name has alternate meaning please rename") } ## get first part of first derivative cell.sum.Y.withname <- ddply(data, .variables=c("cell.name"), function(subdata, colnames){colSums(subdata[colnames])}, dep.var) colnames(cell.sum.Y.withname)[2:ncol(cell.sum.Y.withname)] <- paste("cell.sum", dep.var, sep=".") cell.sum.Y <- cell.sum.Y.withname[2:ncol(cell.sum.Y.withname)] d1e1_top <- as.matrix(cell.sum.Y) d1e1_bot <- t(Z) %*% Y #### Get initial starting values #### if((!is.null(custom.kappa)|!is.null(custom.zeta)) & !(init=="custom")) warning("custom.kappa and custom.zeta are not used without init=custom") if (init == "ols"){ # use estimate.TS once to get estimated kappa and zeta initial.TS <- apc.indiv.estimate.TS(data, dep.var = dep.var, covariates = covariates) kappa <- as.matrix(initial.TS$kappahat) zeta <- as.matrix(initial.TS$zetahat) } else if (init == "zero"){ zeta <- rep(0, ncol(Z)) kappa <- rep(0, nrow(cell.sum.Y)) } else if (init == "custom"){ zeta <- custom.zeta kappa <- custom.kappa } else stop("init must be one of 'ols', 'zero', 'custom'") # set loop iteration counter i <- 1 # get mu, the individual predictor kappa.withname <- cbind(cell.sum.Y.withname[1], kappa) indiv.cell.name <- data["cell.name"] indiv.kappa <- join(indiv.cell.name, kappa.withname, by="cell.name") mu <- indiv.kappa$kappa + Z %*% zeta colnames(mu) <- "mu" # evaluate the log-likelihood at the estimated values of kappa and zeta: lik <- (t(Y) %*% mu) - sum(log(1+exp(mu))) # -------------------------- #### The loop #### while (i <= maxit.loop){ ## Step 1: Elements of updating ## # get Nx1 vector of probabilities pi <- exp(mu)/(1+exp(mu)) colnames(pi) <- "pi" # get matrix W (diagonal matrix) w.element <- pi*(1-pi) colnames(w.element) <- "w.element" ## Step 2: Derivatives ## ## First # second part of first derivative data.w.pi <- cbind(data, pi, w.element) cell.sum.pi.withname <- ddply(data.w.pi, .variables=c("cell.name"), function(subdata, colnames){colSums(subdata[colnames])}, "pi") colnames(cell.sum.pi.withname)[2:ncol(cell.sum.pi.withname)] <- paste("cell.sum", "pi", sep=".") cell.sum.pi <- cell.sum.pi.withname[2:ncol(cell.sum.pi.withname)] d1e2_top <- as.matrix(cell.sum.pi) d1e2_bot <- t(Z) %*% pi # complete first derivative d1_top <- d1e1_top - d1e2_top d1_bot <- d1e1_bot - d1e2_bot ## Second # top-left element of second derivative (diagonal matrix) cell.sum.Welm <- ddply(data.w.pi, .variables = c("cell.name"), function(dfr, colnm){sum(dfr[, colnm])}, "w.element") colnames(cell.sum.Welm)[2] <- paste("cell.sum", "w.element", sep=".") TWT.diag <- cell.sum.Welm[,2] # top-right element of second derivative WelmZ <- t(t(Z) * as.vector(w.element)) colnames(WelmZ) <- paste("weighted", colnames(Z), sep=".") data.pi.Z <- cbind(data.w.pi, WelmZ) cell.sums <- as.data.frame(cell.sum.Welm) for (z in colnames(WelmZ)){ cell.sum.Wcov <- ddply(data.pi.Z, .variables = c("cell.name"), function(dfr, colnm){sum(dfr[, colnm])}, z) colnames(cell.sum.Wcov)[2] <- paste("cell.sum", z, sep=".") cell.sums <- join(cell.sums, cell.sum.Wcov, by="cell.name") } tr2 <- as.matrix(cell.sums[, 3:ncol(cell.sums)]) # bottom-left element of second derivative bl2 <- t(tr2) # bottom-right element of second derivative br2 <- t(WelmZ) %*% Z ## Step 3: Inversion of Second derivative ## TWTinv.diag <- 1/TWT.diag inv.tl2 <- diag(TWTinv.diag) # Schur complement schur <- br2 - bl2 %*% inv.tl2 %*% tr2 inv.schur <- solve(schur, tol=inv.tol) # top-left element of second derivative inverse tl2inv <- inv.tl2 + inv.tl2 %*% tr2 %*% inv.schur %*% bl2 %*% inv.tl2 # top-right element of second derivative inverse tr2inv <- - inv.tl2 %*% tr2 %*% inv.schur # bottom-left element of second derivative inverse bl2inv <- - inv.schur %*% bl2 %*% inv.tl2 # bottom-right element of second derivative inverse br2inv <- inv.schur ## Step 4: Update, including linesearch ## ## Starting values # set linesearch parameter linesearch <- 1 # set linesearch iteration counter j <- 1 ## Loop over linesearch while (j <= maxit.linesearch){ # Kappa (TS parameter) kappa_update <- tl2inv %*% d1_top + tr2inv %*% d1_bot kappa_new <- kappa + linesearch*kappa_update colnames(kappa_new) <- "kappa" # zeta (covariate parameter) zeta_update <- bl2inv %*% d1_top + br2inv %*% d1_bot zeta_new <- zeta + linesearch*zeta_update # mu (individual predictor) kappa_new.withname <- cbind(cell.sum.Y.withname[1], kappa_new) indiv.kappa_new <- join(indiv.cell.name, kappa_new.withname, by="cell.name") mu_new <- indiv.kappa_new$kappa + Z %*% zeta_new # evaluate the likelihood at new values lik_new <- (t(Y) %*% mu_new) - sum(log(1+exp(mu_new))) # compare likelihoods: only to see if need to re-enter linesearch if (lik_new - lik < 0){ # New likelihood is lower, have overstepped. Enter linesearch linesearch <- linesearch/2 j <- j+1 } else break } ## Update pi and first derivative # get Nx1 vector of probabilities pi_new <- exp(mu_new)/(1+exp(mu_new)) colnames(pi_new) <- "pi" # get matrix W (diagonal matrix) w.element_new <- pi_new*(1-pi_new) colnames(w.element_new) <- "w.element" # get new second element of first derivative data.w.pi_new <- cbind(data, pi_new) cell.sum.pi_new.withname <- ddply(data.w.pi_new, .variables=c("cell.name"), function(subdata, colnames){colSums(subdata[colnames])}, "pi") colnames(cell.sum.pi_new.withname)[2:ncol(cell.sum.pi_new.withname)] <- paste("cell.sum", "pi_new", sep=".") cell.sum.pi_new <- cell.sum.pi_new.withname[2:ncol(cell.sum.pi_new.withname)] d1e2_top_new <- as.matrix(cell.sum.pi_new) d1e2_bot_new <- t(Z) %*% pi_new # complete first derivative d1_top_new <- d1e1_top - d1e2_top_new d1_bot_new <- d1e1_bot - d1e2_bot_new d1_new <- c(d1_top_new, d1_bot_new) norm.d1 <- norm(as.matrix(d1_new), type="F") ## Step 5: Compare results of updating ## # Case 1: difference exceeds tolerance, update and re-enter loop if (lik_new - lik > tolerance){ kappa <- kappa_new zeta <- zeta_new mu <- mu_new lik <- lik_new i <- i+1 rm(kappa_new, zeta_new, mu_new, lik_new, pi_new, w.element_new, d1e2_top_new, d1e2_bot_new, d1_top_new, d1_bot_new) result <- "exceed lik tolerance, re-enter loop" r2 <- "at likelihood" } else if (lik_new - lik < 0){ # Case 2: old likelihood higher, linesearch unsuccessful, terminate print(paste("overstepped and linesearch iteration limit reached, exiting after", i-1, "iterations", sep=" ")) i <- maxit.loop+4 # sets i too high for loop to continue result <- "overstep" } else if (lik_new - lik <= tolerance & lik_new - lik >= 0){ # Case 3: likelihood condition satisfied, check first derivative # Case 3a: first derivative condition not satisfied, re-enter loop if (norm.d1 > d1.tol){ kappa <- kappa_new zeta <- zeta_new mu <- mu_new lik <- lik_new i <- i+1 rm(kappa_new, zeta_new, mu_new, lik_new, pi_new, w.element_new, d1e2_top_new, d1e2_bot_new, d1_top_new, d1_bot_new) result <- "exceed d1 tolerance, re-enter loop" r2 <- "at first derivative" } # Case 3b: first derivative condition satisfied, terminate if (norm.d1 <= d1.tol){ print(paste("converged after", i, "iterations", sep=" ")) i <- maxit.loop+4 # sets i too high for loop to continue result <- "converge" } } # end of loop } if (i==maxit.loop+1) print(paste("max iterations exceeded, did not converge", r2, sep=" ")) # at this point one ought to use the final values as the estimates and # also construct the SE based on these, but in the interests of expedience # I don't variance <- c(TWT.diag, diag(br2)) Std.error <- sqrt(variance) coefficients.TS <- matrix(nrow=nrow(kappa), ncol=4) coefficients.TS[,1] <- kappa coefficients.TS[,2] <- Std.error[1:nrow(kappa)] coefficients.TS[,3] <- coefficients.TS[,1]/coefficients.TS[,2] coefficients.TS[,4] <- 2*pnorm(abs(coefficients.TS[,3]), lower.tail=FALSE) colnames(coefficients.TS) <- c("Estimate", "Std. Error", "z value", "Pr(>|z|)") rownames(coefficients.TS) <- cell.sum.Welm[,1] coefficients.covariates <- matrix(nrow=nrow(zeta), ncol=4) coefficients.covariates[,1] <- zeta coefficients.covariates[,2] <- Std.error[(nrow(kappa)+1):nrow(as.matrix(Std.error))] coefficients.covariates[,3] <- coefficients.covariates[,1]/coefficients.covariates[,2] coefficients.covariates[,4] <- 2*pnorm(abs(coefficients.covariates[,3]), lower.tail=FALSE) colnames(coefficients.covariates) <- c("Estimate", "Std. Error", "z value", "Pr(>|z|)") rownames(coefficients.covariates) <- rownames(zeta) } #### Case without covariates #### else{ #### Calculate things that will be used repeatedly #### ## get Y and Z Y <- as.matrix(data[dep.var]) ## get first part of first derivative # add cell.name to data: age-cohort cell identifier (if not already there) if (!"cell.name" %in% colnames(data)) { cell.name <- paste("ik", as.character(data$age), as.character(data$cohort), sep="_") data <- cbind(cell.name, data) } else { warning("Variable cell.name must be unique identifier of age-cohort cells; if preexisting variable cell.name has alternate meaning please rename") } cell.sum.Y.withname <- ddply(data, .variables=c("cell.name"), function(subdata, colnames){colSums(subdata[colnames])}, dep.var) colnames(cell.sum.Y.withname)[2:ncol(cell.sum.Y.withname)] <- paste("cell.sum", dep.var, sep=".") cell.sum.Y <- cell.sum.Y.withname[2:ncol(cell.sum.Y.withname)] d1e1 <- as.matrix(cell.sum.Y) #### Get initial starting values #### if((!is.null(custom.kappa)|!is.null(custom.zeta)) & !(init=="custom")) warning("custom.kappa and custom.zeta are not used without init=custom") if (init == "ols"){ # use estimate.TS once to get estimated kappa and zeta initial.TS <- apc.indiv.estimate.TS(data, dep.var = dep.var) kappa <- as.matrix(initial.TS$kappahat) } else if (init == "zero"){ kappa <- rep(0, nrow(cell.sum.Y)) } else if (init == "custom"){ kappa <- custom.kappa } else stop("init must be one of 'ols', 'zero', 'custom'") # set loop iteration counter i <- 1 # get mu (individual predictor) kappa.withname <- cbind(cell.sum.Y.withname[1], kappa) indiv.cell.name <- data["cell.name"] indiv.kappa <- join(indiv.cell.name, kappa.withname, by="cell.name") mu <- as.matrix(indiv.kappa$kappa) colnames(mu) <- "mu" # evaluate the log-likelihood at the estimated values of kappa and zeta: lik <- (t(Y) %*% mu) - sum(log(1+exp(mu))) # -------------------------- #### The loop #### while (i <= maxit.loop){ ## Step 1: Elements of updating ## # get Nx1 vector of probabilities pi <- exp(mu)/(1+exp(mu)) colnames(pi) <- "pi" # get matrix W (diagonal matrix) w.element <- pi*(1-pi) colnames(w.element) <- "w.element" ## Step 2: Derivatives ## ## First data.w.pi <- cbind(data, pi, w.element) cell.sum.pi.withname <- ddply(data.w.pi, .variables=c("cell.name"), function(subdata, colnames){colSums(subdata[colnames])}, "pi") colnames(cell.sum.pi.withname)[2:ncol(cell.sum.pi.withname)] <- paste("cell.sum", "pi", sep=".") cell.sum.pi <- cell.sum.pi.withname[2:ncol(cell.sum.pi.withname)] d1e2 <- as.matrix(cell.sum.pi) # complete first derivative d1 <- d1e1 - d1e2 ## Second # this is the "top-left" element in the covariate case cell.sum.Welm <- ddply(data.w.pi, .variables = c("cell.name"), function(dfr, colnm){sum(dfr[, colnm])}, "w.element") colnames(cell.sum.Welm)[2] <- paste("cell.sum", "w.element", sep=".") d2.diag <- cell.sum.Welm[, 2] ## Step 3: Inversion of Second derivative ## # Invert top-left of second derivative d2inv.diag <- 1/d2.diag inv.d2 <- diag(d2inv.diag) ## Step 4: Update, including linesearch ## ## Starting values # set linesearch parameter linesearch <- 1 # set linesearch iteration counter j <- 1 ## Loop over linesearch while (j <= maxit.linesearch){ # Kappa (TS parameter) kappa_update <- inv.d2 %*% d1 kappa_new <- kappa + linesearch*kappa_update colnames(kappa_new) <- "kappa" # get mu (individual predictor) kappa_new.withname <- cbind(cell.sum.Y.withname[1], kappa_new) indiv.kappa_new <- join(indiv.cell.name, kappa_new.withname, by="cell.name") mu_new <- as.matrix(indiv.kappa_new$kappa) # evaluate the likelihood at new values lik_new <- (t(Y) %*% mu_new) - sum(log(1+exp(mu_new))) # compare likelihoods: only to see if need to re-enter linesearch if (lik_new - lik < 0){ # New likelihood is lower, have overstepped. Enter linesearch linesearch <- linesearch/2 j <- j+1 } else break } ## Update pi and first derivative # get Nx1 vector of probabilities pi_new <- exp(mu_new)/(1+exp(mu_new)) colnames(pi_new) <- "pi" # get matrix W (diagonal matrix) w.element_new <- pi_new*(1-pi_new) colnames(w.element_new) <- "w.element" # get first derivative data.w.pi_new <- cbind(data, pi_new) cell.sum.pi_new.withname <- ddply(data.w.pi_new, .variables=c("cell.name"), function(subdata, colnames){colSums(subdata[colnames])}, "pi") colnames(cell.sum.pi_new.withname)[2:ncol(cell.sum.pi_new.withname)] <- paste("cell.sum", "pi_new", sep=".") cell.sum.pi_new <- cell.sum.pi_new.withname[2:ncol(cell.sum.pi_new.withname)] d1e2_new <- as.matrix(cell.sum.pi_new) # complete first derivative d1_new <- d1e1 - d1e2_new norm.d1 <- norm(as.matrix(d1_new), type="F") ## Step 5: Compare results of updating ## # Case 1: difference exceeds tolerance, update and re-enter loop if (lik_new - lik > tolerance){ kappa <- kappa_new mu <- mu_new lik <- lik_new i <- i+1 rm(kappa_new, mu_new, lik_new, pi_new, w.element_new, d1e2_new) result <- "exceed lik tolerance, re-enter loop" r2 <- "at likelihood" } else if (lik_new - lik < 0){ # Case 2: old likelihood higher, linesearch unsuccessful, terminate print(paste("overstepped and linesearch iteration limit reached, exiting after", i-1, "iterations", sep=" ")) i <- maxit.loop+4 # sets i too high for loop to continue result <- "overstep" } else if (lik_new - lik <= tolerance & lik_new - lik >= 0){ # Case 3: likelihood condition satisfied, check first derivative # Case 3a: first derivative condition not satisfied, re-enter loop if (norm.d1 > d1.tol){ kappa <- kappa_new mu <- mu_new lik <- lik_new i <- i+1 rm(kappa_new, mu_new, lik_new, pi_new, w.element_new, d1e2_new) result <- "exceed d1 tolerance, re-enter loop" r2 <- "at first derivative" } # Case 3b: first derivative condition satisfied, terminate if (norm.d1 <= d1.tol){ print(paste("converged after", i, "iterations", sep=" ")) i <- maxit.loop+4 # sets i too high for loop to continue result <- "converge" } } } if (i==maxit.loop+1) print(paste("max iterations exceeded, did not converge", r2, sep=" ")) # at this point one ought to use the final values as the estimates and # also construct the SE based on these, but in the interests of expedience # I don't variance <- d2inv.diag Std.error <- sqrt(variance) zeta <- NULL coefficients.TS <- matrix(nrow=nrow(kappa), ncol=4) coefficients.TS[,1] <- kappa coefficients.TS[,2] <- Std.error[1:nrow(kappa)] coefficients.TS[,3] <- coefficients.TS[,1]/coefficients.TS[,2] coefficients.TS[,4] <- 2*pnorm(abs(coefficients.TS[,3]), lower.tail=FALSE) colnames(coefficients.TS) <- c("Estimate", "Std. Error", "z value", "Pr(>|z|)") rownames(coefficients.TS) <- cell.sum.Welm[,1] coefficients.covariates <- NULL } # -------------------------- #### Return objects of interest #### #Define param.dimension <- length(kappa)+length(zeta) aic <- -2*lik + 2*param.dimension #Deleted values if(!exists("lik_new")) lik_new <- NULL if(!exists("kappa_new")) kappa_new <- NULL if(!exists("zeta_new")) zeta_new <- NULL valuables <- list(likelihood = lik, kappa = kappa, zeta = zeta, n.loop.iterations = i-1, n.linesearch.iterations = j-1, new.likelihood = lik_new, new.kappa = kappa_new, new.zeta = zeta_new, param.dimension = param.dimension, aic = aic, std.error = Std.error, coefficients.TS = coefficients.TS, coefficients.covariates = coefficients.covariates, d1_new = d1_new, norm.d1 = norm.d1, result = result) return (valuables) } #### write apc.indiv.LRtest.TS #### apc.indiv.LRtest.TS <- function(data, small.model="APC", dep.var, covariates=NULL, model.family="binomial", unit=1, n.coh.excl.start=0, n.coh.excl.end=0, n.age.excl.start=0, n.age.excl.end=0, n.per.excl.start=0, n.per.excl.end=0, existing.small.model.fit=NULL, existing.big.model.fit=NULL, existing.collinear=NULL, maxit.loop.TS=10, maxit.linesearch.TS=30, init.TS="ols"){ # step 1: fit APC model if(is.null(existing.small.model.fit)){ if(is.null(existing.collinear)){ collinear <- apc.indiv.design.collinear(data, unit = unit, n.coh.excl.start = n.coh.excl.start, n.coh.excl.end = n.coh.excl.end, n.age.excl.start = n.age.excl.start, n.age.excl.end = n.age.excl.end, n.per.excl.start = n.per.excl.start, n.per.excl.end = n.per.excl.end) } else collinear <- existing.collinear design <- apc.indiv.design.model(collinear, model.design = small.model, dep.var=dep.var, covariates=covariates) APC.model <- apc.indiv.fit.model(design, model.family=model.family) } else APC.model <- existing.small.model.fit lik.APC <- -0.5*APC.model$deviance # step 2: fit TS model if(is.null(existing.big.model.fit)){ TS.model <- apc.indiv.logit.TS(data, dep.var=dep.var, covariates=covariates, maxit.loop = maxit.loop.TS, maxit.linesearch = maxit.linesearch.TS, init = init.TS) } else TS.model <- existing.big.model.fit # step 3: run LR test diff.df <- TS.model$param.dimension - length(APC.model$coefficients) diff.lik <- lik.APC - TS.model$likelihood p.value <- pchisq(-2*diff.lik, df = diff.df, lower.tail = FALSE) aic.apc <- APC.model$aic aic.TS <- TS.model$aic NR.report <- list(result = TS.model$result, n.loop.iterations = TS.model$n.loop.iterations, n.linesearch.iterations = TS.model$n.linesearch.iterations, d1_new = TS.model$d1_new, norm.d1 = TS.model$norm.d1) # step 4: extract valuables valuables <- list(LR = -2*diff.lik, df = diff.df, p.value = p.value, aic.small = aic.apc, aic.big = aic.TS, lik.small = lik.APC, lik.big = TS.model$likelihood, NR.report = NR.report) return(valuables) } #### write apc.indiv.LRtest.fullapc #### apc.indiv.LRtest.fullapc <- function(data, big.model="APC", small.model, dep.var, covariates=NULL, model.family="binomial", unit=1, n.coh.excl.start=0, n.coh.excl.end=0, n.age.excl.start=0, n.age.excl.end=0, n.per.excl.start=0, n.per.excl.end=0, existing.big.model.fit=NULL, existing.small.model.fit=NULL, existing.collinear=NULL){ # step 1: fit reduced model if(is.null(existing.small.model.fit)){ if(is.null(existing.collinear)){ collinear.R <- apc.indiv.design.collinear(data, unit = unit, n.coh.excl.start = n.coh.excl.start, n.coh.excl.end = n.coh.excl.end, n.age.excl.start = n.age.excl.start, n.age.excl.end = n.age.excl.end, n.per.excl.start = n.per.excl.start, n.per.excl.end = n.per.excl.end) } else collinear.R <- existing.collinear design.R <- apc.indiv.design.model(apc.indiv.design.collinear = collinear.R, model.design = small.model, dep.var=dep.var, covariates=covariates) reduced.model <- apc.indiv.fit.model(design.R, model.family=model.family) } else reduced.model <- existing.small.model.fit # step 2: fit APC model if(is.null(existing.big.model.fit)){ # note: can use same collinear design matrix for full model as for submodel design.U <- apc.indiv.design.model(apc.indiv.design.collinear = collinear.R, model.design = big.model, dep.var=dep.var, covariates=covariates) APC.model <- apc.indiv.fit.model(design.U, model.family=model.family) } else APC.model <- existing.big.model.fit # step 3: run LR test lik.reduced <- -0.5*reduced.model$deviance lik.APC <- -0.5*APC.model$deviance diff.df <- length(APC.model$coefficients) - length(reduced.model$coefficients) diff.lik <- lik.reduced - lik.APC p.value <- pchisq(-2*diff.lik, df=diff.df, lower.tail=FALSE) aic.reduced <- reduced.model$aic aic.big <- APC.model$aic # step 4: extract valuables valuables <- list(LR = -2*diff.lik, df = diff.df, p.value = p.value, aic.big = aic.big, aic.small = aic.reduced, lik.big = lik.APC, lik.small = lik.reduced) return(valuables) } #### write apc.indiv.LRtable.TS #### apc.indiv.LRtable.TS <- function(data, dep.var, covariates=NULL, unit=1, n.coh.excl.start=0, n.coh.excl.end=0, n.age.excl.start=0, n.age.excl.end=0, n.per.excl.start=0, n.per.excl.end=0, maxit.loop.TS=10, maxit.linesearch.TS=30, init.TS="ols"){ # Step 1: function to write table lines get.table.line <- function(LRtest.TS, LRtest.APC, first=FALSE, second=FALSE){ if(isTRUE(first)) # first line is TS model so not compared to itself/APC submodel line <- c(NA, NA, NA, NA, NA, NA, round(LRtest.TS$aic.big, 3), round(LRtest.TS$lik.big, 3)) else if(isTRUE(second)) # second line is full APC model so not compared to itself line <- c(round(LRtest.TS$LR, 3), LRtest.TS$df, round(LRtest.TS$p.value, 3), NA, NA, NA, round(LRtest.TS$aic.small, 3), round(LRtest.TS$lik.small, 3)) else line <- c(round(LRtest.TS$LR, 3), LRtest.TS$df, round(LRtest.TS$p.value, 3), round(LRtest.APC$LR, 3), LRtest.APC$df, round(LRtest.APC$p.value, 3), round(LRtest.APC$aic.small, 3), round(LRtest.APC$lik.small, 3)) return(line) } ###################### # Step 2: other setup # Submodel list model.design.list <- c("TS", "APC", "AP", "AC", "PC", "Ad", "Pd", "Cd", "A", "P", "C", "t", "tA", "tP", "tC", "1") # Empty table fit.tab <- matrix(nrow=length(model.design.list), ncol=8, data=NA) # Full APC model for comparison fullAPCcollinear <- apc.indiv.design.collinear(data, unit=unit, n.coh.excl.start = n.coh.excl.start, n.coh.excl.end = n.coh.excl.end, n.age.excl.start = n.age.excl.start, n.age.excl.end = n.age.excl.end, n.per.excl.start = n.per.excl.start, n.per.excl.end = n.per.excl.end) fullAPCdesign <- apc.indiv.design.model(fullAPCcollinear, dep.var = dep.var, covariates = covariates) fullAPCfit <- apc.indiv.fit.model(fullAPCdesign, model.family="binomial") # TS model for comparison TSfit <- apc.indiv.logit.TS(data, dep.var=dep.var, covariates=covariates, maxit.loop = maxit.loop.TS, maxit.linesearch = maxit.linesearch.TS, init = init.TS) ###################### # Step 3: Generate table # first line, TS model and second line, full APC model LRtest.full.vs.TS <- apc.indiv.LRtest.TS(existing.small.model.fit = fullAPCfit, existing.big.model.fit = TSfit) first.line <- get.table.line(LRtest.TS=LRtest.full.vs.TS, first=TRUE) second.line <- get.table.line(LRtest.TS=LRtest.full.vs.TS, second=TRUE) fit.tab[1, 1:8] <- first.line fit.tab[2, 1:8] <- second.line # Remaining submodels for (i in 3:length(model.design.list)){ againstTS <- apc.indiv.LRtest.TS(data, small.model=model.design.list[i], dep.var = dep.var,covariates = covariates, unit = unit, n.coh.excl.start = n.coh.excl.start, n.coh.excl.end = n.coh.excl.end, n.age.excl.start = n.age.excl.start, n.age.excl.end = n.age.excl.end, n.per.excl.start = n.per.excl.start, n.per.excl.end = n.per.excl.end, existing.big.model.fit = TSfit, existing.collinear = fullAPCcollinear) againstAPC <- apc.indiv.LRtest.fullapc(data, small.model = model.design.list[i], dep.var=dep.var, covariates=covariates, unit = unit, n.coh.excl.start = n.coh.excl.start, n.coh.excl.end = n.coh.excl.end, n.age.excl.start = n.age.excl.start, n.age.excl.end = n.age.excl.end, n.per.excl.start = n.per.excl.start, n.per.excl.end = n.per.excl.end, existing.big.model.fit = fullAPCfit, existing.collinear = fullAPCcollinear) fit.tab[i, 1:8] <- get.table.line(againstTS, againstAPC) } # Change table column and row names table.colnames <- c("LR-test vs TS", "df", "p-value", "LR-test vs APC", "df", "p-value", "AIC", "Loglihood") colnames(fit.tab) <- table.colnames rownames(fit.tab) <- model.design.list #################### valuables <- list(table = fit.tab, NR.report = againstTS$NR.report) return(valuables) } # ----------------------------------------------------------------------------- # -----------------------------------------------------------------------------