R/aeefit.R

In github-js/acc: Exploring Accelerometer Data

#' @export
#' @importFrom utils head tail 
#' @importFrom plyr ddply
#' @importFrom nleqslv nleqslv
#' @importFrom stats model.matrix
#' @importFrom methods getClass

##############################################################################
# User's Main Function
##############################################################################
aeefit <- function(formula, data, weight=NULL, se="Sandwich", control=list(), boot=NULL) {
  
  method <- "AEE"
  Call <- match.call()
  fm <- formula
  
  #if (is.null(weight)) {
  #  weight <- rep(1,nrow(data[!duplicated(data[,colnames(data)==formula[[2]]$ID]),]))
  #} 
  
  #if (is.null(weight)) {
  #  weight <- rep(1,nrow(X))
  #} 
  
  # A PanelSurv object
  obj <- eval(formula[[2]], data)
  
  # Combine respones data frame and covariate data frame (remove intercept column)
  # Multiple rows per subject
  formula[[2]] <- NULL
  
  if (formula == ~ 1) {
    DF <- cbind(obj$psDF, zero=0)
  } else {
    DF <- cbind(obj$psDF, model.matrix(formula, data))[, -4]
  }
  
  DF <- DF[order(DF$ID, DF$time), ]
  
  # Design matrix, one row per subject
  X <- as.matrix(ddply(DF, "ID", head, n=1)[, -c(1:3)])
  
  # Create an Engine object
  engine.control <- control[names(control) %in% names(attr(getClass(method), "slots"))]
  engine <- do.call("new", c(list(Class=method), engine.control))
  
  if (length(engine@betaInit) == 1 & ncol(X) > 1)
    engine@betaInit <- rep(engine@betaInit, ncol(X))
  if (length(engine@betaInit) > 1 & length(engine@betaInit) != ncol(X))
    stop("Invalid length of initial beta values!")
  
  # Create a StdErr object
  #if(se == "NULL"){
  #  stdErr <- NULL}
  #if(se != "NULL"){
    stdErr.control <- control[names(control) %in% names(attr(getClass(se), "slots"))]
    stdErr <- do.call("new", c(list(Class=se), stdErr.control))
  #}
  
  if(se=="Sandwich"){
  fit <- doPanelFit.AEE.Sandwich(DF=DF, panelMatrix=obj$panelMatrix, timeGrid=obj$timeGrid,
                                 X=X, engine=engine, stdErr ,weight=weight)
  }
    
  if(se=="Bootstrap"){
      fit <- doPanelFit.AEE.Bootstrap(DF=DF, panelMatrix=obj$panelMatrix, timeGrid=obj$timeGrid,
                                     X=X, engine=engine, stdErr ,weight=weight, boot)
  }  
    
  ret = list(formula=fm, beta=fit$beta, 
        baseline=fit$baseline,
        timeGrid=fit$timeGrid,
        lambda=fit$lambda,
        convergence=fit$convergence,
        iter=fit$iter,
        betaSE=fit$betaSE,
        betaVar=fit$betaVar,
        baselineSE=fit$baselineSE)
  
  class(ret) <- "aeefit"
  ret
  
}