R/run_exwas_sex.R

In groovearmada/exwas:

run_exwas_sex <- function(input_data, variable_list, sdq, subcohort = "sex", file_name) {

    exwas_result <- data.frame()

    all_imps <- mids_to_complete(input_data, variable_list)

    # extract cohorts names
    sexes <- levels(all_imps$sex)

    for (element in sexes) {

        # Filter datasets leaving one sex at a time
        sub_sex <- dplyr::filter(all_imps, sex == element) %>%
            droplevels()

        # select one exposure
        for (item in variable_list$Term) {

            # change input dataset to dataframe
            pool_fit <- sub_sex %>%

                # perform regression for each imputed dataset (group by
                # imputation number)
                dplyr::group_by(.imp) %>%

                # run simple regression for sdq ~ exposure + covariates (REMOVE SEX)
                dplyr::do(model = MASS::glm.nb(formula = get(sdq) ~ get(item) +
                                                   cohort + smoking + mother_bmi + mother_age +
                                                   mother_edu + parity + breastf_cat +
                                                   mother_work + child_age, data = .)) %>%

                # transform the regression output to a list
                as.list() %>%

                # pick the list item containing numerical regression output
                # (coefficient and other statistics for the intercept,
                # confounders and exposure), ignore other results provided by
                # the glm.nb function
                .[[-1]] %>%

                # pool the regression outocomes for each imputed dataset
                mice::pool() %>%

                # extract the final coefficient with 95% confidence intervals
                summary(conf.int = TRUE, exponentiate = TRUE)

            tidy_res <- pool_fit %>%
                mutate(sub_cohort = element) %>%
                tidy_exwas_result()

            exwas_result <- rbind(exwas_result, tidy_res)
        }
    }

    exwas_tidy_result <- tidy_exwas_outcomes(input_data, variable_list, exwas_result, subcohort = "sex", file_name)
    return(exwas_tidy_result)
}