From 4923a23d905e970fcbb2a5ace741f434bfa03083 Mon Sep 17 00:00:00 2001 From: Birgen Haest Date: Mon, 16 Jan 2017 17:31:06 +0100 Subject: [PATCH 1/4] Added in functionality in slidingwin, randwin, and basewin for differencing of values after the stat has been applied. This is useful for time series to avoid spurious correlations. --- DESCRIPTION | 2 +- R/otherfunctions.R | 37 ++++++++++++++++++++++++++++++++----- R/randwin.R | 5 +++-- R/slidingwin.R | 6 ++++-- 4 files changed, 40 insertions(+), 10 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 421482a..14f3c0e 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,4 +1,4 @@ -Package: climwin +Package: climwinBH Type: Package Title: Climate Window Analysis Version: 1.1.0 diff --git a/R/otherfunctions.R b/R/otherfunctions.R index db8b5f5..6dc53c8 100644 --- a/R/otherfunctions.R +++ b/R/otherfunctions.R @@ -28,7 +28,7 @@ basewin <- function(exclude, xvar, cdate, bdate, baseline, range, type, stat = "mean", func = "lin", refday, cmissing = FALSE, cinterval = "day", nrandom = 0, k = 0, spatial, upper = NA, lower = NA, binary = FALSE, centre = list(NULL, "both"), - cohort = NULL, fast){ + cohort = NULL, fast, cw_diff = FALSE){ print("Initialising, please wait...") @@ -127,7 +127,7 @@ basewin <- function(exclude, xvar, cdate, bdate, baseline, range, modlist <- list() # dataframes to store ouput baseline <- update(baseline, .~.) nullmodel <- AICc(baseline) - modeldat <- model.frame(baseline) + modeldat <- model.frame(baseline) if(attr(baseline, "class")[1] == "lme"){ @@ -158,8 +158,16 @@ basewin <- function(exclude, xvar, cdate, bdate, baseline, range, func <- "centre" } + # BH - probably needs updating here for cw_diff: + tmp = diff(modeldat$yvar) + modeldat = data.frame(yvar = tmp) + # bdate = bdate[2:length(bdate)] --> !!!NOT!!!! + # cont$bintno = cont$bintno[2:length(cont$bintno)] --> NOT!! + ifelse(class(baseline)[length(class(baseline))]=="coxph", leng<-length(modeldat$yvar[,1]), leng<-length(modeldat$yvar)) - if (leng != length(bdate)){ + if ((leng != length(bdate)) && (!cw_diff)){ + stop("NA values present in biological response. Please remove NA values") + } else if ((leng != (length(bdate)-1)) && (cw_diff)){ stop("NA values present in biological response. Please remove NA values") } @@ -229,6 +237,9 @@ basewin <- function(exclude, xvar, cdate, bdate, baseline, range, cmatrix <- as.matrix(cmatrix[, c(ncol(cmatrix):1)]) + # BH insert: delete the final row for cw_diff = TRUE --> NOT!! + # cmatrix = cmatrix[-nrow(cmatrix),] + if(cmissing == FALSE && length(which(is.na(cmatrix))) > 0){ if(is.null(spatial) == FALSE){ @@ -461,8 +472,16 @@ basewin <- function(exclude, xvar, cdate, bdate, baseline, range, time <- seq(1, n, 1) modeldat$climate <- apply(cmatrix[, windowclose:windowopen], 1, FUN = function(x) coef(lm(x ~ time))[2]) } else { - ifelse (n == 1, modeldat$climate <- cmatrix[, windowclose:windowopen], + # Case: NOT differencing + if (!cw_diff){ + ifelse (n == 1, modeldat$climate <- cmatrix[, windowclose:windowopen], modeldat$climate <- apply(cmatrix[, windowclose:windowopen], 1, FUN = stat)) + # Case: Differencing + } else { + ifelse (n == 1, modeldat$climate <- diff(cmatrix[, windowclose:windowopen]), + modeldat$climate <- diff(apply(cmatrix[, windowclose:windowopen], 1, FUN = stat))) + } + } if (min(modeldat$climate) <= 0 & func == "log" || min(modeldat$climate) <= 0 & func == "inv"){ @@ -758,9 +777,17 @@ basewin <- function(exclude, xvar, cdate, bdate, baseline, range, time <- seq(1, n[1], 1) modeldat$climate <- apply(cmatrix[, windowclose:windowopen], 1, FUN = function(x) coef(lm(x ~ time))[2]) } else { - ifelse (windowopen - windowclose == 0, + # CASE: NOT differencing + if (!cw_diff){ + ifelse (windowopen - windowclose == 0, modeldat$climate <- cmatrix[, windowclose:windowopen], modeldat$climate <- apply(cmatrix[, windowclose:windowopen], 1, FUN = stat)) + # CASE: differencing + } else { + ifelse (windowopen - windowclose == 0, + modeldat$climate <- diff(cmatrix[, windowclose:windowopen]), + modeldat$climate <- diff(apply(cmatrix[, windowclose:windowopen], 1, FUN = stat))) + } } if (is.null(centre[[1]]) == FALSE){ diff --git a/R/randwin.R b/R/randwin.R index ff6f1e6..1082fd8 100644 --- a/R/randwin.R +++ b/R/randwin.R @@ -156,7 +156,8 @@ randwin <- function(exclude = NA, repeats = 5, window = "sliding", xvar, cdate, upper = NA, lower = NA, binary = FALSE, centre = list(NULL, "both"), k = 0, weightfunc = "W", par = c(3, 0.2, 0), control = list(ndeps = c(0.01, 0.01, 0.01)), method = "L-BFGS-B", cutoff.day = NULL, cutoff.month = NULL, - furthest = NULL, closest = NULL, thresh = NULL, cvk = NULL){ + furthest = NULL, closest = NULL, thresh = NULL, cvk = NULL, + cw_diff = FALSE){ fast = FALSE @@ -254,7 +255,7 @@ randwin <- function(exclude = NA, repeats = 5, window = "sliding", xvar, cdate, upper = ifelse(binarylevel == "two" || binarylevel == "upper", allcombos$upper[combo], NA), lower = ifelse(binarylevel == "two" || binarylevel == "lower", allcombos$lower[combo], NA), binary = paste(allcombos$binary[combo]), centre = centre, k = k, spatial = spatial, - cohort = cohort, fast = fast) + cohort = cohort, fast = fast, cw_diff = cw_diff) outputrep$Repeat <- r WeightDist <- sum(as.numeric(cumsum(outputrep$ModWeight) <= 0.95))/nrow(outputrep) diff --git a/R/slidingwin.R b/R/slidingwin.R index 59e924f..cf0ec03 100644 --- a/R/slidingwin.R +++ b/R/slidingwin.R @@ -188,7 +188,8 @@ slidingwin <- function(exclude = NA, xvar, cdate, bdate, baseline, spatial = NULL, cohort = NULL, cutoff.day = NULL, cutoff.month = NULL, furthest = NULL, closest = NULL, - thresh = NULL, cvk = NULL){ + thresh = NULL, cvk = NULL, + cw_diff = FALSE){ fast = FALSE @@ -272,6 +273,7 @@ slidingwin <- function(exclude = NA, xvar, cdate, bdate, baseline, print(allcombos) combined <- list() + for (combo in 1:nrow(allcombos)){ runs <- basewin(exclude = exclude, xvar = xvar[[paste(allcombos[combo, 1])]], cdate = cdate, bdate = bdate, baseline = baseline, range = range, type = paste(allcombos[combo, 2]), refday = refday, stat = paste(allcombos[combo, 3]), func = paste(allcombos[combo, 4]), @@ -279,7 +281,7 @@ slidingwin <- function(exclude = NA, xvar, cdate, bdate, baseline, upper = ifelse(binarylevel == "two" || binarylevel == "upper", allcombos$upper[combo], NA), lower = ifelse(binarylevel == "two" || binarylevel == "lower", allcombos$lower[combo], NA), binary = paste(allcombos$binary[combo]), centre = centre, cohort = cohort, - spatial = spatial, fast = fast) + spatial = spatial, fast = fast, cw_diff = cw_diff) combined[[combo]] <- runs allcombos$DeltaAICc[combo] <- round(runs$Dataset$deltaAICc[1], digits = 2) From 4782f101fcff2219974a596da627beeb39fb8196 Mon Sep 17 00:00:00 2001 From: Birgen Haest Date: Mon, 16 Jan 2017 20:04:24 +0100 Subject: [PATCH 2/4] Fixed introduced bug. --- R/otherfunctions.R | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/R/otherfunctions.R b/R/otherfunctions.R index 6dc53c8..0a93a82 100644 --- a/R/otherfunctions.R +++ b/R/otherfunctions.R @@ -158,11 +158,11 @@ basewin <- function(exclude, xvar, cdate, bdate, baseline, range, func <- "centre" } - # BH - probably needs updating here for cw_diff: - tmp = diff(modeldat$yvar) - modeldat = data.frame(yvar = tmp) - # bdate = bdate[2:length(bdate)] --> !!!NOT!!!! - # cont$bintno = cont$bintno[2:length(cont$bintno)] --> NOT!! + # update in case of differencing + if (cw_diff){ + tmp = diff(modeldat$yvar) + modeldat = data.frame(yvar = tmp) + } ifelse(class(baseline)[length(class(baseline))]=="coxph", leng<-length(modeldat$yvar[,1]), leng<-length(modeldat$yvar)) if ((leng != length(bdate)) && (!cw_diff)){ From c37863499cadca8bd63d66f0a09b79f154c81416 Mon Sep 17 00:00:00 2001 From: Birgen Haest Date: Tue, 17 Jan 2017 16:50:54 +0100 Subject: [PATCH 3/4] Deleted some comments of my previous tests. This commit merely deletes 2 lines of comments that were used in previous tests. --- R/otherfunctions.R | 3 --- 1 file changed, 3 deletions(-) diff --git a/R/otherfunctions.R b/R/otherfunctions.R index 0a93a82..8f030b6 100644 --- a/R/otherfunctions.R +++ b/R/otherfunctions.R @@ -237,9 +237,6 @@ basewin <- function(exclude, xvar, cdate, bdate, baseline, range, cmatrix <- as.matrix(cmatrix[, c(ncol(cmatrix):1)]) - # BH insert: delete the final row for cw_diff = TRUE --> NOT!! - # cmatrix = cmatrix[-nrow(cmatrix),] - if(cmissing == FALSE && length(which(is.na(cmatrix))) > 0){ if(is.null(spatial) == FALSE){ From 2e07df1a4b8866a9aefe89f654eb0418ac44b29a Mon Sep 17 00:00:00 2001 From: Birgen Haest Date: Fri, 31 Mar 2017 11:31:26 +0200 Subject: [PATCH 4/4] Corrected nullmodel AICc (wasn't adjusted before to the differenced yvar values). --- R/otherfunctions.R | 2 ++ 1 file changed, 2 insertions(+) diff --git a/R/otherfunctions.R b/R/otherfunctions.R index 8f030b6..bffdb44 100644 --- a/R/otherfunctions.R +++ b/R/otherfunctions.R @@ -162,6 +162,8 @@ basewin <- function(exclude, xvar, cdate, bdate, baseline, range, if (cw_diff){ tmp = diff(modeldat$yvar) modeldat = data.frame(yvar = tmp) + baseline <- update(baseline, formula = yvar~., data = modeldat) + nullmodel <- AICc(baseline) } ifelse(class(baseline)[length(class(baseline))]=="coxph", leng<-length(modeldat$yvar[,1]), leng<-length(modeldat$yvar))