Dealing with Multiple Imputations
29 October 2019
Abstract
An introductory guide to analysing multiple imputed exposome data with R packagerexposome. The areas covered in this document are: loading the multiple imputations of both exposures and phenotypes from common data.frames, exploration the exposome data, and testing association between exposures and health outcomes.
Package
rexposome 1.8.0
Contents
1 Introduction
1.1 Dummy Imputation with mice
To illustrate how to perform a multiple imputation using mice we start loading both rexposome and mice libraries.
library(rexposome)
library(mice)The we load the txt files includes in rexposome package so we can load the exposures and see the amount of missing data (check vignette Exposome Data Analysis for more information).
The following lines locates where the txt files were installed.
path <- file.path(path.package("rexposome"), "extdata")
description <- file.path(path, "description.csv")
phenotype <- file.path(path, "phenotypes.csv")
exposures <- file.path(path, "exposures.csv")Once the files are located we load them as data.frames:
dd <- read.csv(description, header=TRUE, stringsAsFactors=FALSE)
ee <- read.csv(exposures, header=TRUE)
pp <- read.csv(phenotype, header=TRUE)In order to speed up the imputation process that will be carried in this vignette, we will remove four families of exposures.
dd <- dd[-which(dd$Family %in% c("Phthalates", "PBDEs", "PFOAs", "Metals")), ]
ee <- ee[ , c("idnum", dd$Exposure)]We can check the amount of missing data in both exposures and phenotypes data.frames:
data.frame(
Set=c("Exposures", "Phenotypes"),
Count=c(sum(is.na(ee)), sum(is.na(pp)))
)## Set Count
## 1 Exposures 304
## 2 Phenotypes 5Before running mice, we need to collapse both the exposures and the phenotypes in a single data.frame.
rownames(ee) <- ee$idnum
rownames(pp) <- pp$idnum
dta <- cbind(ee[ , -1], pp[ , -1])
dta[1:3, c(1:3, 52:56)]## DDE DDT HCB birthdate sex age cbmi blood_pre
## id001 NA NA NA 2004-12-29 male 4.2 16.3 120
## id002 1.713577 0.6931915 1.270750 2005-01-05 male 4.2 16.4 121
## id003 2.594590 0.7448906 2.205519 2005-01-05 male 4.2 19.0 120Once this is done, the class of each column needs to be set, so mice will be able to differentiate between continuous and categorical exposures.
for(ii in c(1:13, 18:47, 55:56)) {
dta[, ii] <- as.numeric(dta[ , ii])
}
for(ii in c(14:17, 48:54)) {
dta[ , ii] <- as.factor(dta[ , ii])
}With this data.frame we perform the imputation calling mice functions (for more information about this call, check mice’s vignette). We remove the columns birthdate since it is not necessary for the imputations and carries lots of categories.
imp <- mice(dta[ , -52], pred = quickpred(dta[ , -52], mincor = 0.2,
minpuc = 0.4), seed = 38788, m = 5, maxit = 10, printFlag = FALSE)## Warning: Number of logged events: 226class(imp)## [1] "mids"The created object imp, that is an object of class mids contains 20 data-sets with the imputed exposures and the phenotypes. To work with this information we need to extract each one of these sets and create a new data-set that includes all of them. This new data.frame will be passed to rexposome (check next section to see the requirements).
mice package includes the function complete that allows to extract a single data-set from an object of class mids. We will use this function to extract the sets and join them in a single data.frame.
If we set the argument action of the complete function to 0, it will return the original data:
me <- complete(imp, action = 0)
me[ , ".imp"] <- 0
me[ , ".id"] <- rownames(me)
dim(me)## [1] 109 57summary(me[, c("H_pesticides", "Benzene")])## H_pesticides Benzene
## 0 :68 Min. :-0.47427
## 1 :35 1st Qu.:-0.19877
## NA's: 6 Median :-0.11975
## Mean :-0.12995
## 3rd Qu.:-0.06879
## Max. : 0.13086
## NA's :3If the action number is between 1 and the m value, it will return the selected set.
for(set in 1:5) {
im <- complete(imp, action = set)
im[ , ".imp"] <- set
im[ , ".id"] <- rownames(im)
me <- rbind(me, im)
}
me <- me[ , c(".imp", ".id", colnames(me)[-(97:98)])]
rownames(me) <- 1:nrow(me)
dim(me)## [1] 654 591.2 Data Format
The format of the multiple imputation data for rexposome needs to follow some restrictions:
- Both the exposures and the phenotypes are stored in the same
data.frame. - This
data.framemust have a column called.impindicating the number of imputation. This imputation tagged as0are raw exposures (no imputation). - This
data.framemust have a column called.idindicating the name of samples. This will be converted to character. - A
data.framewith the description with the relation between exposures and families.
1.3 Creating an imExposomeSet
With the exposome data.frame and the description data.frame an object of class imExposomeSet can be created. To this end, the function loadImputed is used:
ex_imp <- loadImputed(data = me, description = dd,
description.famCol = "Family",
description.expCol = "Exposure")The function loadImputed has several arguments:
args(loadImputed)## function (data, description, description.famCol = "family", description.expCol = "exposure",
## exposures.asFactor = 5, warnings = TRUE)
## NULLThe argument data is filled with the data.frame of exposures. The argument decription with the data.frame with the exposures’ description. description.famCol indicates the column on the description that corresponds to the family. description.expCol indicates the column on the description that corresponds to the exposures. Finally, exposures.asFactor indicates that the exposures with less that, by default, five different values are considered categorical exposures, otherwise continuous.
ex_imp## Object of class 'imExposomeSet'
## . exposures description:
## . categorical: 4
## . continuous: 43
## . #imputations: 6 (raw detected)
## . assayData: 47 exposures 109 individuals
## . phenoData: 109 individuals 12 phenotypes
## . featureData: 47 exposures 3 explanationsThe output of this object indicates that we loaded 14 exposures, being 13 continuous and 1 categorical.
1.3.1 Accessing to Exposome Data
The class ExposomeSet has several accessors to get the data stored in it. There are four basic methods that returns the names of the individuals (sampleNames), the name of the exposures (exposureNames), the name of the families of exposures (familyNames) and the name of the phenotypes (phenotypeNames).
head(sampleNames(ex_imp))## [1] "1" "10" "100" "101" "102" "103"head(exposureNames(ex_imp))## [1] "DDE" "DDT" "HCB" "bHCH" "PCB118" "PCB153"familyNames(ex_imp)## [1] "Organochlorines" "Bisphenol A" "Water Pollutants"
## [4] "Cotinine" "Home Environment" "Air Pollutants"
## [7] "Built Environment" "Noise" "Temperature"phenotypeNames(ex_imp)## [1] "whistling_chest" "flu" "rhinitis" "wheezing"
## [5] "sex" "age" "cbmi" "blood_pre"
## [9] ".imp.1" ".id.1"fData will return the description of the exposures (including internal information to manage them).
head(fData(ex_imp), n = 3)## DataFrame with 3 rows and 4 columns
## Family Exposure Name .type
## <character> <character> <character> <character>
## DDE Organochlorines DDE Dichlorodiphenyldichloroethylene numeric
## DDT Organochlorines DDT Dichlorodiphenyltrichloroethane numeric
## HCB Organochlorines HCB Hexachlorobenzene numericpData will return the phenotypes information.
head(pData(ex_imp), n = 3)## DataFrame with 3 rows and 12 columns
## .imp .id whistling_chest flu rhinitis wheezing sex
## <numeric> <character> <factor> <factor> <factor> <factor> <factor>
## 1 0 1 never no no no male
## 2 0 10 1-2 epi yes no no male
## 3 0 100 never no no no male
## age cbmi blood_pre .imp.1 .id.1
## <factor> <numeric> <numeric> <numeric> <character>
## 1 4.2 16.3 120 0 1
## 2 4.2 15.6 116 0 10
## 3 4.1 16.9 117 0 1001.3.2 Exposures Behaviour
The behavior of the exposures through the imputation process can be studies using the plotFamily method. This method will draw the behavior of the exposures in each imputation set in a single chart.
The method required an argument family and it will draw a mosaic with the plots from the exposures within the family. Following the same strategy than using an ExposomeSet, when the exposures are continuous box-plots are used.
plotFamily(ex_imp, family = "Organochlorines")## Warning: Removed 104 rows containing non-finite values (stat_boxplot).
For categorical exposures, the method draws accumulated bar-plot:
plotFamily(ex_imp, family = "Home Environment")
The arguments group and na.omit are not available when plotFamily is used with an imExposomeSet.
1.4 Extracting an ExposomeSet from an imExposomeSet
Once an imExposomeSet is created, an ExposomeSet can be obtained by selecting one of the internal imputed-sets. This is done using the method toES and setting the argument rid with the number of the imputed-set to use:
ex_1 <- toES(ex_imp, rid = 1)
ex_1## Object of class 'ExposomeSet' (storageMode: environment)
## . exposures description:
## . categorical: 4
## . continuous: 43
## . exposures transformation:
## . categorical: 0
## . transformed: 0
## . standardized: 0
## . imputed: 0
## . assayData: 47 exposures 109 individuals
## . element names: exp
## . exposures: AbsPM25, ..., Temp
## . individuals: 1, ..., 98
## . phenoData: 109 individuals 10 phenotypes
## . individuals: 1, ..., 98
## . phenotypes: whistling_chest, ..., .imp.1
## . featureData: 47 exposures 8 explanations
## . exposures: AbsPM25, ..., Temp
## . descriptions: Family, ..., .imp
## experimentData: use 'experimentData(object)'
## Annotation:ex_3 <- toES(ex_imp, rid = 3)
ex_3## Object of class 'ExposomeSet' (storageMode: environment)
## . exposures description:
## . categorical: 4
## . continuous: 43
## . exposures transformation:
## . categorical: 0
## . transformed: 0
## . standardized: 0
## . imputed: 0
## . assayData: 47 exposures 109 individuals
## . element names: exp
## . exposures: AbsPM25, ..., Temp
## . individuals: 1, ..., 98
## . phenoData: 109 individuals 10 phenotypes
## . individuals: 1, ..., 98
## . phenotypes: whistling_chest, ..., .imp.1
## . featureData: 47 exposures 8 explanations
## . exposures: AbsPM25, ..., Temp
## . descriptions: Family, ..., .imp
## experimentData: use 'experimentData(object)'
## Annotation:2 Exposome-Wide Association Studies (ExWAS)
The interesting point on working with multiple imputations is to test the association of the different version of the exposures with a target phenotype. rexposome implements the method exwas to be used with an imExposomeSet.
as_iew <- exwas(ex_imp, formula = blood_pre~sex+age, family = "gaussian")
as_iew## An object of class 'ExWAS'
##
## blood_pre ~ sex+age
##
## Tested exposures: 47
## Threshold for effective tests (TEF): 2.43e-03
## . Tests < TEF: NAAs usual, the object obtained from exwas method can be plotted using plotExwas:
clr <- rainbow(length(familyNames(ex_imp)))
names(clr) <- familyNames(ex_imp)
plotExwas(as_iew, color = clr)
2.1 Extract the exposures over the threshold of effective tests
The method extract allows to obtain a table of P-Values from an ExWAS object. At the same time, the tef method allows to obtain the threshold of effective tests computed at exwas. We can use them combined in order to create a table with the P-Value of the exposures that are beyond the threshold of efective tests.
- First we get the threshold of effective tests
(thr <- tef(as_iew))## [1] 0.002432295- Second we get the table of P-Values
tbl <- extract(as_iew)- Third we filter the table with the threshold
(sig <- tbl[tbl$pvalue <= thr, ])## DataFrame with 8 rows and 4 columns
## pvalue effect X2.5
## <numeric> <numeric> <numeric>
## NO 2.85681908474e-05 10.8324467183189 5.9425775324546
## NOx 2.99671099486254e-05 13.0978494675412 7.17407443122235
## AbsPM25 4.99702952492243e-05 19.8956718227279 10.6043648618683
## NO2 5.68236160709557e-05 14.8195683226825 7.84658043363589
## PM25 9.83771357407903e-05 37.1680055392154 19.0242869871666
## PM25CU 0.000296395370273705 11.4785814253346 5.4142210776783
## Temp 0.000841984485770109 88.8427612292931 37.6998443143644
## PCB153 0.00218231686364478 3.938356331383 1.45754968950307
## X97.5
## <numeric>
## NO 15.7223159041831
## NOx 19.02162450386
## AbsPM25 29.1869787835874
## NO2 21.7925562117291
## PM25 55.3117240912641
## PM25CU 17.5429417729908
## Temp 139.985678144222
## PCB153 6.41916297326293Session info
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.3 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.10-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.10-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] mice_3.6.0 lattice_0.20-38 ggplot2_3.2.1
## [4] rexposome_1.8.0 Biobase_2.46.0 BiocGenerics_0.32.0
## [7] BiocStyle_2.14.0
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-141 bitops_1.0-6 RColorBrewer_1.1-2
## [4] tools_3.6.1 backports_1.1.5 R6_2.4.0
## [7] tmvtnorm_1.4-10 rpart_4.1-15 KernSmooth_2.23-16
## [10] Hmisc_4.2-0 lazyeval_0.2.2 colorspace_1.4-1
## [13] jomo_2.6-10 nnet_7.3-12 withr_2.1.2
## [16] tidyselect_0.2.5 gridExtra_2.3 compiler_3.6.1
## [19] glmnet_2.0-18 lsr_0.5 formatR_1.7
## [22] htmlTable_1.13.2 flashClust_1.01-2 sandwich_2.5-1
## [25] labeling_0.3 bookdown_0.14 caTools_1.17.1.2
## [28] scales_1.0.0 checkmate_1.9.4 mvtnorm_1.0-11
## [31] stringr_1.4.0 digest_0.6.22 foreign_0.8-72
## [34] minqa_1.2.4 rmarkdown_1.16 base64enc_0.1-3
## [37] pkgconfig_2.0.3 htmltools_0.4.0 lme4_1.1-21
## [40] FactoMineR_1.42 htmlwidgets_1.5.1 rlang_0.4.1
## [43] GlobalOptions_0.1.1 rstudioapi_0.10 pryr_0.1.4
## [46] impute_1.60.0 generics_0.0.2 shape_1.4.4
## [49] zoo_1.8-6 gtools_3.8.1 acepack_1.4.1
## [52] dplyr_0.8.3 magrittr_1.5 Formula_1.2-3
## [55] leaps_3.0 Matrix_1.2-17 Rcpp_1.0.2
## [58] munsell_0.5.0 S4Vectors_0.24.0 imputeLCMD_2.0
## [61] lifecycle_0.1.0 scatterplot3d_0.3-41 stringi_1.4.3
## [64] yaml_2.2.0 MASS_7.3-51.4 gplots_3.0.1.1
## [67] plyr_1.8.4 grid_3.6.1 gdata_2.18.0
## [70] ggrepel_0.8.1 mitml_0.3-7 crayon_1.3.4
## [73] splines_3.6.1 circlize_0.4.8 zeallot_0.1.0
## [76] knitr_1.25 pillar_1.4.2 boot_1.3-23
## [79] reshape2_1.4.3 codetools_0.2-16 stats4_3.6.1
## [82] pan_1.6 glue_1.3.1 evaluate_0.14
## [85] latticeExtra_0.6-28 pcaMethods_1.78.0 data.table_1.12.6
## [88] BiocManager_1.30.9 vctrs_0.2.0 nloptr_1.2.1
## [91] foreach_1.4.7 tidyr_1.0.0 gtable_0.3.0
## [94] purrr_0.3.3 norm_1.0-9.5 assertthat_0.2.1
## [97] xfun_0.10 broom_0.5.2 survival_2.44-1.1
## [100] tibble_2.1.3 iterators_1.0.12 gmm_1.6-2
## [103] cluster_2.1.0 corrplot_0.84