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     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

Data Format

The format of the multiple imputation data for rexposome needs to follow some restrictions:

1. Both the exposures and the phenotypes are stored in the same data.frame.
2. This data.frame must have a column called .imp indicating the number of imputation. This imputation tagged as 0 are raw exposures (no imputation).
3. This data.frame must have a column called .id indicating the name of samples. This will be converted to character.
4. A data.frame with the description with the relation between exposures and families.

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) 
## 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.

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"

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

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

plotFamily(ex_imp, family = "Organochlorines")

## Warning: Removed 104 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

plotFamily(ex_imp, family = "Home Environment")

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: 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:

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.

1. First we get the threshold of effective tests

(thr <- tef(as_iew))

## [1] 0.002426352

2. Second we get the table of P-Values

tbl <- extract(as_iew)

3. Third we filter the table with the threshold

(sig <- tbl[tbl$pvalue <= thr, ])

## DataFrame with 11 rows and 4 columns
##              pvalue    effect      X2.5     X97.5
##           <numeric> <numeric> <numeric> <numeric>
## AbsPM25 5.26794e-05  20.09668  10.68607  29.50729
## NOx     7.96448e-05  12.85241   6.69541  19.00941
## NO      8.27436e-05  10.67213   5.55743  15.78682
## NO2     9.90808e-05  15.06292   7.77227  22.35356
## PM25CU  4.00711e-04  11.39992   5.23645  17.56340
## PM25    5.62435e-04  35.76444  16.11183  55.41704
## Temp    1.11891e-03  88.57154  36.30290 140.84017
## PCB180  1.37441e-03   4.00055   1.59345   6.40765
## PCB153  1.67554e-03   3.99374   1.54289   6.44459
## PM10Fe  2.07052e-03  13.57084   5.06769  22.07398
## PCB138  2.21163e-03   4.48011   1.65837   7.30184