R/FLMixedmodel.R

In Fuzzy-Logix/AdapteR: Converts R syntax to SQL to provide fast analytics using in-database C++ library DB Lytix.

## Reading material.
## http://ase.tufts.edu/gsc/gradresources/guidetomixedmodelsinr/mixed%20model%20guide.html
## http://www.bodowinter.com/tutorial/bw_LME_tutorial1.pdf
##http://www.bodowinter.com/tutorial/bw_LME_tutorial2.pdf 

#' Fit Linear Mixed-Effects Models
#'
#' \code{lmer} performs mixed model  on FLTable objects.
#'
#' The DB Lytix function called is FLLinMixedModel. The mixed model extends the linear
#' model by allowing correlation and heterogeneous variances in the
#' covariance matrix of the residuals. Fit a linear mixed-effects model (LMM) to data,
#' FLLinMixedModel estimates the coefficients and covariance matrix of the mixed model
#' via the expectations.
#'
#' @seealso \code{\link[lme4]{lmer}} for R reference implementation.
#' @param formula A symbolic description of model to be fitted
#' @param data An object of class FLTable.
#' @slot results cache list of results computed
#' @slot table Input data object
#' @method residuals FLMix 
#' @method plot FLMix
#' @method summary FLMix
#' @method predict FLMix
#' @method print FLMix
#' @method coefficients FLMix
#' @method AIC FLMix
#' @method logLik FLMix
#' @return \code{lmer} returns an object of class \code{FLMix} or \code{FLMixUDT}
#' @section Constraints:
#' DBLytix currently supports one Fixed and upto 2 Random Effect variables with intercept.
#' For non-categorical Fixed Effect, only one Random Effect is supported.
#' Specify if Fixed Effect is categorical using categoricalFixedEffect logical input to the function.
#' The intercept coefficient returned is the mean value for different Random Effects.
#' @examples
#' ## One Random Effect.
#' fltbl  <- FLTable(getTestTableName("tblMixedModel"), "ObsID")
#' colnames(fltbl) <- tolower(colnames(fltbl))
#' flmod <- lmer(yval ~ fixval + (1 | ranval), data = fltbl)
#' flpred <- predict(flmod)
#' ## One Categorical Fixed Effect and Two Random Effects.
#' fltbl  <- FLTable(getTestTableName("tblMixedModelInt"), "ObsID")
#' colnames(fltbl) <- tolower(colnames(fltbl))
#' flmod <- lmer(yval ~ fixval + (1 | ranval1) + (1 | ranval2 ), fltbl, categoricalFixedEffect=TRUE)
#' flpred <- predict(flmod)
#' ## One Categorical Fixed Effect and One Random Effect.
#' flmod <- lmer(yval ~ fixval + (1 | ranval1), data = fltbl, categoricalFixedEffect=TRUE)
#' flpred <- predict(flmod)
#' @export
lmer <- function (formula,data=list(),...) {
    UseMethod("lmer", data)
}

#' @export
lmer.default <- function (formula,data=list(),...) {
    if (!requireNamespace("lme4", quietly = TRUE)){
        stop("lme4 package needed for mixed models. Please install it.",
             call. = FALSE)
    }
    else
        return(lme4::lmer(formula,data=data,...))
}

#' @export
setClass(
    "FLMix",
    slots = list(formula = "formula",
                 scoreTable = "character",
                 results = "list",
                 table = "FLTable"
                 ) )



#' @export
lmer.FLTable <- function(formula, data, fetchID = TRUE,
                        maxiter = 10,categoricalFixedEffect=FALSE,...)
{
    ## Call mixUDT if the fixed effect variable is a factor
    if(categoricalFixedEffect){
      return(mixUDT(formula=formula,data=data,
                    fetchID=fetchID,...))
    }
    vcallObject <- match.call()
    vform <- as.character(formula)
    vform <- paste0(vform[2],vform[1],vform[3])
    # vform <- as.character(vcallObject)[2]
    vreg <- gsub(pattern = "[\\|\\)]", replacement = "", x = regmatches(vform, gregexpr("\\|.*?\\)", vform))[[1]])
    Rvar <- gsub("[[:space:]]", "", vreg)
    Dvar <- gsub(pattern = " ", replacement = "", x = gsub(x = vform, pattern = "~.*", replacement = ""))
    Dvar <- gsub("[[:space:]]", "", Dvar)
    Fvar <- regmatches(vform, gregexpr("\\~.*?\\(", vform))
    Fvar <- gsub("[~.*?(]", "",Fvar)
    Fvar <- gsub("[[:space:]]", "", Fvar)
    Fvar <- strsplit(Fvar, split = "+", fixed = TRUE)
    Fvar <- unlist(Fvar)
    myformula <- as.formula(paste0(Dvar,"~ ",paste0(Fvar,collapse =  " + ")," + ",paste0(Rvar,collapse =  " + ")))
    deeptblname <- gen_unique_table_name("deepmixlin")
    vArgs <- c(Dvar, Fvar, unlist(Rvar))
    if(!isDeep(data))
        {
    FLdeep <- prepareData(formula         = myformula ,
                          data            = data,
                          outDeepTable    = deeptblname,
                          makeDataSparse  = 1,
                          performVarReduc = 0,
                          minStdDev       = .01,
                          maxCorrel       = .8,
                          fetchIDs        = FALSE)}
    vmap <- FLdeep$vmapping
    outtblname <- gen_unique_table_name("mixedtbl")
    
    data <- setAlias(data,"")
    functionName <- "FLLinMixedModel"
    vlen <- 1 + length(Fvar) + length(Rvar)
    vchr <- c("DEPENDENT", rep("FIXED", length(Fvar)), rep("RANDOM", length(Rvar)))
    vterm <- c(1,1:length(Fvar),1:length(Rvar))
    vclass <- c(0,rep(0, length(Fvar)), rep(1, length(Rvar)))

    ## use inserintotbl function
    var <- lapply(1:vlen, function(i){
        paste0("INSERT INTO fzzlLinMixedModelSpec VALUES (",fquote(outtblname)," , ",fquote(vchr[[i]]),", ",vterm[[i]],", ",vmap[[vArgs[[i]]]]," , ",vclass[[i]],");")})

    sqlSendUpdate(connection, var)
    vinputcols <- c(TableName = deeptblname ,
                    ObsIDCol = FLdeep$deepx@select@variables$obs_id_colname,
                    VarIDCol = FLdeep$deepx@select@variables$var_id_colname,
                    ValueCol = FLdeep$deepx@select@variables$cell_val_colname,
                    SpecID = outtblname,
                    MaxIter =maxiter ,
                    Note = "Mixed model for AdapteR")


    ret <- sqlStoredProc(connection,
                         functionName,
                         pInputParams = vinputcols,
                         outputParameter = c(OutTable = 'a')
                         )
    vAnalysisID <- ret[[1]]

    vin <-paste0("SELECT VarType,CoeffVal AS coeff, StdErr,TStat FROM fzzlLinMixedModelCoeffs WHERE AnalysisID = '",vAnalysisID,"' AND VarType LIKE ANY('%INTERCEPT%', '%FIXED%') ORDER BY VarType ;")

    vin <- sqlQuery(connection, vin)
    colnames(vin) <- tolower(colnames(vin))
    str <- paste0("SELECT ParamName, ParamVal FROM fzzlLinMixedModelStats WHERE AnalysisID = '",vAnalysisID,"' ORDER BY ParamName")
    vdf <- sqlQuery(connection, str)
    colnames(vdf) <- tolower(colnames(vdf))
    
    return(new("FLMix",
               formula=formula,
               table=data,
               results=list(call=vcallObject,
                            AnalysisID=vAnalysisID,
                            pArgs = list(Dvar=Dvar, Fvar = Fvar, Rvar = Rvar ), 
                            deeptbl = FLdeep$deepx,
                            vspec = outtblname,
                            vin = vin,
                            vdf = vdf,
                            scoreTable="",
                            vformula=vform
                            )))   
}




#' @export
`$.FLMix`<-function(object,property){
                                        #parentObject <- deparse(substitute(object))
    parentObject <- unlist(strsplit(unlist(strsplit(as.character(sys.call()),"(",fixed=T))[2],",",fixed=T))[1]

    if(property == "AIC"){
        df <- object@results$vdf
        return(df[df$paramname == "AIC",]$paramval) }


    if(property == "logLik"){
        df <- object@results$vdf
        return(df[df$paramname == "LogLikeliHood",]$paramval) }

    if(property == "CovarErr"){
        df <- object@results$vdf
        return(df[df$paramname == "CovarErr",]$paramval) }

    if(property == "CovarRandom"){
        df <- object@results$vdf
        return(df[df$paramname == "CovarRan",]$paramval)
        }

    if(property == "u"){
        str <- paste0("SELECT CoeffVal as coeffval, ClassVal as classval FROM fzzlLinMixedModelCoeffs WHERE AnalysisID = '",
                      object@results$AnalysisID,"' AND VarType LIKE '%RANDOM%' ORDER BY ClassVal")
        df <- sqlQuery(connection, str)
        vres <- df$coeffval
        names(vres) <- paste0("Random ",vres$classval)
        return(vres)}

    if(property == "fixedcoef"){
        str <- paste0("SELECT CoeffVal as coeffval, ClassVal as classval FROM fzzlLinMixedModelCoeffs WHERE AnalysisID = '",
                      object@results$AnalysisID,"' AND VarType LIKE '%FIXED%' ORDER BY ClassVal")
        df <- sqlQuery(connection, str)
        return(df$coeffval)
    }   
}

#' @export
setMethod("names", signature("FLMix"), function(x) {c("AIC","logLik",
                                                          "CovarErr","CovarRandom",
                                                          "u",
                                                          "fixedcoef" )})



#' @export
predict.FLMix <- function(object,
                          newdata = object@results$deeptbl,
                          scoreTable = "")
{
    parentObject <- unlist(strsplit(unlist(strsplit(
        as.character(sys.call()),"(",fixed=T))[2],")",fixed=T))[1]
    if(!any(names(object@results) == "pred"))
{
    scoretbl <- gen_score_table_name("mixed")
    vinputcols <- list(InTable = newdata@select@table_name,
                       ObsIDCol = newdata@select@variables$obs_id_colname,
                       VarIDCol = newdata@select@variables$var_id_colname,
                       ValueCol = newdata@select@variables$cell_val_colname,
                       SpecID = object@results$vspec,                       
                       pAnalysisID = as.character(object@results$AnalysisID),
                       ScoreTable = scoretbl)
    
    functionName <- "FLLinMixedModelScore"

    ret <- sqlStoredProc(connection,
                         functionName,
                         pInputParams = vinputcols,
                         outputParameter = c(OutTable = 'a')
                         )


    
    val <- new("FLSimpleVector",
               select = new("FLSelectFrom",
                            table_name = scoretbl,
                            connectionName = getFLConnectionName(),
                            variables = list(ObsID = object@table@select@variables$obs_id_colname,
                                             pred = "predval"),
                            whereconditions = "",
                            order = ""),
               dimColumns = c("obsid", "pred"),
               Dimnames = list(row = newdata@Dimnames[[1]]),
               dims = as.integer(nrow(newdata)),
               type = "integer"
               )
    object@results <- c(object@results,list(pred = val))
    assign(parentObject,object,envir=parent.frame())
    return(val)}
    else
        return(object@results$pred)
    }


#' @export
AIC.FLMix <- function(object, ...){
    return(object$AIC)}


#' @export
logLik.FLMix <- function(object, ...){
    return(object$logLik)}

#' @export
residuals.FLMix <- function(object,newdata = object@results$deeptbl, ...){
    parentObject <- unlist(strsplit(unlist(strsplit(
        as.character(sys.call()),"(",fixed=T))[2],")",fixed=T))[1]
    if(any(names(object@results) == "pred"))
        flpred <- object@results$pred
    else
        flpred <- predict(object)
    tbl <- newdata@select@table_name
    vob <- newdata@select@variables$obs_id_colname
    str <- paste0("SELECT (a.pred -b.",object@results$pArgs$Dvar,") AS res , a.",
                vob," AS obsid FROM (",constructSelect(flpred),") AS a , ",
                object@table@select@table_name," AS b WHERE a.",vob," = b.",vob," ")

    tblfunqueryobj <- new("FLTableFunctionQuery",
                          connectionName = getFLConnectionName(),
                          variables = list(
                              obs_id_colname = "obsid",
                              cell_val_colname = "res"),
                          whereconditions="",
                          order = "",
                          SQLquery=str)
    val <- new("FLSimpleVector",
               select = tblfunqueryobj,
               dimColumns = c("obsid", "res"),
               Dimnames = list(row = newdata@Dimnames[[1]]),
               dims = as.integer(nrow(newdata)),
               type = "integer"
               )
    object@results <- c(object@results,list(res = val))
    assign(parentObject,object,envir=parent.frame())

    return(val)
}


#' @export
plot.FLMix <- function(object,limit = 1000, ...){
    vres <- residuals(object)
    vpred <- predict(object)
    p <- min(limit,vpred@dims)/(vpred@dims)
    
    str1 <- paste0("SELECT  b.pred AS pred, a.res AS res FROM (",constructSelect(vpred),") AS b, (",
                constructSelect(vres),") AS a WHERE a.ObSID = b.ObsID AND FLSimUniform(RANDOM(1,10000), 0.0, 1.0) < ",p," ")
    df <- sqlQuery(connection, str1)
    return(plot(df$pred, df$res))}


#' @export
print.FLMix <- function(object, ...)
{
    cat("Linear mixed model fit by Expectation Maximization method \n\n")
    cat("Formula: ")
    print(object@results$vformula)
    cat("Data: \n")
    val <- c("AIC" = object$AIC, "logLik" = object$logLik, "CovErr" = object$CovarErr, "df.residual" = (object@results$vdf$vobs -1))
    print.default(val, digits = 3, print.gap = 2L, quote = FALSE)

    cat("Random Effects: \n")
    val <- c("Groups" = object@results$pArgs$Rvar , "Name" = "(Intercept)", "Std.Dev" = object$CovarRandom^.5)
    print(val, digits = 2, quote = FALSE)
    cat("         Residual           ", "     ", object$CovarErr^.5, "          \n\n")
    cat("Number of obs: ",object@results$vdf[7,2],", groups: ", object@results$pArgs$Rvar, ", ",object@results$vdf[8,2],"\n")
    
    cat("Fixed Effects: \n")
    cat("(Intercept) ",object@results$pArgs$Fvar,"\n", sep = "  ")
    cat(object@results$vin[object@results$vin$vartype == "Intercept", ]$coeff, 
      object@results$vin[object@results$vin$vartype == "FIXED", ]$coeff,"\n" , sep = "  ")
    ##print.default(val, digits = 3, print.gap = 2L, quote = FALSE)
    ##val <- c("FIXED" =)
    ##print.default(val, digits = 3, print.gap = 2L, quote = FALSE)
}

#' @export
setMethod("show","FLMix",function(object){print.FLMix(object)})


#' @export
summary.FLMix <- function(object, ...){
    return(print(object))
}