Skip to content

Commit 42bb6c8

Browse files
msuchardmsuchard
committed
first working test of proportionality with strata and weights
1 parent 4958256 commit 42bb6c8

4 files changed

Lines changed: 106 additions & 97 deletions

File tree

‎R/Residuals.R‎

Lines changed: 19 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -72,11 +72,10 @@ residuals.cyclopsFit <- function(object, parm = NULL, type = "schoenfeld", ...)
7272
#' @param object A Cyclops model fit object
7373
#' @param parm A specification of which parameters require a proportionality test,
7474
#' either a vector of numbers or covariateId names
75-
#' @param transformedTimes Vector of transformed time. If NULL, then the default "identity"
76-
#' transform of `cox.zph` is used
75+
#' @param transform Same as `transform` in `cox.zph`
7776
#'
7877
#' @export
79-
testProportionality <- function(object, parm = NULL, transformedTimes = NULL) {
78+
testProportionality <- function(object, parm = NULL, transform = "km") {
8079

8180
.checkInterface(object$cyclopsData, testOnly = TRUE)
8281

@@ -89,17 +88,26 @@ testProportionality <- function(object, parm = NULL, transformedTimes = NULL) {
8988
stop("Only single-covariate models are currently implemented")
9089
}
9190

92-
if (is.null(transformedTimes)) {
93-
times <- getTimeVector(object$cyclopsData)
94-
y <- getYVector(object$cyclopsData)
95-
transformedTimes <- times - mean(times[y == 1])
91+
times <- getTimeVector(object$cyclopsData)
92+
y <- getYVector(object$cyclopsData)
93+
survY <- survival::Surv(time = times, event = y)
94+
if (is.character(transform)) {
95+
tname <- transform
96+
ttimes <- switch(transform, identity = times, rank = rank(times),
97+
log = log(times), km = {
98+
temp <- survfitKM(factor(rep(1L, nrow(survY))), survY,
99+
se.fit = FALSE)
100+
indx <- findInterval(times, temp$time, left.open = TRUE)
101+
1 - c(1, temp$surv)[indx + 1]
102+
}, stop("Unrecognized transform"))
96103
} else {
97-
if (getNumberOfRows(object$cyclopsData) != length(transformedTimes)) {
98-
stop("Incorrect 'transformedTime' length")
104+
tname <- deparse(substitute(transform))
105+
if (length(tname) > 1) {
106+
tname <- "user"
99107
}
100-
101-
transformedTimes <- transformedTimes[object$cyclopsData$sortOrder]
108+
ttimes <- transform(times)
102109
}
110+
transformedTimes <- ttimes - mean(ttimes[y == 1])
103111

104112
res <- .cyclopsTestProportionality(object$interface, NULL, transformedTimes)
105113
nCovariates <- 1 # TODO Remove

‎man/testProportionality.Rd‎

Lines changed: 2 additions & 3 deletions
Some generated files are not rendered by default. Learn more about customizing how changed files appear on GitHub.

‎src/cyclops/engine/ModelSpecifics.hpp‎

Lines changed: 29 additions & 74 deletions
Original file line numberDiff line numberDiff line change
@@ -774,7 +774,7 @@ void ModelSpecifics<BaseModel,RealType>::getSchoenfeldResidualsImpl(int index,
774774
IteratorType it(hX, index);
775775

776776
RealType resNumerator = static_cast<RealType>(0);
777-
RealType resDenominator = static_cast<RealType>(0);
777+
// RealType resDenominator = static_cast<RealType>(0);
778778
RealType scoreNumerator1 = static_cast<RealType>(0);
779779
RealType scoreNumerator2 = static_cast<RealType>(0);
780780
RealType scoreDenominator = static_cast<RealType>(0);
@@ -793,16 +793,10 @@ void ModelSpecifics<BaseModel,RealType>::getSchoenfeldResidualsImpl(int index,
793793

794794
auto processRow = [&](int i, RealType x) {
795795

796-
// std::cerr << "row " << i << "\n";
797-
798-
// std::cerr << "r" << *reset << " ";
799-
800-
if (*reset <= i) {
801-
802-
// std::cerr << *reset << " ";
796+
if (*reset <= hPid[i]) {
803797

804798
resNumerator = static_cast<RealType>(0);
805-
resDenominator = static_cast<RealType>(0);
799+
// resDenominator = static_cast<RealType>(0);
806800
scoreNumerator1 = static_cast<RealType>(0);
807801
scoreNumerator2 = static_cast<RealType>(0);
808802
scoreDenominator = static_cast<RealType>(0);
@@ -812,25 +806,27 @@ void ModelSpecifics<BaseModel,RealType>::getSchoenfeldResidualsImpl(int index,
812806

813807
const auto expXBeta = offsExpXBeta[i]; // std::exp(hXBeta[i]);
814808

815-
resNumerator += expXBeta * x;
816-
resDenominator += expXBeta;
809+
resNumerator += hKWeight[i] * expXBeta * x;
810+
// resDenominator += expXBeta;
817811

818812
if (hY[i] == 1) {
819813
// std::cerr << " " << x << " for " << resNumerator << " / " << resDenominator << " - " << ( x - resNumerator / resDenominator) << "\n";
820814
// std::cerr << denomPid[i] << "\n";
821815
if (hasResiduals) {
822-
if (i < N && i < *reset &&
816+
if (i < (K - 1) && hPid[i] < *reset &&
823817
(hOffs[i] == hOffs[i+ 1]) && hY[i + 1] == 1) {
824818
deferredResidualX.push_back(x);
825819
} else {
826820

821+
const auto denom = accDenomPid[hPid[i]];
822+
827823
while (deferredResidualX.size() > 0) {
828-
const auto residual = deferredResidualX.front() - resNumerator / resDenominator;
824+
const auto residual = deferredResidualX.front() - resNumerator / denom;
829825
residuals->push_back(residual);
830826
deferredResidualX.pop_front();
831827
}
832828

833-
const auto residual = x - resNumerator / resDenominator;
829+
const auto residual = x - resNumerator / denom;
834830
residuals->push_back(residual);
835831
}
836832
}
@@ -845,64 +841,21 @@ void ModelSpecifics<BaseModel,RealType>::getSchoenfeldResidualsImpl(int index,
845841
}
846842

847843
if (hasScore) {
848-
849-
// std::cerr << "i = " << i << "\n";
850-
const auto weight = covariate[i];
851-
// MAS does not believe reweighing is correct, but is matching cox.zph
852-
853-
// const auto xt = x * covariate[i];
854-
// MAS believes covariate should be adjusted
855-
const auto numerator1 = expXBeta * x; // TODO not xt?
856-
const auto numerator2 = expXBeta * x * x; // TODO not xt?
844+
const auto cov = covariate[i];
845+
const auto numerator1 = hKWeight[i] * expXBeta * x;
846+
const auto numerator2 = hKWeight[i] * expXBeta * x * x;
857847

858848
if (hY[i] == 1) {
859-
uGradient2 += x;
860-
wGradient2 += x * weight; // TODO not xt and no weight?
849+
uGradient2 += hKWeight[i] * x;
850+
wGradient2 += hKWeight[i] * x * cov;
861851
}
862852

863853
scoreNumerator1 += numerator1;
864854
scoreNumerator2 += numerator2;
865855

866-
// if (hY[i] == 1) {
867-
868-
if (hY[i] == 1 && i < (K - 1) && i < *reset &&
869-
(hOffs[i] == hOffs[i+ 1]) && hY[i + 1] == 1) {
870-
871-
// deferredScore1.push_back(scoreNumerator1);
872-
// deferredScore2.push_back(scoreNumerator2);
873-
// // deferredScore3.push_back(hNWeight[hPid[i]]);
874-
// deferredScore3.push_back(0.0);
875-
// deferredScore4.push_back(i);
876-
877-
} else {
878-
879-
// while (deferredScore1.size() > 0) {
880-
//
881-
// const auto popScoreNumerator1 = deferredScore1.front();
882-
// const auto popScoreNumerator2 = deferredScore2.front();
883-
// const auto popNWeight = deferredScore3.front();
884-
// const auto popI = deferredScore4.front();
885-
//
886-
// const auto denom = accDenomPid[hPid[popI]];
887-
//
888-
// const auto t = scoreNumerator1 / denom;
889-
// const auto gradient = popNWeight * t;
890-
// const auto hessian = popNWeight * (scoreNumerator2 / denom - t * t);
891-
//
892-
// uGradient += gradient;
893-
// wGradient += gradient * weight;
894-
//
895-
// // std::cerr<< "newP" << popI << " " << hPid[popI] << " " << scoreNumerator1 << "/" << denom << " " << hY[popI] << " " << uGradient << " " << popNWeight << "\n";
896-
//
897-
// uHessian += hessian;
898-
// wHessian += hessian * weight * weight;
899-
// xHessian += hessian * weight;
900-
//
901-
// deferredScore1.pop_front();
902-
// deferredScore2.pop_front();
903-
// deferredScore3.pop_front();
904-
// deferredScore4.pop_front();
905-
// }
856+
if (i < (K - 1) && hPid[i] == hPid[i + 1]) {
857+
// zero contribution under Breslow
858+
} else {
906859

907860
const auto denom = accDenomPid[hPid[i]];
908861

@@ -911,22 +864,19 @@ void ModelSpecifics<BaseModel,RealType>::getSchoenfeldResidualsImpl(int index,
911864
const auto hessian = hNWeight[hPid[i]] * (scoreNumerator2 / denom - t * t);
912865

913866
uGradient += gradient;
914-
wGradient += gradient * weight;
867+
wGradient += gradient * cov;
915868

916-
// std::cerr<< "newA" << i << " " << hPid[i] << " " << scoreNumerator1 << "/" << denom << " " << hY[i] << " " << uGradient << " " << hNWeight[hPid[i]] << "\n";
869+
// std::cerr << "newA" << i << " " << hPid[i] << " " << scoreNumerator1 << "/" << denom << " " << hY[i] << " " << uGradient << " " << hNWeight[hPid[i]] << "\n";
917870

918871
uHessian += hessian;
919-
wHessian += hessian * weight * weight;
920-
xHessian += hessian * weight;
872+
wHessian += hessian * cov * cov;
873+
xHessian += hessian * cov;
921874
}
922875
}
923876
};
924877

925878
// main loop
926-
927-
// for (int i = 0; i <= N; ++i) {
928879
for (int i = 0; i < K; ++i) {
929-
// std::cerr << i << " ";
930880
if (i == it.index()) {
931881
processRow(i, it.value());
932882
++it;
@@ -1183,6 +1133,9 @@ void ModelSpecifics<BaseModel,RealType>::computeGradientAndHessianImpl(int index
11831133
w, // Signature-only, for iterator-type specialization
11841134
&gradient, &hessian, accNumerPid, accNumerPid2,
11851135
accDenomPid[i], hNWeight[i], 0.0, hXBeta[i], hY[i]); // When function is in-lined, compiler will only use necessary arguments
1136+
1137+
// std::cerr << "oldA" << i << " " << accNumerPid << "/" << accDenomPid[i] << " " << hY[i] << " " << gradient << " " << hNWeight[i] << "\n";
1138+
11861139
++it;
11871140

11881141
if (IteratorType::isSparse) {
@@ -1201,6 +1154,8 @@ void ModelSpecifics<BaseModel,RealType>::computeGradientAndHessianImpl(int index
12011154
w, // Signature-only, for iterator-type specialization
12021155
&gradient, &hessian, accNumerPid, accNumerPid2,
12031156
accDenomPid[i], hNWeight[i], static_cast<RealType>(0), hXBeta[i], hY[i]); // When function is in-lined, compiler will only use necessary arguments
1157+
1158+
// std::cerr << "oldB" << i << " " << accNumerPid << "/" << accDenomPid[i] << " " << hY[i] << " " << gradient << " " << hNWeight[i] << "\n";
12041159
}
12051160
}
12061161
}
@@ -2312,7 +2267,7 @@ void ModelSpecifics<BaseModel,RealType>::setPidForAccumulationImpl(const AnyReal
23122267

23132268
if (nextPid != lastPid) { // start new strata
23142269
pid++;
2315-
accReset.push_back(pid + nTies);
2270+
accReset.push_back(pid);
23162271
nTies = 0;
23172272
lastPid = nextPid;
23182273
} else {
@@ -2333,7 +2288,7 @@ void ModelSpecifics<BaseModel,RealType>::setPidForAccumulationImpl(const AnyReal
23332288
}
23342289
}
23352290
pid++;
2336-
accReset.push_back(pid + nTies); // TODO BUG HERE????
2291+
accReset.push_back(pid);
23372292

23382293
// Save number of denominators
23392294
N = pid;

‎tests/testthat/test-residuals.R‎

Lines changed: 56 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,6 @@ library("survival")
44
suppressWarnings(RNGversion("3.5.0"))
55

66
test_that("Check Schoenfeld residuals and PH test, no strata", {
7-
skip("residuals not yet implemented")
87
gfit <- coxph(Surv(futime, fustat) ~ age,
98
data=ovarian, method = "breslow")
109
gres <- residuals(gfit, "schoenfeld")
@@ -28,7 +27,6 @@ test_that("Check Schoenfeld residuals and PH test, no strata", {
2827
})
2928

3029
test_that("Check Schoenfeld residuals and PH test, with strata", {
31-
skip("residuals not yet implemented")
3230
gfit <- coxph(Surv(futime, fustat) ~ age + strata(ecog.ps),
3331
data=ovarian, method = "breslow")
3432
gres <- residuals(gfit, "schoenfeld")
@@ -50,8 +48,6 @@ test_that("Check Schoenfeld residuals and PH test, with strata", {
5048
})
5149

5250
test_that("Check Schoenfeld residuals and PH test, with sparse covariates", {
53-
skip("residuals not yet implemented")
54-
5551
test <- read.table(header=T, sep = ",", text = "
5652
start, length, event, x1, x2
5753
0, 4, 1,0.2,0
@@ -105,8 +101,58 @@ start, length, event, x1, x2
105101
# show sch_x1
106102
})
107103

104+
test_that("Check residuals on large Cox regression with weighting",{
105+
tolerance <- 1E-4
106+
set.seed(123)
107+
sim <- simulateCyclopsData(nstrata=1000, #1
108+
ncovars=1,
109+
nrows=10000, #30
110+
effectSizeSd=0.5,
111+
eCovarsPerRow=2,
112+
model="survival")
113+
sim$outcomes$weights <- 1/sim$outcomes$rr
114+
# sim$outcomes$weights <- rep(1, length(sim$outcomes$weights))
115+
sim$outcomes$weights <- runif(length(sim$outcomes$weights))
116+
117+
# Gold standard
118+
covariates <- sim$covariates
119+
ncovars <- max(covariates$covariateId)
120+
nrows <- nrow(sim$outcomes)
121+
m <- matrix(0,nrows,ncovars)
122+
for (i in 1:nrow(covariates)){
123+
m[covariates$rowId[i],covariates$covariateId[i]] <- 1
124+
}
125+
data <- as.data.frame(m)
126+
127+
data$rowId <- 1:nrow(data)
128+
data <- merge(data,sim$outcomes)
129+
data <- data[order(data$stratumId,data$rowId),]
130+
formula <- as.formula(paste(c("Surv(time,y) ~ strata(stratumId)",paste("V",1:ncovars,sep="")),collapse=" + "))
131+
fitCoxph <- survival::coxph(formula, data = data, weights = data$weights, ties = "breslow")
132+
133+
# Cyclops
134+
cyclopsData <- convertToCyclopsData(outcomes = sim$outcomes,
135+
covariates = sim$covariates,
136+
modelType = "cox")
137+
fitCyclops <- fitCyclopsModel(cyclopsData = cyclopsData)
138+
139+
expect_equivalent(coef(fitCyclops), coef(fitCoxph), tolerance = tolerance)
140+
141+
resCoxph <- residuals(fitCoxph, type = "schoenfeld")
142+
resCyclops <- residuals(fitCyclops, type = "schoenfeld")
143+
144+
expect_equivalent(resCoxph, resCyclops)
145+
146+
gtest1 <- cox.zph(fitCoxph, transform = "km", global = FALSE)
147+
ctest1 <- testProportionality(fitCyclops, transform = "km")
148+
expect_equivalent(gtest1$table, ctest1$table)
149+
150+
gtest2 <- cox.zph(fitCoxph, transform = "identity", global = FALSE)
151+
ctest2 <- testProportionality(fitCyclops, transform = "identity")
152+
expect_equivalent(gtest2$table, ctest2$table)
153+
})
154+
108155
test_that("Check Schoenfeld residuals and PH test, with sparse covariates", {
109-
skip("residuals not yet implemented")
110156
test <- read.table(header=T, sep = ",", text = "
111157
start, length, event, x1, x2
112158
0, 4, 1,0,0
@@ -121,11 +167,12 @@ start, length, event, x1, x2
121167
gfit <- coxph(Surv(length, event) ~ x1 + strata(x2), test, ties = "breslow")
122168
gres <- residuals(gfit, "schoenfeld")
123169

124-
data <- createCyclopsData(Surv(length, event) ~ x1+ strata(x2),
125-
# sparseFormula = ~ x1,
170+
data <- createCyclopsData(Surv(length, event) ~ strata(x2),
171+
sparseFormula = ~ x1,
126172
data = test, modelType = "cox")
127173

128174
cfit <- fitCyclopsModel(data)
129-
cres <- residuals(cfit, "schoenfeld") # TODO broken (and not even sparse yet)
130-
# expect_equivalent(cres, gres) # TODO
175+
cres <- residuals(cfit, "schoenfeld")
176+
177+
expect_equivalent(cres, gres)
131178
})

0 commit comments

Comments
 (0)