R/get.coef.R

In mdbrown/TreatmentSelection: Evaluate Treatment Selection Biomarkers

get.coef <-
function(formula, treatment.name, data, study.design, rho, family, ...){ 
  

  
  if(family$family == "risks_provided"){
    return(NULL)
  }else if(family$family == "time-to-event"){
    mycoxph <- coxph(formula, data = data)
    coefficients <- summary(mycoxph)$coefficients
    
  }else{
  #binary outcome
    event = data[[as.character(formula[[2]])]]
  myglm<- glm(formula, data = data, family=family)
  #myglm <- glm.fit(cbind(1, trt, marker, trt*marker), event, family=binomial(link = link))
  mycoef <- myglm$coefficients

  
  if(substr(study.design, 1, 4) == "nest") {
    mycoef[1] <- mycoef[1] - logit(mean(event)) + logit(rho[3])
    }#adjust the intercept
  else if( substr(study.design, 1,5) == "strat" ) {
    #extract the treatment
   trt = data[[treatment.name]]
   #need to adjust the interecept (mycoef[1]) and a1 (mycoef[2])
    p0 <- rho[3]/(rho[2]+rho[3])
    p1 <- rho[5]/(rho[4]+rho[5])
    #adjust intercept 
    mycoef[1] <- mycoef[1] - logit(mean(event[trt==0])) + logit(p0)
    #adjust treatment main effect
    mycoef[[treatment.name]] <- mycoef[[treatment.name]] + logit(mean(event[trt==0])) - logit(mean(event[trt==1]))  - logit(p0)+ logit(p1)
                                                                                                         
  }
 # browser()
  coefficients <- summary(myglm)$coefficients
  coefficients[,1] <- mycoef
  
  }
  return(coefficients)
}