closes #14 closes #17 and passes all the current version tests.

shug0131 · shug0131 · commit f30679360a67 · 2019-11-26T17:00:17.000Z
diff --git a/R/cox.zph.R b/R/cox.zph.R
@@ -0,0 +1,26 @@
+#' Test the proportional hazards assumption of an RPSFTM/Cox Regression
+#' @param fit the result of fitting a rpsftm model using \code{\link[survival]{coxph}} as the inner estimation tool.
+#' @param ... any other arguments to pass to \code{\link[survival]{cox.zph}}.
+#' 
+#' @description If the the \code{fit} inherits *both* \code{rpsftm} and \code{coxph} then this pulls out the genuine 
+#' \code{survival::coxph} object that is deeply nested in the object, and then runs \code{survival::cox.zph} on it.
+#' Or it avoids overwriting the function from \code{survival} by calling \code{survival::cox.zph} directly if the object 
+#' does not inherit \code{rpsftm}. Or it fails.
+#' 
+#' @note This does rely on the order of loading packages. The \code{rpsftm} package must be loaded after \code{survival}, 
+#' if both are required, to avoid the masking of synonymous functions causing errors. 
+#' 
+#' @seealso  \code{\link[survival]{cox.zph}}
+#' @export
+
+
+cox.zph <- function(fit,...){
+  if(inherits(fit, "rpsftm")){
+    if( !inherits(fit, "coxph")){
+      stop("This model did not use coxph within g-estimation")
+    } else {
+      fit <- attr(fit$ans$f.root,"fit")
+    }
+  }
+  survival::cox.zph(fit,...=...)
+}
diff --git a/R/est_eqn.R b/R/est_eqn.R
@@ -24,9 +24,12 @@
 est_eqn <- function(psi, response, data, formula_list, treatment_matrix, rand_matrix,  target = 0, test = "survdiff", 
                    autoswitch, ...) {
   
+  # ensure the psi is now always a matrix, with rows equal to n.
   if ("(treat_modifier)" %in% names(data)) {
-      psi <- psi * data[, "(treat_modifier)"]
-  }
+      psi <-  data[, "(treat_modifier)"] %*%diag(psi, length(psi), length(psi))
+  } else{ 
+    psi <- matrix(psi, nrow=nrow(treatment_matrix), ncol=length(psi), byrow=TRUE)
+    }
   
   #response <- model.response(data)
   #treatment_matrix <- model.matrix(treatment, data=data)
diff --git a/R/print.rpsftm.R b/R/print.rpsftm.R
@@ -43,51 +43,52 @@ print.rpsftm <- function(x,...) {
 #' @keywords internal
 
 print.coxph <- function (x, digits = max(options()$digits - 4, 3), ...){
-    if (!is.null(x$fail)) {
-      cat(" Coxph failed.", x$fail, "\n")
-      return()
-    }
-    savedig <- options(digits = digits)
-    on.exit(options(savedig))
-    coef <- x$coefficients
-    rand_names <- colnames(model.matrix(x$formula_list$randomise, data=x$data))
-    rand_names <- rand_names[rand_names!="(Intercept)"]
-    arm_index <- which(names(coef) %in% rand_names)
-    coef <- coef[-arm_index, drop=FALSE]
-    se <- sqrt(diag(x$var[-arm_index, -arm_index, drop=FALSE]))
-    if (is.null(coef) | is.null(se)) 
-      stop("Input is not valid")
-    if (is.null(x$naive.var)) {
-      tmp <- cbind(coef, exp(coef), se, coef/se, 1 - stats::pchisq((coef/se)^2, 
-                                                            1))
-      dimnames(tmp) <- list(names(coef), c("coef", "exp(coef)", 
-                                           "se(coef)", "z", "p"))
-    }
-    else {
-      nse <- sqrt(diag(x$naive.var[-arm_index, -arm_index, drop=FALSE]))
-      tmp <- cbind(coef, exp(coef), nse, se, coef/se, 1 - stats::pchisq((coef/se)^2, 
+  if (!is.null(x$fail)) {
+    cat(" Coxph failed.", x$fail, "\n")
+    return()
+  }
+ 
+  savedig <- options(digits = digits)
+  on.exit(options(savedig))
+  coef <- x$coefficients
+  if (is.null(coef))  stop("Input is not valid")
+  rand_names <- colnames(model.matrix(x$formula_list$randomise, data=x$data))
+  rand_names <- rand_names[rand_names!="(Intercept)"]
+  arm_index <- which(names(coef) %in% rand_names)
+  coef <- coef[-arm_index, drop=FALSE]
+  se <- sqrt(diag(x$var[-arm_index, -arm_index, drop=FALSE]))
+  if (is.null(se))  stop("Input is not valid")
+  if (is.null(x$naive.var)) {
+    tmp <- cbind(coef, exp(coef), se, coef/se, 1 - stats::pchisq((coef/se)^2, 
                                                                  1))
-      dimnames(tmp) <- list(names(coef), c("coef", "exp(coef)", 
-                                           "se(coef)", "robust se", "z", "p"))
-    }
-    if(length(coef)>0){ 
-      stats::printCoefmat(tmp, signif.stars = FALSE, P.values = TRUE, 
-                 has.Pvalue = TRUE)}
-    #logtest <- -2 * (x$loglik[1] - x$loglik[2])
-    #if (is.null(x$df)) 
-    #  df <- sum(!is.na(coef))
-    #else df <- round(sum(x$df), 2)
-    #cat("\n")
-    #cat("Likelihood ratio test=", format(round(logtest, 2)), 
-    #    "  on ", df, " df,", " p=", format(1 - pchisq(logtest, 
-    #                                                  df)), "\n", sep = "")
-    omit <- x$na.action
-    cat("n=", x$n)
-    if (!is.null(x$nevent)) cat(", number of events=", x$nevent)
-    cat("\n")
-    if (length(omit)) 
-      cat("   (", stats::naprint(omit), ")\n", sep = "")
-    invisible(x)
+    dimnames(tmp) <- list(names(coef), c("coef", "exp(coef)", 
+                                         "se(coef)", "z", "p"))
+  }
+  else {
+    nse <- sqrt(diag(x$naive.var[-arm_index, -arm_index, drop=FALSE]))
+    tmp <- cbind(coef, exp(coef), nse, se, coef/se, 1 - stats::pchisq((coef/se)^2, 
+                                                                      1))
+    dimnames(tmp) <- list(names(coef), c("coef", "exp(coef)", 
+                                         "se(coef)", "robust se", "z", "p"))
+  }
+  if(length(coef)>0){ 
+    stats::printCoefmat(tmp, signif.stars = FALSE, P.values = TRUE, 
+                        has.Pvalue = TRUE)}
+  #logtest <- -2 * (x$loglik[1] - x$loglik[2])
+  #if (is.null(x$df)) 
+  #  df <- sum(!is.na(coef))
+  #else df <- round(sum(x$df), 2)
+  #cat("\n")
+  #cat("Likelihood ratio test=", format(round(logtest, 2)), 
+  #    "  on ", df, " df,", " p=", format(1 - pchisq(logtest, 
+  #                                                  df)), "\n", sep = "")
+  omit <- x$na.action
+  cat("n=", x$n)
+  if (!is.null(x$nevent)) cat(", number of events=", x$nevent)
+  cat("\n")
+  if (length(omit)) 
+    cat("   (", stats::naprint(omit), ")\n", sep = "")
+  invisible(x)
 }
 
 #'modified version of print.survreg
diff --git a/R/untreated.R b/R/untreated.R
@@ -25,23 +25,26 @@ untreated <- function(psi, response,treatment_matrix, rand_matrix, censor_time,
   
  
   nontreatment <- 1-rowSums(treatment_matrix) #apply(treatment_matrix,1,sum)
- # won't work with treatment_modifier at present!!!
-  #treatment_matrix <- sweep(treatment_matrix,2, exp(psi), FUN="*")
-  #treatment <- rowSums(treatment_matrix)#apply(treatment_matrix,1, sum)
-  treatment <- treatment_matrix %*% exp(psi)
+  treatment <- rowSums(treatment_matrix*exp(psi))
   
   u <- time * (nontreatment + treatment)
   
   
   #make use of setting censor_time=Inf to avoid recensoring, and implimenting Autoswitch
   
-  c_star <- censor_time* min(1, exp(psi) ) #apply(cbind( censor_time , censor_time %o% exp(psi)),1, min)
-  if( autoswitch){
+  
+  
+  c_star <- censor_time* apply(cbind(1, exp(psi)),1, min)#min(1, exp(psi) ) #apply(cbind( censor_time , censor_time %o% exp(psi)),1, min)
+  if( autoswitch & ncol(treatment_matrix)==1){
+    check <- aggregate( treatment_matrix[,1], by=list(rand_matrix[,-1]), FUN=sd)
+    switch_values <- check[check[,2]==0,1]
+    c_star <- ifelse( rand_matrix[,-1] %in% switch_values, Inf, c_star)
+    
     #if( all(rx[arm==1]==1)){ c_star <- ifelse(arm==1, Inf, c_star)}
     #if( all(rx[arm==0]==0)){ c_star <- ifelse(arm==0, Inf, c_star)}
   }
   
-  t_star <-  (u < c_star)*(u-c_star) + c_star  #pmin(u, c_star)
+  t_star <-  pmin(u, c_star) # doesn't work if c_star=Inf  : (u < c_star)*(u-c_star) + c_star 
   #only change delta if necessary
   delta_star <-  (u < c_star)*delta #ifelse( c_star < u, 0, delta)
   output <- cbind("time"=t_star, "status"=delta_star)
diff --git a/build_dashboard.R b/build_dashboard.R
@@ -0,0 +1,37 @@
+
+rcmdcheck::rcmdcheck()
+
+
+devtools::build_vignettes(clean=FALSE)
+devtools::clean_vignettes()
+
+
+
+.libPaths("U:/My Documents/R/win-library/3.6")
+devtools::build_vignettes()
+devtools::document(roclets = c('rd', 'collate', 'namespace'))
+devtools::build()
+
+#Check you've got the right version number
+install.packages("../rpsftm_1.2.6.tar.gz",
+                 repos = NULL, type = "source")
+devtools::build("../rpsftm_1.2.6.tar.gz", binary=TRUE)
+
+install.packages("V:/STATISTICS/NON STUDY FOLDER/Software/R/R Code Library/cctu_0.9.0.zip",
+                 repos = NULL, type="win.binary")
+
+#devtools::build( binary=TRUE)
+# modify news.md, cran-comments.md README.rmd
+devtools::check_win_release()
+rcmdcheck::rcmdcheck()
+devtools::check_win_devel()
+devtools::check_rhub()
+spelling::spell_check_package()
+#check github is totally up to date
+devtools::release()
+devtools::submit_cran()
+#add a tag to github
+#do some publicity , edit the website news, CTU stats email list.
+#Add the .9000 suffix to the Version field in the DESCRIPTION to
+#indicate that this is a development version.
+#Create a new heading in NEWS.md and commit the changes.
diff --git a/tests/testthat/test_rpsftm.R b/tests/testthat/test_rpsftm.R
@@ -6,8 +6,8 @@ context("Test the rpsftm() function")
 
 test_that("check of test that with mixed data sources",{
   propX <- with(immdef,I(1-xoyrs/progyrs))
-  fit <- lm(progyrs~propX, data=immdef)
-  expect_is(fit$coefficients, "numeric")
+  fit <- rpsftm(Surv(progyrs,prog)~rand(propX~imm), data=immdef,censyrs)
+  expect_is(fit$psi, "numeric")
 })
 
 test_that("first basict fit",{
