An introductory guide to analysing multiple imputed exposome data with R package rexposome
. The areas covered in this document are: loading the multiple imputations of both exposures and phenotypes from common data.frame
s, exploration the exposome data, and testing association between exposures and health outcomes.
rexposome 1.28.0
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 5
Before 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 120
Once 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: 223
class(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 57
summary(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 :3
If 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 59
The format of the multiple imputation data for rexposome
needs to follow some restrictions:
data.frame
.data.frame
must have a column called .imp
indicating the number of imputation. This imputation tagged as 0
are raw exposures (no imputation).data.frame
must have a column called .id
indicating the name of samples. This will be converted to character.data.frame
with the description with the relation between exposures and families.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)
## NULL
The 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 explanations
The output of this object indicates that we loaded 14 exposures, being 13 continuous and 1 categorical.
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] "id001" "id002" "id003" "id004" "id005" "id006"
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 Dichlorodiphenyldich.. numeric
## DDT Organochlorines DDT Dichlorodiphenyltric.. numeric
## HCB Organochlorines HCB Hexachlorobenzene numeric
pData
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 id001 never no no no male
## 2 0 id002 never no no no male
## 3 0 id003 7-12 epi no no yes male
## age cbmi blood_pre .imp.1 .id.1
## <factor> <numeric> <numeric> <numeric> <character>
## 1 4.2 16.3 120 0 id001
## 2 4.2 16.4 121 0 id002
## 3 4.2 19.0 120 0 id003
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 outside the scale range
## (`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
.
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: id001, ..., id108
## . phenoData: 109 individuals 10 phenotypes
## . individuals: id001, ..., id108
## . 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: id001, ..., id108
## . phenoData: 109 individuals 10 phenotypes
## . individuals: id001, ..., id108
## . phenotypes: whistling_chest, ..., .imp.1
## . featureData: 47 exposures 8 explanations
## . exposures: AbsPM25, ..., Temp
## . descriptions: Family, ..., .imp
## experimentData: use 'experimentData(object)'
## Annotation:
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: 7
## Robust standar errors: Computed
As 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)