# Testing non-linear effects

### Categories and continuous variables

Developed by [Gabriel Hoffman](http://gabrielhoffman.github.io/)

#### Run on 2024-06-09 19:05:05

# 1 Introduction

Typical analysis using regression models assumes a linear affect of the covariate on the response. Here we consider testing non-linear effects in the case of 1) continuous and 2) ordered categorical variables.

We demonstrate this feature on a lightly modified analysis of PBMCs from 8 individuals stimulated with interferon-β (Kang, et al, 2018, Nature Biotech).

# 2 Standard processing

Here is the code from the main vignette:

    library(dreamlet)
    library(muscat)
    library(ExperimentHub)
    library(scater)

    # Download data, specifying EH2259 for the Kang, et al study
    eh <- ExperimentHub()
    sce <- eh[["EH2259"]]

    # only keep singlet cells with sufficient reads
    sce <- sce[rowSums(counts(sce) > 0) > 0, ]
    sce <- sce[, colData(sce)$multiplets == "singlet"]

    # compute QC metrics
    qc <- perCellQCMetrics(sce)

    # remove cells with few or many detected genes
    ol <- isOutlier(metric = qc$detected, nmads = 2, log = TRUE)
    sce <- sce[, !ol]

    # set variable indicating stimulated (stim) or control (ctrl)
    sce$StimStatus <- sce$stim

    sce$id <- paste0(sce$StimStatus, sce$ind)

    # Create pseudobulk
    pb <- aggregateToPseudoBulk(sce,
      assay = "counts",
      cluster_id = "cell",
      sample_id = "id",
      verbose = FALSE
    )

# 3 Continuous variable

Consider the continuous variable `Age`. Typical analysis only considers linear effects using a single regression coefficient, but we also want to consider the non-linear effects of age. We can peform a basis expansion using splines instead use 3 coefficients to model the age effect.

    # Simulate age between 18 and 65
    pb$Age <- runif(ncol(pb), 18, 65)

    # formula included non-linear effects of Age
    # by using a natural spline of degree 3
    # This corresponds to using 3 coefficients instead of 1
    form <- ~ splines::ns(Age, 3)

    # Normalize and apply voom/voomWithDreamWeights
    res.proc <- processAssays(pb, form, min.count = 5)

    # Differential expression analysis within each assay
    res.dl <- dreamlet(res.proc, form)

    # The spline has degree 3, so there are 3 coefficients
    # estimated for Age effects
    coefNames(res.dl)

    ## [1] "(Intercept)"          "splines::ns(Age, 3)1" "splines::ns(Age, 3)2"
    ## [4] "splines::ns(Age, 3)3"

    # Jointly test effects of the 3 spline components
    # The test of the 3 coefficients is performed with an F-statistic
    topTable(res.dl, coef = coefNames(res.dl)[2:4], number = 3)

    ## DataFrame with 3 rows and 9 columns
    ##             assay             ID splines..ns.Age..3.1 splines..ns.Age..3.2
    ##       <character>    <character>            <numeric>            <numeric>
    ## 1     CD8 T cells          COTL1             1.023185             -1.12893
    ## 2 CD14+ Monocytes C7orf55-LUC7L2             1.032476             -1.51160
    ## 3 CD14+ Monocytes          CBWD1             0.266558              1.48442
    ##   splines..ns.Age..3.3   AveExpr         F     P.Value adj.P.Val
    ##              <numeric> <numeric> <numeric>   <numeric> <numeric>
    ## 1              2.01841   8.04228   10.6116 0.000268956  0.859166
    ## 2             -2.35106   3.78659   11.6871 0.000310285  0.859166
    ## 3              2.84975   5.23680   11.4090 0.000351426  0.859166

# 4 Ordered categorical

We can also test non-linear effects in the case of categorical variables with a natural ordering to the categories. Consider time course data with 4 time points. Each time point is a category and has a natural ordering from first to last.

We have multiple options to model the time course.

- Continuous: Modeling time point as a continuous variable uses a single regression coefficient to model the linear effects of the time course. This is simple, models the order of the time points, but ignores non-linear effects

  Model using `as.numeric(TimePoint)`

- Categorical: Including time point as a typical categorical variable uses estimated the mean response value for each category. So it estimates 4 coefficients. While this can be useful for comparing two categories, it ignores the order of the time points.

  Model using `factor(TimePoint)`

- Ordered categorical: Here, the trend across ordered time points is modled using orthogonal polynomials. The trend is decomposed into independent linear, quadratic, etc., effects that can be tested either jointly or by themselves.

  Model using:

      ord <- c("time_1", "time_2", "time_3", "time_4")
      ordered(factor(TimePoint), ord)

Here we simulated 4 time points, and perform differential expression analysis.

    # Consider data generated across 4 time points
    # While there are no time points in the real data
    # we can add some for demonstration purposes
    pb$TimePoint <- ordered(paste0("time_", rep(1:4, 4)))

    # examine the ordering
    pb$TimePoint

    ##  [1] time_1 time_2 time_3 time_4 time_1 time_2 time_3 time_4 time_1 time_2
    ## [11] time_3 time_4 time_1 time_2 time_3 time_4
    ## Levels: time_1 < time_2 < time_3 < time_4

    # Use formula including time point
    form <- ~TimePoint

    # Normalize and apply voom/voomWithDreamWeights
    res.proc <- processAssays(pb, form, min.count = 5)

    # Differential expression analysis within each assay
    res.dl <- dreamlet(res.proc, form)

    # Examine the coefficient estimated
    # for TimePoint it estimates
    # linear (i.e. L)
    # quadratic (i.e. Q)
    # and cubic (i.e. C) effects
    coefNames(res.dl)

    ## [1] "(Intercept)" "TimePoint.L" "TimePoint.Q" "TimePoint.C"

    # Test only linear effect
    topTable(res.dl, coef = "TimePoint.L", number = 3)

    ## DataFrame with 3 rows and 9 columns
    ##         assay          ID     logFC   AveExpr         t     P.Value adj.P.Val
    ##   <character> <character> <numeric> <numeric> <numeric>   <numeric> <numeric>
    ## 1 CD4 T cells        DCXR -0.671645   6.52867  -5.42880 4.78666e-05  0.393058
    ## 2 CD4 T cells        GGA2 -0.890112   4.98987  -5.26372 6.69370e-05  0.393058
    ## 3 CD8 T cells        FTH1 -0.759875  14.55380  -4.99199 8.13672e-05  0.393058
    ##           B     z.std
    ##   <numeric> <numeric>
    ## 1  1.686147  -5.42880
    ## 2  0.936906  -5.26372
    ## 3  1.627753  -4.99199

    # Test linear, quadratic and cubic effcts
    coefs <- c("TimePoint.L", "TimePoint.Q", "TimePoint.C")
    topTable(res.dl, coef = coefs, number = 3)

    ## DataFrame with 3 rows and 9 columns
    ##         assay                   ID TimePoint.L TimePoint.Q TimePoint.C
    ##   <character>          <character>   <numeric>   <numeric>   <numeric>
    ## 1 CD8 T cells                 CD52    0.903032   -1.528456   -0.984065
    ## 2 CD8 T cells CCL5_ENSG00000161570    0.626201   -0.986502   -0.645468
    ## 3 CD8 T cells                  CD2    0.887578   -1.152415   -0.195616
    ##     AveExpr         F     P.Value adj.P.Val
    ##   <numeric> <numeric>   <numeric> <numeric>
    ## 1   8.69458   16.0946 1.90868e-05 0.0853298
    ## 2  11.88811   16.0433 1.94982e-05 0.0853298
    ## 3   9.67181   15.7793 2.17777e-05 0.0853298

## 4.1 Sample filtering

Due to variation in cell and read count for each sample, `processAssays()` filters out some sample. This filtering is summarized here:

    details(res.dl)

    ##               assay n_retain    formula formDropsTerms n_genes n_errors
    ## 1           B cells       16 ~TimePoint          FALSE    1961        0
    ## 2   CD14+ Monocytes       16 ~TimePoint          FALSE    3087        0
    ## 3       CD4 T cells       16 ~TimePoint          FALSE    5262        0
    ## 4       CD8 T cells       16 ~TimePoint          FALSE    1030        0
    ## 5   Dendritic cells       13 ~TimePoint          FALSE     164        0
    ## 6 FCGR3A+ Monocytes       16 ~TimePoint          FALSE    1160        0
    ## 7    Megakaryocytes       13 ~TimePoint          FALSE     172        0
    ## 8          NK cells       16 ~TimePoint          FALSE    1656        0
    ##   error_initial
    ## 1         FALSE
    ## 2         FALSE
    ## 3         FALSE
    ## 4         FALSE
    ## 5         FALSE
    ## 6         FALSE
    ## 7         FALSE
    ## 8         FALSE

Whle all 16 samples are detained in B cells, only 9 are retained for megakaryocytes. This can result in a time point being dropped, and so the polynomial expansion for some cell types can have a lower degree. The combined results will then have `NA` values for these coefficients. For example, for `TIMP1` in `Megakaryocytes` above there is not enought data to fit the cubic term, so `TimePoint.C` is `NA`.

# 5 Session Info

    ## R version 4.4.0 (2024-04-24 ucrt)
    ## Platform: x86_64-w64-mingw32/x64
    ## Running under: Windows Server 2022 x64 (build 20348)
    ##
    ## Matrix products: default
    ##
    ##
    ## locale:
    ## [1] LC_COLLATE=C
    ## [2] LC_CTYPE=English_United States.utf8
    ## [3] LC_MONETARY=English_United States.utf8
    ## [4] LC_NUMERIC=C
    ## [5] LC_TIME=English_United States.utf8
    ##
    ## time zone: America/New_York
    ## tzcode source: internal
    ##
    ## attached base packages:
    ## [1] stats4    stats     graphics  grDevices utils     datasets  methods
    ## [8] base
    ##
    ## other attached packages:
    ##  [1] GO.db_3.19.1                org.Hs.eg.db_3.19.1
    ##  [3] AnnotationDbi_1.66.0        zenith_1.6.0
    ##  [5] muscData_1.18.0             scater_1.32.0
    ##  [7] scuttle_1.14.0              ExperimentHub_2.12.0
    ##  [9] AnnotationHub_3.12.0        BiocFileCache_2.12.0
    ## [11] dbplyr_2.5.0                muscat_1.18.0
    ## [13] dreamlet_1.2.1              SingleCellExperiment_1.26.0
    ## [15] SummarizedExperiment_1.34.0 Biobase_2.64.0
    ## [17] GenomicRanges_1.56.0        GenomeInfoDb_1.40.1
    ## [19] IRanges_2.38.0              S4Vectors_0.42.0
    ## [21] BiocGenerics_0.50.0         MatrixGenerics_1.16.0
    ## [23] matrixStats_1.3.0           variancePartition_1.34.0
    ## [25] BiocParallel_1.38.0         limma_3.60.2
    ## [27] ggplot2_3.5.1               BiocStyle_2.32.0
    ##
    ## loaded via a namespace (and not attached):
    ##   [1] bitops_1.0-7              httr_1.4.7
    ##   [3] RColorBrewer_1.1-3        doParallel_1.0.17
    ##   [5] Rgraphviz_2.48.0          numDeriv_2016.8-1.1
    ##   [7] tools_4.4.0               sctransform_0.4.1
    ##   [9] backports_1.5.0           utf8_1.2.4
    ##  [11] R6_2.5.1                  metafor_4.6-0
    ##  [13] mgcv_1.9-1                GetoptLong_1.0.5
    ##  [15] withr_3.0.0               prettyunits_1.2.0
    ##  [17] gridExtra_2.3             cli_3.6.2
    ##  [19] sandwich_3.1-0            labeling_0.4.3
    ##  [21] sass_0.4.9                KEGGgraph_1.64.0
    ##  [23] SQUAREM_2021.1            mvtnorm_1.2-5
    ##  [25] blme_1.0-5                mixsqp_0.3-54
    ##  [27] parallelly_1.37.1         invgamma_1.1
    ##  [29] RSQLite_2.3.7             generics_0.1.3
    ##  [31] shape_1.4.6.1             gtools_3.9.5
    ##  [33] dplyr_1.1.4               Matrix_1.7-0
    ##  [35] metadat_1.2-0             ggbeeswarm_0.7.2
    ##  [37] fansi_1.0.6               abind_1.4-5
    ##  [39] lifecycle_1.0.4           multcomp_1.4-25
    ##  [41] yaml_2.3.8                edgeR_4.2.0
    ##  [43] mathjaxr_1.6-0            gplots_3.1.3.1
    ##  [45] SparseArray_1.4.8         grid_4.4.0
    ##  [47] blob_1.2.4                crayon_1.5.2
    ##  [49] lattice_0.22-6            beachmat_2.20.0
    ##  [51] msigdbr_7.5.1             annotate_1.82.0
    ##  [53] KEGGREST_1.44.0           magick_2.8.3
    ##  [55] pillar_1.9.0              knitr_1.47
    ##  [57] ComplexHeatmap_2.20.0     rjson_0.2.21
    ##  [59] boot_1.3-30               estimability_1.5.1
    ##  [61] corpcor_1.6.10            future.apply_1.11.2
    ##  [63] codetools_0.2-20          glue_1.7.0
    ##  [65] data.table_1.15.4         vctrs_0.6.5
    ##  [67] png_0.1-8                 Rdpack_2.6
    ##  [69] gtable_0.3.5              assertthat_0.2.1
    ##  [71] cachem_1.1.0              xfun_0.44
    ##  [73] mime_0.12                 rbibutils_2.2.16
    ##  [75] S4Arrays_1.4.1            Rfast_2.1.0
    ##  [77] coda_0.19-4.1             survival_3.7-0
    ##  [79] iterators_1.0.14          tinytex_0.51
    ##  [81] statmod_1.5.0             TH.data_1.1-2
    ##  [83] nlme_3.1-165              pbkrtest_0.5.2
    ##  [85] bit64_4.0.5               filelock_1.0.3
    ##  [87] progress_1.2.3            EnvStats_2.8.1
    ##  [89] bslib_0.7.0               TMB_1.9.11
    ##  [91] irlba_2.3.5.1             vipor_0.4.7
    ##  [93] KernSmooth_2.23-24        colorspace_2.1-0
    ##  [95] rmeta_3.0                 DBI_1.2.3
    ##  [97] DESeq2_1.44.0             tidyselect_1.2.1
    ##  [99] emmeans_1.10.2            curl_5.2.1
    ## [101] bit_4.0.5                 compiler_4.4.0
    ## [103] graph_1.82.0              BiocNeighbors_1.22.0
    ## [105] DelayedArray_0.30.1       bookdown_0.39
    ## [107] scales_1.3.0              caTools_1.18.2
    ## [109] remaCor_0.0.18            rappdirs_0.3.3
    ## [111] stringr_1.5.1             digest_0.6.35
    ## [113] minqa_1.2.7               rmarkdown_2.27
    ## [115] aod_1.3.3                 XVector_0.44.0
    ## [117] RhpcBLASctl_0.23-42       htmltools_0.5.8.1
    ## [119] pkgconfig_2.0.3           lme4_1.1-35.3
    ## [121] sparseMatrixStats_1.16.0  highr_0.11
    ## [123] mashr_0.2.79              fastmap_1.2.0
    ## [125] rlang_1.1.4               GlobalOptions_0.1.2
    ## [127] UCSC.utils_1.0.0          DelayedMatrixStats_1.26.0
    ## [129] farver_2.1.2              jquerylib_0.1.4
    ## [131] zoo_1.8-12                jsonlite_1.8.8
    ## [133] BiocSingular_1.20.0       RCurl_1.98-1.14
    ## [135] magrittr_2.0.3            GenomeInfoDbData_1.2.12
    ## [137] munsell_0.5.1             Rcpp_1.0.12
    ## [139] babelgene_22.9            viridis_0.6.5
    ## [141] EnrichmentBrowser_2.34.1  RcppZiggurat_0.1.6
    ## [143] stringi_1.8.4             zlibbioc_1.50.0
    ## [145] MASS_7.3-60.2             plyr_1.8.9
    ## [147] listenv_0.9.1             parallel_4.4.0
    ## [149] ggrepel_0.9.5             Biostrings_2.72.1
    ## [151] splines_4.4.0             hms_1.1.3
    ## [153] circlize_0.4.16           locfit_1.5-9.9
    ## [155] reshape2_1.4.4            ScaledMatrix_1.12.0
    ## [157] BiocVersion_3.19.1        XML_3.99-0.16.1
    ## [159] evaluate_0.23             RcppParallel_5.1.7
    ## [161] BiocManager_1.30.23       nloptr_2.0.3
    ## [163] foreach_1.5.2             tidyr_1.3.1
    ## [165] purrr_1.0.2               future_1.33.2
    ## [167] clue_0.3-65               scattermore_1.2
    ## [169] ashr_2.2-63               rsvd_1.0.5
    ## [171] broom_1.0.6               xtable_1.8-4
    ## [173] fANCOVA_0.6-1             viridisLite_0.4.2
    ## [175] truncnorm_1.0-9           tibble_3.2.1
    ## [177] lmerTest_3.1-3            glmmTMB_1.1.9
    ## [179] memoise_2.0.1             beeswarm_0.4.0
    ## [181] cluster_2.1.6             globals_0.16.3
    ## [183] GSEABase_1.66.0