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

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

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 = 1, 
                       description.expCol = 2)

The function loadImputed has several arguments:

args(loadImputed)

## function (data, description, description.famCol = "family", description.expCol = "exposure", 
##     exposures.asFactor = 5) 
## 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 Dichlorodiphenyldichloroethylene     numeric
## DDT Organochlorines         DDT  Dichlorodiphenyltrichloroethane     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       120         0       id003

plotFamily(ex_imp, family = "Organochlorines")

## Warning: Removed 104 rows containing non-finite values (stat_boxplot).

plotFamily(ex_imp, family = "Home Environment")

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

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.

1. First we get the threshold of effective tests

(thr <- tef(as_iew))

## [1] 0.002328967

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 9 rows and 4 columns
##                       pvalue           effect             X2.5
##                    <numeric>        <numeric>        <numeric>
## NOx     2.33754752367865e-05   13.30742440529 7.37157756537237
## NO2     3.26808155453051e-05 15.3801234451541 8.38299620908511
## AbsPM25  5.0551376725938e-05 20.0096273070629 10.6620990311518
## NO      8.65201746469424e-05 10.4578390363152 5.41017414036896
## PM25    9.07079664942412e-05 37.6132930641664   19.35876102632
## PM25CU  0.000445844607924073 11.2249431211944 5.10144233877516
## Temp     0.00101202678105872 89.2065903105215 37.0494381796216
## PCB138   0.00121485138786048 4.72037409613723 1.91411802337969
## PCB153   0.00124860796707127 4.21281048217499 1.69929217781251
##                    X97.5
##                <numeric>
## NOx     19.2432712452077
## NO2     22.3772506812232
## AbsPM25  29.357155582974
## NO      15.5055039322615
## PM25    55.8678251020127
## PM25CU  17.3484439036136
## Temp    141.363742441421
## PCB138  7.52663016889477
## PCB153  6.72632878653748