Generating data

set.seed(2022)
{
  trials<-6
  data<-vector("list", trials)
  mylF<-vector("list", trials)
  blocklF<-vector("list", trials)
  Z<-vector("list", trials)
  bZ<-vector("list", trials)
  b.trt<-vector("list", trials)
  b.block<-vector("list", trials)


  for(i in 1:trials){
      data[[i]]<- data.frame(expand.grid(trt=1:8,block=1:4)) # create a data frame with 8 treatments and 4 blocks
      data[[i]]$site=i  # add a site id variable
      data[[i]]$y <- sample(30:70,1) + rnorm(nrow(data[[i]]),sd=8) # generate a random sample from 30 to 70 and add error Y~N(0,64)
      
      myformula <- "y ~ block + (1|trt)"  
      mylF[[i]] <- lFormula(eval(myformula), data = data[[i]]) # extract design matrix
      Z[[i]] <- mylF[[i]]$reTrms$Zt %>% as.matrix() %>% t()   # extract design matrix
      b.trt[[i]] <- rnorm(8, 2, sd=2)  # add a treatment effect. Note that effect varies between treatments but it is constant between replicates and random across trials
      
      myformula <- "y ~ (1|block)"  # add a block effect
      blocklF[[i]] <- lFormula(eval(myformula), data = data[[i]])
      bZ[[i]] <- blocklF[[i]]$reTrms$Zt %>% as.matrix() %>% t()
      b.block[[i]] <- rnorm(4, 1, sd=0.5)
      
      data[[i]]$Y <- data[[i]]$y + Z[[i]] %*% b.trt[[i]] + bZ[[i]] %*% b.block[[i]] # add everything together
      data[[i]]<-data.frame(data[[i]])
    }
}  

As we can see, the first generated data set has 32 observations, comprising of 8 treatments, 4 blocks, the column site (site identifier), y which is the outcome with no treatment or block effect added, Y which is the outcome after treatment and block effects have been structurally added to each trial. The choice of treatment, block, and overall trial effect sizes simulates, at least partially, some scenarios of yield responses from trials I am used to analyzing. Note that the treatment effect size is the same across all replicates but varies across trials. We have set things up in this way because I believe this would be the most pessimistic scenario, which treatment 1 varies across trials. Do you agree?

data[[1]]

##    trt block site        y        Y
## 1    1     1    1 68.02345 72.83317
## 2    2     1    1 65.88010 68.55153
## 3    3     1    1 67.77784 71.13164
## 4    4     1    1 51.17438 56.27865
## 5    5     1    1 70.57995 75.00631
## 6    6     1    1 57.03774 56.94811
## 7    7     1    1 65.38912 69.57787
## 8    8     1    1 55.69978 58.48249
## 9    1     2    1 65.29101 70.54203
## 10   2     2    1 73.32275 76.43549
## 11   3     2    1 66.41765 70.21274
## 12   4     2    1 56.25631 61.80188
## 13   5     2    1 77.11794 81.98561
## 14   6     2    1 87.62045 87.97212
## 15   7     2    1 71.00630 75.63636
## 16   8     2    1 67.81710 71.04112
## 17   1     3    1 63.08412 67.88555
## 18   2     3    1 53.07169 55.73483
## 19   3     3    1 66.00283 69.34833
## 20   4     3    1 70.78151 75.87748
## 21   5     3    1 63.82672 68.24479
## 22   6     3    1 65.49232 65.39439
## 23   7     3    1 54.77697 58.95744
## 24   8     3    1 67.61152 70.38594
## 25   1     4    1 73.38302 77.94689
## 26   2     4    1 60.45917 62.88475
## 27   3     4    1 62.32650 65.43444
## 28   4     4    1 66.54877 71.40719
## 29   5     4    1 63.51244 67.69295
## 30   6     4    1 56.41025 56.07476
## 31   7     4    1 65.74459 69.68749
## 32   8     4    1 72.32717 74.86403

Data wrangling and missingness generation

Next, we will be doing some data wrangling and creating a completely random missing pattern (clustered by site of course, although some may say this is not MCAR) in the response variable Y. That means that a treatment is omitted from block 1 but may still be present in blocks 2 and 3. This is important to note as this pattern will not be repeated in the second part of the exercise.

dt = do.call(rbind,data) %>% # combine trials 
  select(-y) %>%
  mutate_at(vars(block,trt),as.factor) %>% # format variables (some of the imputation methods below require dt to be numeric)
  unite("block",c(site,block),remove = FALSE) %>% # create unique id for blocks
  mutate(block = as.factor(as.numeric(as.factor(block)))) # create unique id for blocks

dt$Y=round(as.numeric(dt$Y),1)

We will be generating 10%, 20%, and 40% of missing observations per site, randomly performed across treatments. The mechanism for missingness is INDEPENDENT for each missing severity (i.e. 10%, 20%, or 40%). In other words, the missingness is not increasingly added to the previous “treated” dataset. With such design, comparisons BETWEEN missingness severities will be unfair, right? So, we focus on comparing methods within the same missingness severity.

# formula to produce NAs based on site

na.producer = function(dataset,na_percentage){
  dataset %>%
        group_by(site) %>%
        mutate(Y =  replace(Y, sample(row_number(),size = ceiling(na_percentage * n()), replace = FALSE), NA)) %>%  #Create mistiness in the dt by site
        as.data.frame() 
}  
set.seed(2022)
mis.10 = na.producer(dt,na_percentage = 0.1)
mis.20 = na.producer(dt,na_percentage = 0.2)
mis.40 = na.producer(dt,na_percentage = 0.4)

Checking the number of missing data points

In total, 4, 7 or 8, and 12 or 13 observations are lost at random per site for 10%, 20% and 40% missingness data sets, respectively.

We also can visually explore missing data patterns in this simulated data sets.

The plot above display the pattern of missingness. Note that for 10% missingness, at site 1, one observation were lost for treatment 7, 6, 4, and 3. The remaining treatments are present.

For 40% missingness at site 1, nearly all treatments except 6 have lost one or more observations. All other treatments occur at least once per site.

1) missForest

missForest is a R package used to impute missing values particularly in the case of mixed-type data. It uses a random forest trained on the observed values of a dt matrix to predict the missing values. It can be used to impute continuous and/or categorical dt including complex interactions and non-linear relations. It yields an out-of-bag (OOB) imputation error estimate without the need of a test set or elaborate cross-validation.

imp.for.10=missForest(mis.10,maxiter = 20,ntree = 500)

imp.for.20=missForest(mis.20,maxiter = 20,ntree = 500)

imp.for.40=missForest(mis.40,maxiter = 20,ntree = 500)

2) mice

The mice package in R, helps you imputing missing values with plausible data values. These plausible values are drawn from a distribution specifically designed for each missing data point. The package creates multiple imputations (replacement values) for multivariate missing data. The method is based on Fully Conditional Specification, where each incomplete variable is imputed by a separate model. The MICE algorithm can impute mixes of continuous, binary, unordered categorical and ordered categorical dt. MICE operates under the assumption that given the variables used in the imputation procedure, the missing dt are Missing At Random (MAR), which means that the probability that a value is missing depends only on observed values and not on unobserved values

predMatrix = make.predictorMatrix(data = mis.10) # matrix applies for all the missing data sets ()10%, 20%, 40%)
impMethod = make.method(data = mis.10)
impMethod["Y"] = "2l.lmer"

# specify cluster indicator and random effect
predMatrix["site","Y"] = -2 # site needs to be numeric
predMatrix["block","Y"] = 2

# ... specify levels of higher-level variables
level <- character(ncol(mis.10))
names(level) <- colnames(mis.10)

level["block"] <- "site"

# specify cluster indicators (as list)
cluster = list()

cluster[["block"]] = c("site")
cluster[["Y"]] = c("block", "site")

imp.lmer.10 = mice(mis.10, method = impMethod, predictorMatrix = predMatrix, maxit = 20, 
            m = 1, plevels_id = cluster, variables_levels = level,seed=500)

imp.lmer.20 = mice(mis.20, method = impMethod, predictorMatrix = predMatrix, maxit = 20, 
            m = 1, plevels_id = cluster, variables_levels = level,seed=500)

imp.lmer.40 = mice(mis.40, method = impMethod, predictorMatrix = predMatrix, maxit = 20, 
            m = 1, plevels_id = cluster, variables_levels = level,seed=500)

predMatrix = make.predictorMatrix(data = mis.10)
impMethod = make.method(data = mis.10)
impMethod["Y"] = "norm.predict"

# specify cluster indicator
predMatrix["site","Y"] = -2
predMatrix["block","Y"] = 2

imp.norm.pred.10 = mice(mis.10, method = impMethod,predictorMatrix = predMatrix, maxit = 1,seed=500)

imp.norm.pred.20 = mice(mis.20, method = impMethod,predictorMatrix = predMatrix, maxit = 1,seed=500)

imp.norm.pred.40 = mice(mis.40, method = impMethod,predictorMatrix = predMatrix, maxit = 1,seed=500)

3) VIM

VIM introduces new tools for the visualization of missing and/or imputed values, which can be used for exploring the data and the structure of the missing and/or imputed values. Depending on this structure of the missing values, the corresponding methods may help to identify the mechanism generating the missings and allows to explore the data including missing values. In addition, the quality of imputation can be visually explored using various univariate, bivariate, multiple and multivariate plot methods.

imp.ranger.10 = rangerImpute(Y ~ site + trt + block, mis.10)

imp.ranger.20 = rangerImpute(Y ~ site + trt + block, mis.20)

imp.ranger.40 = rangerImpute(Y ~ site + trt + block, mis.40)

imp.reg.10 = regressionImp(Y ~ site + trt + block, mis.10)

imp.reg.20 = regressionImp(Y ~ site + trt + block, mis.20)

imp.reg.40 = regressionImp(Y ~ site + trt + block, mis.40)

4) mixgb

‘mixgb’ is a scalable multiple imputation framework based on XGBoost, bootstrapping and predictive mean matching. The proposed framework is implemented in an R package mixgb. We have shown that our framework obtains less biased estimates and reflects appropriate imputation variability, while achieving high computational efficiency.

mixgb.10 = mixgb(data = mis.10, m =1,nrounds = 10)

mixgb.20 = mixgb(data = mis.20, m =1,nrounds = 10)

mixgb.40 = mixgb(data = mis.40, m =1,nrounds = 10)

Results

Below we can see the performance of imputation methods via: plot density and RMSE.

Not surprisingly, as the amount of missing data increases, the uncertainty of imputation increases and more extreme values are imputed.

We can see that methods have similar performance under low missingness severity (i.e. 10%) but larger differences arise as more and more data is missing. When accounting for visual differences solely, the imputaion method seems to matter more under higher missigness severity.

rmse = function(imp, mis_dt, truedt) {
  sqrt(sum((imp - truedt)^2)/sum(is.na(mis_dt)))
}

Missingness generation

We also will change the mechanism of data missingness. Before, observations were lost at random, but now we will add a structure to data deletion. We will omit specific treatments from half of the sites and other treatments from the remaining ones.

In other words, treatments 1 and 2 will be removed from the first three sites and similarly, treatments 7 and 8 will be removed from sites 5 and 6, simulating a scenario of MAR missing mechanism.

# formula to produce NAs based on site and treatment

mis = dt %>%
  mutate(Y = replace(Y, site %in% c(1:3) & trt %in% c(1:2),NA)) %>%
  mutate(Y = replace(Y, site %in% c(5:6) & trt %in% c(7:8),NA))

Checking the number of missing data points

Let’s take a look on the missing data.

As it was said before, treatments 1 and 2 are missing in site 1, 2, and 3. Additionally, treatments 7 and 8 are missing in sites 5 and 6.

1) missForest

Imputation is carried out through the following formula:

imp.for=missForest(mis,maxiter = 20,ntree = 500)

2) mice

Imputation is carried out through the following formula:

predMatrix = make.predictorMatrix(data = mis) 
impMethod = make.method(data = mis)
impMethod["Y"] = "2l.lmer"

# specify cluster indicator and random effect
predMatrix["site","Y"] = -2 # site needs to be numeric
predMatrix["block","Y"] = 2

# ... specify levels of higher-level variables
level <- character(ncol(mis))
names(level) <- colnames(mis)

level["block"] <- "site"
level["rain"] <- "site"

# specify cluster indicators (as list)
cluster = list()

cluster[["block"]] = c("site")
cluster[["Y"]] = c("block", "site")
cluster[["rain"]] = c("site")

imp.lmer = mice(mis, method = impMethod, predictorMatrix = predMatrix, maxit = 20, 
            m = 1, plevels_id = cluster, variables_levels = level,seed=500)

predMatrix = make.predictorMatrix(data = mis)
impMethod = make.method(data = mis)
impMethod["Y"] = "norm.predict"

# specify cluster indicator
predMatrix["site","Y"] = -2
predMatrix["block","Y"] = 2

imp.norm.pred = mice(mis, method = impMethod,predictorMatrix = predMatrix, maxit = 1,seed=500)

3) VIM

Imputation is carried out through the following formula:

imp.ranger = rangerImpute(Y ~ site + trt + block + rain, mis)

imp.reg = regressionImp(Y ~ site + trt + block + rain, mis)

4) mixgb

Imputation is carried out through the following formula:

mixgb = mixgb(data = mis, m =1,nrounds = 10)

Results