R/linear-fit.R

In ShouyeLiu/metaboliteUtility:

linear_fit_old <- function(data,
                       response,
                       variable,
                       isSeperate =T,
                       random=NULL,
                       family="gaussian")
{
    if (is.null(random)){ # non random effect involved begin
        message(paste("Use glm to model",family," data"));
        conf.high <- conf.low <- estimate <- NULL
        formulas <- as.formula(paste0(response," ~ ",paste(variable,collapse = ' + ')))
        model <- stats::glm(formulas,data = data,family=family)
        browser()
        # use broom to get all required statistics
        coef <- broom::tidy(model)
        hypothesis <- broom::glance(model)

        # for glm, calculate confidence interval
        suppressMessages(coef_conf <- confint(model))
        colnames(coef_conf) <- c("conf.low","conf.high")
        attr(coef,"row.names") <- coef$term
        coef <-  merge(coef,coef_conf,by="row.names")
        rownames(coef) <- gsub(coef$term, pattern = "`", replacement = "")
        parameter <- subset(coef, select=-c(Row.names,term))

        output <- list(parameter =parameter,hypothesis=hypothesis )
        #  non random effect involved end
    } else {

        if(family == "gaussian"){
            message(paste("Use lme4::lmer to model",family," data"));
            formulas <- as.formula(paste0(response," ~ ",paste(variable,collapse = ' + '),"+(", random,")"))
            model<- lme4::lmer(formulas,data=data);
            # use broom to get all required statistics
            suppressMessages(parameter.fixed <- broom::tidy(model, effects = "fixed", conf.int=TRUE, conf.method="profile"))
            parameter.rand <- broom::tidy(model, effects = "ran_modes", conf.int=TRUE)
            hypothesis <- broom::glance(model)
            output <- list(parameter.fixed =parameter.fixed,
                           parameter.rand = parameter.rand,
                           hypothesis=hypothesis )

        } else {
            message(paste("Use lme4::glmer to model",family," data"));
            formulas <- as.formula(paste0(response," ~ ",paste(variable,collapse = ' + '),"+(", random,")"))
            model<- lme4::glmer(formulas,data=data,family = family);

            # use broom to get all required statistics
            suppressMessages(parameter.fixed <- broom::tidy(model, effects = "fixed", conf.int=TRUE, conf.method="profile"))
            hypothesis <- broom::glance(model)
            output <- list(parameter.fixed =parameter.fixed,
                           hypothesis=hypothesis )
        }

    }

    return (output)
}

#'  Regression for data, including fixed linear model, mixed linear model and generalized mixed linear model
#'
#'
#' @param data  A data frame
#' @param response Response variables
#' @param predict A list of predict variables
#' @param random A string of random effects,Random-effects terms are distinguished by vertical bars (|) separating expressions for design matrices from grouping factors. Two vertical bars (||) can be used to specify multiple uncorrelated random effects for the same grouping variable. (Because of the way it is implemented, the ||-syntax works only for design matrices containing numeric (continuous) predictors; to fit models with independent categorical effects, see dummy or the lmer_alt function from the afex package.)
#' @param family a GLM family, see \code{\link{glm}} and \code{\link{family}}
#' @param stepwise Select a formula-based model by AIC, can be one of "both", "backward", or "forward"
#' @examples
#'    counts <- c(18,17,15,20,10,20,25,13,12)
#'    outcome <- gl(3,1,9)
#'    treatment <- gl(3,3)
#'    print(d.AD <- data.frame(treatment, outcome, counts))
#'    regression(data =d.AD,response= c("counts"), predict = c("outcome","treatment"))
#'    regression(data =d.AD,response= c("counts"), predict = c("outcome","treatment"),stepwise="backward")
#'    ## Dobson (1990) Page 93: Randomized Controlled Trial :
#'    ct <- c(18,17,15,20,10,20,25,13,12)
#'    out <- rnorm(9,1,9)
#'    treat <- rnorm(9,3,3)
#'    print(d <- data.frame(treat, out, ct))
#'    regression(data =d,response= c("ct"), predict = c("out","treat"))
#'    counts <- c(18,17,15,20,10,20,25,13,12)
#'    outcome <- gl(3,1,9)
#'    treatment <- gl(3,3)
#'    regression(data =d.AD,response= c("counts"), predict = c("outcome","treatment"),family = "poisson")
#'    a<- regression(data =d.AD,response= c("counts"),   predict= c("outcome","treatment"))
#'    library(HSAUR2)
#'    regression(response= c("outcome"), predict = c("treatment","visit"),
#'           family = "binomial",
#'           data= toenail,random = ("1|patientID"))
#'
#'
#' @name regression
#' @rdname regression
#' @export


regression <- function(data,
                       response,
                       predict,
                       random=NULL,
                       family="gaussian",
                       stepwise=NULL)
{
    response<-paste0("`",response,"`")
    predict<-paste0("`",predict,"`")
    model<-c()
    if (is.null(random)){ # non random effect involved begin
        message(paste("Use glm to model",family," data"));
        conf.high <- conf.low <- estimate <- NULL
        formulas <- as.formula(paste0(response," ~ ",paste(predict,collapse = ' + ')))
        model <- stats::glm(formulas,data = data,family=family)
        # for glm, calculate confidence interval
        suppressMessages(coef_conf <- confint.default(model))
        fit<-summary(model)
        estimate<-fit$coefficients[,c(1,2,4)]
        estimate<-cbind(estimate,coef_conf)
        #  non random effect involved end
    } else {
        if(family == "gaussian"){
            message(paste("Use lme4::lmer to model",family," data"));
            formulas <- as.formula(paste0(response," ~ ",paste(predict,collapse = ' + '),"+(", random,")"))
            model<-lme4::lmer(formulas,data=data)
            fit<-summary(model)
            # extract coefficients
            estimate <- data.frame(coef(summary(model)))
            # use normal distribution to approximate p-value
            estimate$p.z <- 2 * (1 - pnorm(abs(estimate$t.value)))

            coef_conf<-suppressMessages(confint(model,rownames(fit$coefficients),method="Wald"))
            estimate<-cbind(estimate[,c(1,2,4)],coef_conf)
        } else {
            message(paste("Use lme4::glmer to model",family," data"));
            formulas <- as.formula(paste0(response," ~ ",paste(predict,collapse = ' + '),"+(", random,")"))
            model<- lme4::glmer(formulas,data=data,family = family);
            fit<-summary(model)
            coef_conf<-suppressMessages(confint(model,rownames(fit$coefficients)))
            estimate<-fit$coefficients[,c(1,2,4)]
            estimate<-cbind(estimate,coef_conf)
        }

    }
    colnames(estimate)<-c( "Estimate","Std. Error","P.Value","Conf.low","Conf.high")
    rownames(estimate) <- gsub("`(.*)`","\\1",rownames(estimate))
    estimate<-as.data.frame(estimate)
    out=list(estimate=estimate)
    if(!is.null(stepwise))
    {
        step_result= step(model,direction = stepwise)
        out$step = step_result
    }
    return (out)
}