SOCR ≫ DSPA ≫ DSPA2 Topics ≫

This is Part 1 of the larger DSPA Visualization Chapter, which is difficult to render in a single browser window due to extreme memory demands. Visualization Chapter Part 2 includes exploratory data analytics (EDA), probability distributions, and mixture distribution modeling.

In this chapter, we will present a number of complementary strategies for data wrangling, harmonization, manipulation, aggregation, visualization, and graphical exploration. Specifically, we will discuss alternative methods for loading and saving computable data objects, importing and exporting different data structures, measuring sample statistics for quantitative variables, plotting sample histograms and model distribution functions, and scraping data from websites. In addition, we will cover exploratory data analytical (EDA) techniques, handling of incomplete (missing) data, and cohort-rebalancing of imbalanced groups.

1 Data Handling

In this section, we will discuss strategies to import data and export results. Also, we are going to learn the basic tricks we need to know about processing different types of data. Specifically, we will illustrate common R data structures and strategies for loading (ingesting) and saving (regurgitating) data. In addition, we will (1) present some basic statistics, e.g., for measuring central tendency (mean, median, mode) or dispersion (variance, quartiles, range), (2) explore simple plots, (3) demonstrate the uniform and normal distributions, (4) contrast numerical and categorical types of variables, (5) present strategies for handling incomplete (missing) data, and (6) show the need for cohort-rebalancing when comparing imbalanced groups of subjects, cases or units.

1.1 Saving and Loading R Data Structures

Let’s start by extracting Edgar Anderson’s Iris Data from the package datasets. The iris dataset quantifies morphologic shape variations of 50 Iris flowers of three related genera - Iris setosa, Iris virginica and Iris versicolor. Four shape features were measured from each sample - length and the width of the sepals and petals (in centimeters). These data were used by Ronald Fisher in his 1936 linear discriminant analysis paper.

data()
data(iris)
class(iris)
## [1] "data.frame"

As an I/O (input/output) demonstration, after we load the iris data and examine its class type, we can save it into a file named “myData.RData” and then reload it back into R.

save(iris, file="myData.RData")
load("myData.RData")

1.2 Importing and Saving Data from CSV Files

Importing the data from "CaseStudy07_WorldDrinkingWater_Data.csv" from these case-studies and saving it into the R dataset named “water”. The variables in the dataset are as follows:

  • Time: Years (1990, 1995, 2000, 2005, 2010, 2012)
  • Demographic: Country (across the world)
  • Residence Area Type: Urban, rural, or total
  • WHO Region
  • Population using improved drinking-water sources: The percentage of the population using an improved drinking water source.
  • Population using improved sanitation facilities: The percentage of the population using an improved sanitation facility.

Generally, the separator of a CSV file is comma. By default, we have optionsep=", " in the command read.csv(). Also, we can use colnames() to rename the column variables. Let’s use CaseStudy07_WorldDrinkingWater_Data.csv from out Canvas Data Archive as an example. This code loads CSV files that already include a header line containing the names of the variables. If we don’t have a header in the dataset, we can use the header = FALSE option to read the first row in the file as data. In such cases, R will assign default names to the column variables of the dataset.

water <- read.csv(
  'https://umich.instructure.com/files/399172/download?download_frd=1',
  header=TRUE, fileEncoding = "UTF-8") #, fileEncoding = "UTF-8")
water[1:3, ]
##   Year..string. WHO.region..string. Country..string.
## 1          1990              Africa          Algeria
## 2          1990              Africa           Angola
## 3          1990              Africa            Benin
##   Residence.Area.Type..string.
## 1                        Rural
## 2                        Rural
## 3                        Rural
##   Population.using.improved.drinking.water.sources......numeric.
## 1                                                             88
## 2                                                             42
## 3                                                             49
##   Population.using.improved.sanitation.facilities......numeric.
## 1                                                            77
## 2                                                             7
## 3                                                             0
colnames(water)<-c("year", "region", "country", "residence_area", "improved_water", "sanitation_facilities")
water[1:3, ]
##   year region country residence_area improved_water sanitation_facilities
## 1 1990 Africa Algeria          Rural             88                    77
## 2 1990 Africa  Angola          Rural             42                     7
## 3 1990 Africa   Benin          Rural             49                     0
which.min(water$improved_water)
## [1] 7
# rowMeans(water[,5:6])
mean(water[,6], trim=0.08, na.rm=T)
## [1] 20.4

To save a data frame to CSV files, we could use the write.csv() function. The option file = "a/local/file/path" allows us to specify the output file name and location.

write.csv(water, file = "C:/Users/water.csv")

1.3 Importing Data from ZIP and SAV Files

This example demonstrates data import from a compressed (ZIP) SPSS (SAV) file. In this case, we utilize DSPA Case-Study 25: National Ambulatory Medical Care Survey (NAMCS).

# install.packages("foreign")
library("foreign")
pathToZip <- tempfile() 
download.file("https://umich.instructure.com/files/8111611/download?download_frd=1", pathToZip, mode = "wb")
dataset <- read.spss(unzip(pathToZip, files = "namcs2015-spss.sav", list = F, overwrite = TRUE), to.data.frame=TRUE)
dim(dataset)
## [1] 28332  1096
## [1] 28332  1096
# str(dataset)
# View(dataset)  
unlink(pathToZip)

1.4 Exploring the Structure of Data

We can use the command str() and describe() to explore the structure of a dataset (in this case the CaseStudy07_WorldDrinkingWater_Data).

str(water)
## 'data.frame':    11 obs. of  6 variables:
##  $ year                 : int  1990 1990 1990 1990 1990 1990 1990 1990 1990 1990 ...
##  $ region               : chr  "Africa" "Africa" "Africa" "Africa" ...
##  $ country              : chr  "Algeria" "Angola" "Benin" "Botswana" ...
##  $ residence_area       : chr  "Rural" "Rural" "Rural" "Rural" ...
##  $ improved_water       : int  88 42 49 86 39 67 34 46 37 83 ...
##  $ sanitation_facilities: int  77 7 0 22 2 42 27 12 4 11 ...
Hmisc::describe(water)
## water 
## 
##  6  Variables      11  Observations
## --------------------------------------------------------------------------------
## year 
##        n  missing distinct     Info     Mean      Gmd 
##       11        0        1        0     1990        0 
##                
## Value      1990
## Frequency    11
## Proportion    1
## --------------------------------------------------------------------------------
## region 
##        n  missing distinct    value 
##       11        0        1   Africa 
##                  
## Value      Africa
## Frequency      11
## Proportion      1
## --------------------------------------------------------------------------------
## country 
##        n  missing distinct 
##       11        0       11 
## 
## lowest : Algeria                  Angola                   Benin                    Botswana                 Burkina Faso            
## highest: C                        Cameroon                 Central African Republic Chad                     Comoros                 
## --------------------------------------------------------------------------------
## residence_area 
##        n  missing distinct    value 
##       10        1        1    Rural 
##                 
## Value      Rural
## Frequency     10
## Proportion     1
## --------------------------------------------------------------------------------
## improved_water 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##       10        1       10        1     57.1    25.04    35.35    36.70 
##      .25      .50      .75      .90      .95 
##    39.75    47.50    79.00    86.20    87.10 
##                                                   
## Value       34  37  39  42  46  49  67  83  86  88
## Frequency    1   1   1   1   1   1   1   1   1   1
## Proportion 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
## 
## For the frequency table, variable is rounded to the nearest 0
## --------------------------------------------------------------------------------
## sanitation_facilities 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##       10        1       10        1     20.4     25.2     0.90     1.80 
##      .25      .50      .75      .90      .95 
##     4.75    11.50    25.75    45.50    61.25 
##                                                   
## Value        0   2   4   7  11  12  22  27  42  77
## Frequency    1   1   1   1   1   1   1   1   1   1
## Proportion 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
## 
## For the frequency table, variable is rounded to the nearest 0
## --------------------------------------------------------------------------------

We can see that this World Drinking Water dataset has 3331 observations and 6 variables. The output also includes the class of each variable and first few elements in the variable. The dimension of the other dataset (Case-Study 25: National Ambulatory Medical Care Survey) is much larger, \(28,332\times 1,096\).

1.5 Exploring Numeric Variables

Summary statistics for numeric variables in the dataset could be accessed by using the command summary().

library(plotly)
summary(water$sanitation_facilities)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    0.00    4.75   11.50   20.40   25.75   77.00       1
summary(water[c("improved_water", "sanitation_facilities")])
##  improved_water  sanitation_facilities
##  Min.   :34.00   Min.   : 0.00        
##  1st Qu.:39.75   1st Qu.: 4.75        
##  Median :47.50   Median :11.50        
##  Mean   :57.10   Mean   :20.40        
##  3rd Qu.:79.00   3rd Qu.:25.75        
##  Max.   :88.00   Max.   :77.00        
##  NA's   :1       NA's   :1
# plot(density(water$improved_water,na.rm = T))  # no need to be continuous, we can still get intuition about the variable distribution

fit <- density(as.numeric(water$improved_water),na.rm = T)
plot_ly(x = fit$x, y = fit$y, type = "scatter", mode = "lines", 
              fill = "tozeroy", name = "Density") %>% 
    layout(title='Density of (%) Improved Water Quality', 
           xaxis = list (title = 'Percent'), yaxis = list (title = 'Density'))

The six summary statistics and NA’s (missing data) are reported in the output.

1.6 Measuring Central Tendency - mean, median, mode

Mean and median are two frequent measurements of the central tendency. Mean is “the sum of all values divided by the number of values”. Median is the number in the middle of an ordered list of values. In R, mean() and median() functions can provide us with these two measurements.

vec1<-c(40, 56, 99)
mean(vec1)
## [1] 65
mean(c(40, 56, 99))
## [1] 65
median(vec1)
## [1] 56
median(c(40, 56, 99))
## [1] 56
# install.packages("psych"); 
library("psych")
geometric.mean(vec1, na.rm=TRUE)
## [1] 60.52866

The mode is the value that occurs most often in the dataset. It is often used in categorical data, where mean and median are inappropriate measurements.

We can have one or more modes. In the water dataset, we have “Europe” and “Urban” as the modes for region and residence area respectively. These two variables are unimodal, which has a single mode. For the year variable, we have two modes 2000 and 2005. Both of the categories have 570 counts. The year represent an example of a multimodal variable that has two, or more, modes.

Mode is one of the measures for the central tendency. The best way to use it is to compare the mode to other values in the data. This helps us determine whether one or several categories dominate all others in the data. In numeric datasets, we could think mode as the highest bin in the histogram, since it is unlikely to have many repeated measurements for continuous variables. In this way, we can also examine if the numeric data is multimodal.

More information about measures of centrality is available here.

1.7 Measuring Spread - variance, quartiles and the five-number summary

The five-number summary describes the spread of a dataset. They are:

  • Minimum (Min.), representing the smallest value in the data
  • First quartile/Q1 (1st Qu.), representing the \(25^{th}\) percentile, which splits off the lowest 25% of data from the highest 75%
  • Median/Q2 (Median), representing the \(50^{th}\) percentile, which splits off the lowest 50% of data from the top 50%
  • Third quartile/Q3 (3rd Qu.), representing the \(75^{th}\) percentile, which splits off the lowest 75% of data from the top 25%
  • Maximum (Max.), representing the largest value in the data.

Min and Max can be obtained by using min() and max() respectively.

The difference between maximum and minimum is known as range. In R, range() function gives us both the minimum and maximum. A combination of range() and diff() could do the trick of getting the actual range value. To avoid problems with missing values, we will ignore them using the option na.rm=TRUE.

range(water$improved_water, na.rm=TRUE)
## [1] 34 88
diff(range(water$improved_water, na.rm=TRUE))
## [1] 54

Q1 and Q3 are the 25th and 75th percentiles of the data. Median (Q2) is right in the middle of Q1 and Q3. The difference between Q3 and Q1 is called the interquartile range (IQR). Within the IQR lies half of our data that has no extreme values.

In R, we use the IQR() to calculate the interquartile range. If we use IQR() for a data with NA’s, the NA’s are ignored by the function while using the option na.rm=TRUE.

IQR(water$improved_water, na.rm=TRUE)
## [1] 39.25
summary(water$improved_water)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   34.00   39.75   47.50   57.10   79.00   88.00       1

Just like the command summary() that we have talked about earlier in this chapter. A similar function quantile() could be used to obtain the five-number summary.

quantile(water$improved_water, na.rm = TRUE)
##    0%   25%   50%   75%  100% 
## 34.00 39.75 47.50 79.00 88.00

We can also calculate specific percentiles in the data. For example, if we want the 20th and 60th percentiles, we can do the following.

quantile(water$improved_water, probs = c(0.2, 0.6), na.rm = TRUE)
##  20%  60% 
## 38.6 56.2

When we include the seq() function, generating percentiles of evenly-spaced values is available.

quantile(water$improved_water, seq(from=0, to=1, by=0.2), na.rm = TRUE)
##   0%  20%  40%  60%  80% 100% 
## 34.0 38.6 44.4 56.2 83.6 88.0

Let’s re-examine the five-number summary for the improved_water variable. When we ignore the NA’s, the difference between minimum and Q1 is 74 while the difference between Q3 and maximum is only 1. The interquartile range is 22%. Combining these facts, the first quarter is more widely spread than the middle 50 percent of values. The last quarter is the most condensed one that has only two percentages 99% and 100%. Also, we can notice that the mean is smaller than the median. The mean is more sensitive to the extreme values than the median. Having some of very small values may spread out the first quartile, skew the distribution to the left and make the mean less than the median.

Distribution models offer a way to characterize data using only a few parameters. For example, the normal distribution can be defined by only two parameters - center and spread, statistically speaking, mean and standard deviation.

The mean value is obtained by arithmetic averaging of all data points.

\[mean(X)=\mu=\frac{1}{n}\sum_{i=1}^{n} x_i\]

The geometric mean of a set of n non-negative numbers \(\{a_1,a_2,\cdots,a_n\}\) is given by \[G = \left( \prod_{i=1}^{n} a_i \right)^{1/n}.\]

The standard deviation is the square root of the variance. And the variance is the average sum of square deviation from the mean.

\[Var(X)=\sigma^2=\frac{1}{n-1}\sum^{n}_{i=1} (x_i-\mu)^2\] \[StdDev(X)=\sigma=\sqrt{Var(X)}\]

Since the water dataset is not close to normal, in this example, we will use MLB baseball players dataset to illustrate normal distribution properties. The MLB dataset (01a_data.txt) in our class file data archive has following variables - Name, Team, Position, Height, Weight, and Age.

We can use histograms to visually assess approximate normality of baseball players’ Height and Weight.

baseball<-read.table("https://umich.instructure.com/files/330381/download?download_frd=1", header=T)
# hist(baseball$Weight, main = "Histogram for Baseball Player's Weight", xlab="weight")
# hist(baseball$Height, main = "Histogram for Baseball Player's Height", xlab="height")

x <- rnorm(10000, mean=mean(baseball$Weight, na.rm=T), sd=sd(baseball$Weight, na.rm=T)) 
fit <- density(x, bw=10)    

plot_ly(x=~baseball$Weight, type = "histogram", name = "Weight Histogram", histnorm = "probability") %>%    
    add_trace(x =~fit$x, y =~5*fit$y, type = "scatter", mode = "lines", opacity=0.1,    
              fill = "tozeroy", name = "Normal Density") %>%    
    layout(title='Baseball Weight Histogram & Model Normal Distribution',   
           xaxis = list(title = "Weight"), yaxis = list(title = "relative frequency/density"),  
           legend = list(orientation = 'h'))    
x <- rnorm(10000, mean=mean(baseball$Height, na.rm=T), sd=sd(baseball$Height, na.rm=T)) 
fit <- density(x, bw=1) 

plot_ly(x=~baseball$Height, type = "histogram", name = "Height Histogram", histnorm = "probability") %>%    
    add_trace(x=~fit$x, y=~fit$y, type = "scatter", mode = "lines", opacity=0.1,    
              fill = "tozeroy", name = "Normal Density") %>%    
    layout(title='Baseball Height Histogram & Model Normal Distribution',   
           xaxis = list(title = "Height"), yaxis = list(title = "relative frequency/density"),  
           legend = list(orientation = 'h'))    

TWe could also report the mean and standard deviation of the weight and height variables.

mean(baseball$Weight)
## [1] 201.7166
mean(baseball$Height)
## [1] 73.69729
var(baseball$Weight)
## [1] 440.9913
sd(baseball$Weight)
## [1] 20.99979
var(baseball$Height)
## [1] 5.316798
sd(baseball$Height)
## [1] 2.305818

Larger standard deviation, or variance, suggests the data is more spread out from the mean. Therefore, for MLB players, weights appear to be more spread than heights.

Given the first two moments (mean and standard deviation), we can easily estimate how extreme a specific value is. Assuming we have a normal distribution, the values follow a \(68-95-99.7\) rule. This means 68% of the data lies within the interval \([\mu-\sigma, \mu+\sigma]\);95% of the data lies within the interval \([\mu-2*\sigma, \mu+2*\sigma]\) and 99.7% of the data lies within the interval \([\mu-3*\sigma, \mu+3*\sigma]\). The following graph plotted by R illustrates the \(68-95-99.7\) rule.

#  hist(x, probability=T, 
#    col='lightblue', xlab=' ', ylab=' ', axes = F, 
#    main='68-95-99.7 Rule')
# lines(density(x, bw=0.4), col='red', lwd=3)
# axis(1, at=c(-3, -2, -1, 0, 1, 2, 3), labels = expression(mu-3*sigma, mu-2*sigma, mu-sigma, mu, mu+sigma, mu+2*sigma, mu+3*sigma))
# abline(v=-1, lwd=3, lty=2)
# abline(v=1, lwd=3, lty=2)
# abline(v=-2, lwd=3, lty=2)
# abline(v=2, lwd=3, lty=2)
# abline(v=-3, lwd=3, lty=2)
# abline(v=3, lwd=3, lty=2)
# text(0, 0.2, "68%")
# segments(-1, 0.2, -0.3, 0.2, col = 'red', lwd=2)
# segments(1, 0.2, 0.3, 0.2, col = 'red', lwd=2)
# text(0, 0.15, "95%")
# segments(-2, 0.15, -0.3, 0.15, col = 'red', lwd=2)
# segments(2, 0.15, 0.3, 0.15, col = 'red', lwd=2)
# text(0, 0.1, "99.7%")
# segments(-3, 0.1, -0.3, 0.1, col = 'red', lwd=2)
# segments(3, 0.1, 0.3, 0.1, col = 'red', lwd=2)

N<- 1000
norm <- rnorm(N, 0, 1)
#  hist(x, probability=T, 
#    col='lightblue', xlab=' ', ylab=' ', axes=F, 
#    main='Normal Distribution')
# lines(density(x, bw=0.4), col='red', lwd=3)
normDensity <- density(norm, bw=0.5)
dens <- data.frame(x = normDensity$x, y = normDensity$y)
miny <- 0
maxy <- max(dens$y)

xLabels <- c("&mu;-3&#963;","&mu;-2&#963;", "&mu;-&#963;", "&mu;", "&mu;+&#963;", "&mu;+2&#963;", "&mu;+3&#963;")
labelColors <- c("green", "red", "orange", "black", "orange", "red", "green")
xLocation <- c(-3, -2, -1, 0, 1, 2, 3)
yLocation <- 0.2
data <- data.frame(xLabels, xLocation, yLocation)

plot_ly(dens) %>% 
    add_histogram(x = norm, name="Normal Histogram") %>% 
    add_lines(data = dens, x = ~x, y = ~y+0.05, yaxis = "y2", 
              line = list(width = 3), name="N(0,1)") %>% 
    add_annotations(x = ~xLocation, y = ~yLocation, type = 'scatter', ax = 20, ay = 20,
                    mode = 'text', text = ~xLabels, textposition = 'middle right',
                    textfont = list(color = labelColors, size = 16)) %>%
    add_segments(x=-3, xend=-3, y=0, yend=100, name="99.7%", line=list(dash="dash", color="green")) %>%
    add_segments(x=-2, xend=-2, y=0, yend=90, name="95%", line=list(dash="dash", color="red")) %>%
    add_segments(x=-1, xend=-1, y=0, yend=80, name="68%", line=list(dash="dash", color="orange")) %>%
    add_segments(x=1, xend=1, y=0, yend=80, name="68%", line = list(dash = "dash", color="orange")) %>%
    add_segments(x=2, xend=2, y=0, yend=90, name="95%", line=list(dash="dash", color="red")) %>%
    add_segments(x=3, xend=3, y=0, yend=100, name="99.7%", line=list(dash="dash", color="green")) %>%
    add_segments(x=-3, xend=3, y=100, yend=100, name="99.7%", line=list(dash="dash", color="green")) %>%
    add_segments(x=-2, xend=2, y=90, yend=90, name="95%", line=list(dash="dash", color="red")) %>%
    add_segments(x=-1, xend=1, y=80, yend=80, name="68%", line=list(dash="dash", color="orange")) %>%
    layout(bargap=0.1, xaxis=list(name=""), yaxis=list(title="density/frequency"),
           yaxis2 = list(overlaying = "y", side = "right", # title="prob",
                         range = c(miny, maxy+0.1), showgrid = F, zeroline = F),
           legend = list(orientation = 'h'), title="Normal 68-95-99.7% Rule")

Applying the 68-95-99.7 rule to our baseball weight variable, we know that 68% of our players weighted between 180.7168 pounds and 222.7164 pounds; 95% of the players weighted between 159.7170 pounds and 243.7162 pounds; And 99.7% of the players weighted between 138.7172 pounds and 264.7160 pounds.

1.8 Visualizing Numeric Variables - boxplots

We can visualize the five-number summary by a boxplot (box-and-whiskers plot). With the boxplot() function we can manage the title (main="") and labels for x (xlab="") and y (ylab="") axis.

# boxplot(water$improved_water, main="Boxplot for Percent improved_water", ylab="Percentage")
plot_ly(y = ~water$improved_water, type = "box", name="improved water qual") %>% 
  add_trace(y = ~water$sanitation_facilities, name ="sanitation") %>% 
    layout(title='Boxplots of Improved Water Quality and Sanitation Facilities', 
           yaxis = list (title = 'Percent'))

In the boxplot we have five horizontal lines each representing the corresponding value in the five-number summary. The box in the middle represents the middle 50 percent of values. The bold line in the box is the median. Mean value is not illustrated on the graph.

Boxplots only allow the two ends to extend to a minimum or maximum of 1.5 times the IQR. Therefore, any value that falls outside of the \(3\times IQR\) range will be represented as circles or dots. They are considered outliers. We can see that there are a lot of outliers with small values on the low ends of the graph.

1.9 Visualizing Numeric Variables - histograms

Histograms offer another way to show the distribution spread of numeric variables. They require a specification of a number of bins, value containers, to divide and stratify the original data. The heights of the bins indicate the observed frequencies within each bin.

# hist(water$improved_water, main = "Histogram of  Percent improved_water", xlab="Percentage")
# hist(water$sanitation_facilities, main = "Histogram of  Percent sanitation_facilities", xlab = "Percentage")

plot_ly(x = ~water$improved_water, type = "histogram", name="improved_water") %>% 
  add_trace(x = ~water$sanitation_facilities, type = "histogram", name="sanitation_facilities") %>% 
  layout(bargap=0.1, title='Histograms', legend = list(orientation = 'h'),
         xaxis = list(title = 'Percent'), yaxis = list (title = 'Frequency'))

We could see that the shape of two graphs are somewhat similar. They both appear to have left skewed patterns (\(mean \lt median\)). Other common skew patterns are shown in the following graph.

N <- 10000
x <- rnbinom(N, 5, 0.1)
#  hist(x, 
#  xlim=c(min(x), max(x)), probability=T, nclass=max(x)-min(x)+1, 
#    col='lightblue', xlab=' ', ylab=' ', axes=F, 
#    main='Right Skewed')
# lines(density(x, bw=1), col='red', lwd=3)
fit <- density(x)

plot_ly(x = x, type = "histogram", name = "Data Histogram") %>% 
    add_trace(x = fit$x, y = fit$y, type = "scatter", mode = "lines", opacity=0.3,
              fill = "tozeroy", yaxis = "y2", name = "Density (rnbinom(N, 5, 0.1))") %>% 
    layout(title='Right Skewed Process', yaxis2 = list(overlaying = "y", side = "right"),
           legend = list(orientation = 'h'))
N <- 10000
x <- rnorm(N, 15, 3.7)
#  hist(x, 
#  xlim=c(min(x), max(x)), probability=T, nclass=max(x)-min(x)+1, 
#    col='lightblue', xlab=' ', ylab=' ', axes=F, 
#    main='Right Skewed')
# lines(density(x, bw=1), col='red', lwd=3)
fit <- density(x)

plot_ly(x = x, type = "histogram", name = "Data Histogram") %>% 
    add_trace(x = fit$x, y = fit$y, type = "scatter", mode = "lines", opacity=0.3,
              fill = "tozeroy", yaxis = "y2", name = "Density (rnorm(N, 15, 3.7))") %>% 
    layout(title='Symmetric Process', yaxis2 = list(overlaying = "y", side = "right"),
           legend = list(orientation = 'h'))
# N <- 10000
# xNu <- extraDistr::rlaplace(N, mu = 0, sigma = 0.4)
# yNu <- density(xNu, bw=0.2)
# xMu <- extraDistr::rlaplace(N, mu = 0, sigma = 0.5)
# yMu <- density(xMu, bw=0.2)
# # correct second Laplace Density (mu) to ensure absolute continuity, nu << mu
# yMu$y <- 2*yMu$y
# plot_ly(x = x, type = "histogram", name = "Data Histogram") %>% 
#     add_trace(x = yNu$x, y = yNu$y, type = "scatter", mode = "lines", opacity=0.3,
#               fill = "tozeroy", yaxis = "y2", name = "nu, Laplace(N,0,0.4) Density") %>% 
#     add_trace(x = yMu$x, y = yMu$y, type="scatter", mode="lines", opacity=0.3,
#               fill = "tozeroy", yaxis = "y2", name = "mu, Laplace(N,0,0.5) Density") %>% 
#     layout(title="Absolutely Continuous Laplace Distributions, nu<<mu", 
#            yaxis2 = list(overlaying = "y", side = "right"),
#            xaxis = list(range = list(-pi, pi)),
#            legend = list(orientation = 'h'))
# integrate(approxfun(yNu), -pi, pi)
# integrate(approxfun(yMu), -pi, pi)

You can learn more about Probability Distributions in the SOCR EBook and see the density plots of over 80 different probability distributions using the SOCR Java Distribution Calculators or the Distributome HTML5 Distribution Calculators.

For each probability distribution defined in R, there are four functions that provide the density (e.g., dnorm), the cumulative probability (e.g., pnorm), the inverse cumulative distribution (quantile) function (e.g., qnorm), and the random sampling (simulation) function (e.g., rnorm). The plots below show the standard normal density, cumulative probability and the quantile functions. As the density is very small outside of the interval \((-4,4)\), the plots are restricted to this domain.

z<-seq(-4, 4, 0.1)  # points from -4 to 4 in 0.1 steps
q<-seq(0.001, 0.999, 0.001)  # probability quantile values from 0.1% to 99.9% in 0.1% steps

dStandardNormal <- data.frame(Z=z, Density=dnorm(z, mean=0, sd=1), Distribution=pnorm(z, mean=0, sd=1))  

qStandardNormal <- data.frame(Q=q, Quantile=qnorm(q, mean=0, sd=1))  
head(dStandardNormal)
##      Z      Density Distribution
## 1 -4.0 0.0001338302 3.167124e-05
## 2 -3.9 0.0001986555 4.809634e-05
## 3 -3.8 0.0002919469 7.234804e-05
## 4 -3.7 0.0004247803 1.077997e-04
## 5 -3.6 0.0006119019 1.591086e-04
## 6 -3.5 0.0008726827 2.326291e-04
# plot(z, dStandardNormal$Density, main="Normal Density Curve", type = "l", xlab = "critical values", ylab="density", lwd=4, col="blue")
# polygon(z, dStandardNormal$Density, col="red", border="blue")
# plot(z, dStandardNormal$Distribution, main="Normal Distribution", type = "l", xlab = "critical values", ylab="Cumulative Distribution", lwd=4, col="blue")
# plot(q, qStandardNormal$Quantile, main="Normal Quantile Function (Inverse CDF)", type = "l", xlab = "p-values", ylab="Critical Values", lwd=4, col="blue")

plot_ly(x = z, y= dStandardNormal$Density, name = "Normal Density Curve",
        mode = 'lines') %>% 
    layout(title='Normal Density Curve', 
           xaxis = list(title = 'critical values'),
           yaxis = list(title ="Density"),
           legend = list(orientation = 'h'))
plot_ly(x = z, y= dStandardNormal$Distribution,
        name = "Normal Density Curve", mode = 'lines') %>% 
    layout(title='Normal Distribution', 
           xaxis = list(title = 'critical values'),
           yaxis = list(title ="Cumulative Distribution"),
           legend = list(orientation = 'h'))
plot_ly(x = q, y= qStandardNormal$Quantile,
        name = "Normal Quantile Function (Inverse CDF)", mode = 'lines') %>%
    layout(title='Normal Distribution', 
           xaxis = list(title = 'probability values'),
           yaxis = list(title ="Critical Values"),
           legend = list(orientation = 'h'))

1.10 Understanding Numeric Data - uniform and normal distributions

If the data follows a uniform distribution, then all values are equally likely to occur. The histogram for a uniformly distributed data would have equal heights for each bin like the following graph.

Often, but not always, real world processes may appear as normally distributed data. A normal distribution would have a higher frequency for middle values and lower frequency for more extreme values. It has a symmetric and bell-curved shape just like the following diagram generated by R. Many parametric-based statistical approaches assume normality of the data. In cases where this parametric assumption is violated, variable transformations or distribution-free tests may be more appropriate.

N<- 1000
norm <- rnorm(N, 0, 1)
#  hist(x, probability=T, 
#    col='lightblue', xlab=' ', ylab=' ', axes=F, 
#    main='Normal Distribution')
# lines(density(x, bw=0.4), col='red', lwd=3)
normDensity <- density(norm, bw=0.5)
dens <- data.frame(x = normDensity$x, y = normDensity$y)
miny <- 0
maxy <- max(dens$y)

plot_ly(dens) %>% 
  add_histogram(x = norm, name="Normal Histogram") %>% 
  add_lines(data = dens, x = ~x, y = ~y, yaxis = "y2", 
            line = list(width = 3), name="N(0,1)") %>% 
  layout(bargap=0.1, yaxis2 = list(overlaying = "y",  side = "right", 
            range = c(miny, maxy), showgrid = F, zeroline = F),
         legend = list(orientation = 'h'), title="Normal(0,1)")

1.11 Exploring Categorical Variables

Back to our water dataset, we can treat the year variable as categorical rather than a numeric variable. Since the year variable only has six distinctive values, it is rational to treat it as a categorical variable where each value is a category that could apply to multiple WHO regions. Moreover, region and residence area variables are also categorical.

Different from numeric variables, the categorical variables are better examined by tables rather than summary statistics. One-way table represents a single categorical variable. It gives us the counts of different categories. table() function can create one-way tables for our water dataset:

water <- read.csv('https://umich.instructure.com/files/399172/download?download_frd=1', header=T, stringsAsFactors=FALSE, fileEncoding="latin1")

colnames(water)<-c("year", "region", "country", "residence_area", "improved_water", "sanitation_facilities")

table(water$year)
## 
## 1990 1995 2000 2005 2010 2012 
##  520  561  570  570  556  554
table(water$region)
## 
##                Africa              Americas Eastern Mediterranean 
##                   797                   613                   373 
##                Europe       South-East Asia       Western Pacific 
##                   910                   191                   447
table(water$residence_area)
## 
## Rural Total Urban 
##  1095  1109  1127

Given that we have a total of 3331 observations, the WHO region table tells us that about 27% (910/3331) of the areas examined in the study are in Europe.

R can directly give us table proportions when using the prop.table() function. The proportion values can be transformed into percentage form and edit number of digits.

year_table<-table(water$year)
prop.table(year_table)
## 
##      1990      1995      2000      2005      2010      2012 
## 0.1561093 0.1684179 0.1711198 0.1711198 0.1669168 0.1663164
year_pct<-prop.table(year_table)*100
round(year_pct, digits=1)
## 
## 1990 1995 2000 2005 2010 2012 
## 15.6 16.8 17.1 17.1 16.7 16.6

1.12 Exploring Relationships Between Variables

So far, the methods and statistics that we have gone through are at univariate level. Sometimes we want to examine the relationship between two or multiple variables. For example, does the percentage of population that uses improved drinking-water sources increase over time? To address these problems we need to look at bivariate or multivariate relationships.

1.12.1 Visualizing Relationships - scatterplots

Let’s look at the bivariate case first. A scatterplot is a good way to visualize bivariate relationships. We have x axis and y axis each representing one of the variables. Each observation is illustrated on the graph by a glyph, e.g., a solid point. If the graph shows a clear pattern, rather than a random scatter of points or a horizontal line, the two variables may be correlated with each other.

In R we can use the plot() function to create scatterplots. We have to define the variables for x-axis and y-axis. The labels in the graph are editable.

# plot.window(c(400,1000), c(500,1000))
# plot(x=water$year, y=water$improved_water, 
#      main= "Scatterplot of Year vs. Improved_water", 
#      xlab= "Year", 
#      ylab= "Percent of Population Using Improved Water")

plot_ly(x = ~water$sanitation_facilities, y = ~water$improved_water, type = "scatter", 
        mode = "markers") %>% 
    layout(title='Scatterplot: Improved Water Quality vs. Sanitation Facilities', 
           xaxis = list (title = 'Water Quality'), yaxis = list (title = 'Sanitation'))

We can see from the scatterplot that there is an increasing pattern. In later years, the percentages are more centered around one hundred. Especially, in 2012, none of the regions had less than 20% of people using improved water sources while there used to be some regions that had such low percentages in the early years.

1.12.2 Examining Relationships - two-way cross-tabulations

Scatterplot is a useful tool to examine the relationship between two variables where at least one of them is numeric. When both variables are nominal, two-way cross-tabulation would be a better choice (also named as crosstab or contingency table).

The function CrossTable() is available in R under the package gmodels. Let’s install it first.

#install.packages("gmodels", repos = "http://cran.us.r-project.org")
library(gmodels)

We are interested in investigating the relationship between WHO region and residence area type in the water study. We might want to know if there is a difference in terms of residence area type between the African WHO region and all other WHO regions.

To address this problem we need to create an indicator variable for the African WHO region first.

water$africa<-water$region=="Africa"

Let’s revisit the table() function to see how many WHO regions are in Africa.

table(water$africa)
## 
## FALSE  TRUE 
##  2534   797

Now, let’s create a two-way cross-tabulation using CrossTable().

CrossTable(x=water$residence_area, y=water$africa)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## | Chi-square contribution |
## |           N / Row Total |
## |           N / Col Total |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  3331 
## 
##  
##                      | water$africa 
## water$residence_area |     FALSE |      TRUE | Row Total | 
## ---------------------|-----------|-----------|-----------|
##                Rural |       828 |       267 |      1095 | 
##                      |     0.030 |     0.096 |           | 
##                      |     0.756 |     0.244 |     0.329 | 
##                      |     0.327 |     0.335 |           | 
##                      |     0.249 |     0.080 |           | 
## ---------------------|-----------|-----------|-----------|
##                Total |       845 |       264 |      1109 | 
##                      |     0.002 |     0.007 |           | 
##                      |     0.762 |     0.238 |     0.333 | 
##                      |     0.333 |     0.331 |           | 
##                      |     0.254 |     0.079 |           | 
## ---------------------|-----------|-----------|-----------|
##                Urban |       861 |       266 |      1127 | 
##                      |     0.016 |     0.050 |           | 
##                      |     0.764 |     0.236 |     0.338 | 
##                      |     0.340 |     0.334 |           | 
##                      |     0.258 |     0.080 |           | 
## ---------------------|-----------|-----------|-----------|
##         Column Total |      2534 |       797 |      3331 | 
##                      |     0.761 |     0.239 |           | 
## ---------------------|-----------|-----------|-----------|
## 
## 

Each cell in the table contains five numbers. The first one N gives us the count that falls into its corresponding category. The Chi-square contribution provides us information about the cell’s contribution in the Pearson’s Chi-squared test for independence between two variables. This number measures the probability that the differences in cell counts are due to chance alone.

The number of most interest is the N/ Col Total or the counts over column total. In this case, these numbers represent the distribution for residence area type among African regions and the regions in the rest of the world. We can see the numbers are very close between African and non-African regions for each type of residence area. Therefore, we can conclude that African WHO regions do not have a difference in terms of residence area types compared to the rest of the world.

1.13 Missing Data

In the previous sections, we simply ignored the missing observations in our water dataset (na.rm = TRUE). Is this an appropriate strategy to handle incomplete data? Could the missingness pattern of those incomplete observations be important? It is possible that the arrangement of the missing observations may reflect an important factor that was not accounted for in our statistics or our models.

Missing Completely at Random (MCAR) is an assumption about the probability of missingness being equal for all cases; Missing at Random (MAR) assumes the probability of missingness has a known but random mechanism (e.g., different rates for different groups); Missing not at Random (MNAR) suggest a missingness mechanism linked to the values of predictors and/or response, e.g., some participants may drop out of a drug trial when they have side-effects.

There are a number of strategies to impute missing data. The expectation maximization (EM) algorithm provides one example for handling missing data. The SOCR EM tutorial, activity, and documentations provides the theory, applications and practice for effective (multidimensional) EM parameter estimation.

The simplest way to handle incomplete data is to substitute each missing value with its (feature or column) average. When the missingness proportion is small, the effect of substituting the means for the missing values will have little effect on the mean, variance, or other important statistics of the data. Also, this will preserve those non-missing values of the same observation or row.

m1<-mean(water$improved_water, na.rm = T)
m2<-mean(water$sanitation_facilities, na.rm = T)
water_imp<-water
for(i in 1:3331){
  if(is.na(water_imp$improved_water[i])){
    water_imp$improved_water[i] <- m1
  }
  if(is.na(water_imp$sanitation_facilities[i])){
    water_imp$sanitation_facilities[i] <- m2
  }
}
summary(water_imp)
##       year         region            country          residence_area    
##  Min.   :1990   Length:3331        Length:3331        Length:3331       
##  1st Qu.:1995   Class :character   Class :character   Class :character  
##  Median :2005   Mode  :character   Mode  :character   Mode  :character  
##  Mean   :2002                                                           
##  3rd Qu.:2010                                                           
##  Max.   :2012                                                           
##  improved_water  sanitation_facilities   africa       
##  Min.   :  3.0   Min.   :  0.00        Mode :logical  
##  1st Qu.: 77.0   1st Qu.: 44.00        FALSE:2534     
##  Median : 93.0   Median : 79.00        TRUE :797      
##  Mean   : 84.9   Mean   : 68.87                       
##  3rd Qu.: 99.0   3rd Qu.: 97.00                       
##  Max.   :100.0   Max.   :100.00

A more sophisticated way of resolving missing data is to use a model (e.g., linear regression) to predict the missing feature and impute its missing values. This is called the predictive mean matching approach. This method is good for data with multivariate normality. However, a disadvantage of it is that it can only predict one value at a time, which is very time consuming. Also, the multivariate normality assumption might not be satisfied and there may be important multivariate relations that are not accounted for. We are using the mi package for the predictive mean matching procedure.

Let’s install the mi package first.

# install.packages("mi")
library(mi)

Then we need to get the missing information matrix. We are using the imputation method pmm(predictive mean matching approach) for both missing variables.

mdf<-missing_data.frame(water)
head(mdf)
##   year region      country residence_area improved_water sanitation_facilities
## 1 1990 Africa      Algeria          Rural             88                    77
## 2 1990 Africa       Angola          Rural             42                     7
## 3 1990 Africa        Benin          Rural             49                     0
## 4 1990 Africa     Botswana          Rural             86                    22
## 5 1990 Africa Burkina Faso          Rural             39                     2
## 6 1990 Africa      Burundi          Rural             67                    42
##   africa missing_improved_water missing_sanitation_facilities
## 1   TRUE                  FALSE                         FALSE
## 2   TRUE                  FALSE                         FALSE
## 3   TRUE                  FALSE                         FALSE
## 4   TRUE                  FALSE                         FALSE
## 5   TRUE                  FALSE                         FALSE
## 6   TRUE                  FALSE                         FALSE
show(mdf)
## Object of class missing_data.frame with 3331 observations on 7 variables
## 
## There are 3 missing data patterns
## 
## Append '@patterns' to this missing_data.frame to access the corresponding pattern for every observation or perhaps use table()
## 
##                                        type missing method  model
## year                             continuous       0   <NA>   <NA>
## region                unordered-categorical       0   <NA>   <NA>
## country               unordered-categorical       0   <NA>   <NA>
## residence_area        unordered-categorical       0   <NA>   <NA>
## improved_water                   continuous      32    ppd linear
## sanitation_facilities            continuous     135    ppd linear
## africa                               binary       0   <NA>   <NA>
## 
##                         family     link transformation
## year                      <NA>     <NA>    standardize
## region                    <NA>     <NA>           <NA>
## country                   <NA>     <NA>           <NA>
## residence_area            <NA>     <NA>           <NA>
## improved_water        gaussian identity    standardize
## sanitation_facilities gaussian identity    standardize
## africa                    <NA>     <NA>           <NA>
mdf<-change(mdf, y="improved_water", what = "imputation_method", to="pmm")
mdf<-change(mdf, y="sanitation_facilities", what = "imputation_method", to="pmm")
  • Notes:

  • Converting the input data.frame to a missing_data.frame allows us to include in the DF enhanced metadata about each variable, which is essential for the subsequent modeling, interpretation and imputation of the initial missing data.

  • show() displays all missing variables and their class-labels (e.g., continuous), along with meta-data. The missing_data.frame constructor suggests the most appropriate classes for each missing variable, however, the user often needs to correct, modify or change these meta-data, using change().

  • Use the change() function to change/correct many meta-data in the constructed missing_data.frame object which are incorrect when using show(mfd).

  • To get a sense of the raw data, look at the summary, image, or hist of the missing_data.frame.

  • The mi vignettes provide many useful examples of handling missing data.

We can perform the initial imputation. Here we imputed three times, which will create three different (complete) datasets, three chains, with slightly different imputed values.

imputations <- mi(mdf, n.iter=10, n.chains=3, verbose=T)

Next, we need to extract several multiply imputed data.frames from imputations object. Finally, we can compare the summary stats between the original dataset and the imputed datasets.

library(mi)
data.frames <- complete(imputations, 3)
summary(water)
##       year         region            country          residence_area    
##  Min.   :1990   Length:3331        Length:3331        Length:3331       
##  1st Qu.:1995   Class :character   Class :character   Class :character  
##  Median :2005   Mode  :character   Mode  :character   Mode  :character  
##  Mean   :2002                                                           
##  3rd Qu.:2010                                                           
##  Max.   :2012                                                           
##                                                                         
##  improved_water  sanitation_facilities   africa       
##  Min.   :  3.0   Min.   :  0.00        Mode :logical  
##  1st Qu.: 77.0   1st Qu.: 42.00        FALSE:2534     
##  Median : 93.0   Median : 81.00        TRUE :797      
##  Mean   : 84.9   Mean   : 68.87                       
##  3rd Qu.: 99.0   3rd Qu.: 97.00                       
##  Max.   :100.0   Max.   :100.00                       
##  NA's   :32      NA's   :135
summary(data.frames[[1]])
##       year                        region                   country    
##  Min.   :1990   Africa               :797   Albania            :  18  
##  1st Qu.:1995   Americas             :613   Algeria            :  18  
##  Median :2005   Eastern Mediterranean:373   Andorra            :  18  
##  Mean   :2002   Europe               :910   Angola             :  18  
##  3rd Qu.:2010   South-East Asia      :191   Antigua and Barbuda:  18  
##  Max.   :2012   Western Pacific      :447   Argentina          :  18  
##                                             (Other)            :3223  
##  residence_area improved_water   sanitation_facilities   africa    
##  Rural:1095     Min.   :  3.00   Min.   :  0.00        FALSE:2534  
##  Total:1109     1st Qu.: 77.00   1st Qu.: 43.00        TRUE : 797  
##  Urban:1127     Median : 93.00   Median : 81.00                    
##                 Mean   : 84.83   Mean   : 69.28                    
##                 3rd Qu.: 99.00   3rd Qu.: 97.00                    
##                 Max.   :100.00   Max.   :100.00                    
##                                                                    
##  missing_improved_water missing_sanitation_facilities
##  Mode :logical          Mode :logical                
##  FALSE:3299             FALSE:3196                   
##  TRUE :32               TRUE :135                    
##                                                      
##                                                      
##                                                      
## 
mySummary <- lapply(data.frames, summary)
mySummary$`chain:1`   # report just the summary of the first chain.
##       year                        region                   country    
##  Min.   :1990   Africa               :797   Albania            :  18  
##  1st Qu.:1995   Americas             :613   Algeria            :  18  
##  Median :2005   Eastern Mediterranean:373   Andorra            :  18  
##  Mean   :2002   Europe               :910   Angola             :  18  
##  3rd Qu.:2010   South-East Asia      :191   Antigua and Barbuda:  18  
##  Max.   :2012   Western Pacific      :447   Argentina          :  18  
##                                             (Other)            :3223  
##  residence_area improved_water   sanitation_facilities   africa    
##  Rural:1095     Min.   :  3.00   Min.   :  0.00        FALSE:2534  
##  Total:1109     1st Qu.: 77.00   1st Qu.: 43.00        TRUE : 797  
##  Urban:1127     Median : 93.00   Median : 81.00                    
##                 Mean   : 84.83   Mean   : 69.28                    
##                 3rd Qu.: 99.00   3rd Qu.: 97.00                    
##                 Max.   :100.00   Max.   :100.00                    
##                                                                    
##  missing_improved_water missing_sanitation_facilities
##  Mode :logical          Mode :logical                
##  FALSE:3299             FALSE:3196                   
##  TRUE :32               TRUE :135                    
##                                                      
##                                                      
##                                                      
## 

This is just a brief introduction for handling incomplete datasets. In later chapters, we will discuss more about missing data with different imputation methods and how to evaluate the complete imputed results.

1.13.1 Overview of R Missing Data Packages

Imputation is the process of filling in missing values with plausible estimates. Modern imputation methods go far beyond simple mean or median replacement. Examples of recent and powerful R packages for handling missing data are included below.

  • naniar: a “tidy” package to missing data. It integrates well with the tidyverse and provides excellent functions and visualizations for exploring missing data patterns, such as gg_miss_var() to see the proportion of missing values in each variable.

  • VIM: Visualization and Imputation of Missing Values: a classiccal missing data package offering a variety of graphical representations of missing data, including aggr() for aggregated missingness patterns and matrix plots to visualize missingness across variables.

  • mice: Another classical multiple imputation R package. The Multivariate Imputation by Chained Equations package is a comprehensive and flexible framework for creating multiple plausible imputed datasets, especially under the assumption of “Missing at Random” (MAR).

  • missForest: For non-parametric imputation, it uses a random forest algorithm to predict missing values based on the other variables in the dataset. This is particularly effective for mixed-type data (e.g., a mix of continuous, categorical, and ordinal variables) and for non-linear relationships.

  • missRanger: provides a fast and efficient alternative to missForest. It uses the ranger package for a speedy implementation of the random forest algorithm for imputation. It also includes the option for predictive mean matching, which helps to retain the original data distribution and avoid imputing values not present in the original dataset.

  • imputomics: A more specialized package designed for missing value imputation in metabolomics data. It serves as a wrapper for over 40 different imputation algorithms and includes a web application for user-friendly access. While specialized, it demonstrates the cutting edge of imputation for specific scientific domains.

  • inlamemi: Integrated Nested Laplace Approximations (INLA) handles missing data and measurement error in a Bayesian framework. It’s particularly useful for hierarchical models and is a good option for analysts who are less experienced with the R-INLA package itself, as it provides a more accessible interface.

  • missMDA: offers multiple imputation specifically for multivariate data analysis, such as Principal Component Analysis (PCA) and Multiple Correspondence Analysis (MCA). It’s a key tool when your primary goal is dimension reduction or visualization of high-dimensional data with missing values.

  • jomo: designed for multilevel multiple imputation, which is necessary when your data has a hierarchical or clustered structure (e.g., students nested within schools). It’s an important tool for avoiding bias in multilevel models with missing data.

  • rMIDAS: multiple imputation using denoising auto-encoders, a modern machine learning approach.

Multiple imputation methods like mice work best for data that is “Missing at Random” (MAR), while machine learning-based approaches like missForest can handle more complex relationships.

1.13.2 Simulate some real multivariate data

Suppose we would like to generate a synthetic dataset: \[sim\_data=\{y, x_1, x_2, x_3, x_4, x_5, x_6, x_7, x_8, x_9, x_{10}\}.\]

Then, we can introduce a method that takes a dataset and a desired proportion of missingness and wipes out the same proportion of the data, i.e., introduces random patterns of missingness. Note that there are already R functions that automate the introduction of missingness, e.g., missForest::prodNA(), however writing such a method from scratch is also useful.

set.seed(123)
# create MCAR missing-data generator
create.missing <- function (data, pct.mis = 10) 
{
    n <- nrow(data)
    J <- ncol(data)
    if (length(pct.mis) == 1) {
        if(pct.mis>= 0 & pct.mis <=100) {
            n.mis <- rep((n * (pct.mis/100)), J)
        }
        else {
          warning("Percent missing values should be an integer between 0 and 100! Exiting"); break
        }
      }
    else {
        if (length(pct.mis) < J) 
            stop("The length of the missing-vector is not equal to the number of columns in the data! Exiting!")
            n.mis <- n * (pct.mis/100)
    }
    for (i in 1:ncol(data)) {
        if (n.mis[i] == 0) { # if the column has no missing values, do nothing
            data[, i] <- data[, i]
        }
        else {
            data[sample(1:n, n.mis[i], replace = FALSE), i] <- NA
              # For each given column (i), sample the row indices (1:n), 
              # a number of indices to replace as "missing", n.mis[i], "NA", 
          # without replacement
        }
    }
    return(as.data.frame(data))
}

Next, let’s synthetically generate (simulate) \(1,000\) cases including all 11 features in the data (\(\{y, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10\}\)).

n <- 1000; u1 <- rbinom(n, 1, .5); v1 <- log(rnorm(n, 5, 1)); x1 <- u1*exp(v1)
u2 <- rbinom(n, 1, .5); v2 <- log(rnorm(n, 5, 1)); x2 <- u2*exp(v2)
x3 <- rbinom(n, 1, prob=0.45); x4 <- ordered(rep(seq(1, 5), n)[sample(1:n, n)])
x5 <- rep(letters[1:10], n)[sample(1:n, n)]; x6 <- trunc(runif(n, 1, 10))
x7 <- rnorm(n); x8 <- factor(rep(seq(1, 10), n)[sample(1:n, n)])
x9 <- runif(n, 0.1, .99); x10 <- rpois(n, 4)
y <- x1 + x2 + x7 + x9 + rnorm(n)

# package the simulated data as a data frame object
sim_data <- cbind.data.frame(y, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10)

# randomly create missing values
sim_data_30pct_missing <- create.missing(sim_data, pct.mis=30); 
# head(sim_data_30pct_missing); summary(sim_data_30pct_missing)

# install.packages("DT")
library("DT")
library(dplyr)
df_raw <- sim_data %>%  mutate_if(is.numeric, round, digits = 2)
datatable(df_raw)
df_miss <- sim_data_30pct_missing %>%  mutate_if(is.numeric, round, digits = 2)
datatable(df_miss)
# install.packages("mi")
# install.packages("betareg")
library("betareg"); library("mi")

# get show the missing information matrix           
mdf <- missing_data.frame(sim_data_30pct_missing) 
# show(mdf)
df_mdf <- as.data.frame(mdf) %>%  mutate_if(is.numeric, round, digits = 2)
datatable(df_mdf)
# mdf@patterns   # to get the textual missing pattern
image(mdf)   # remember the visual pattern of this MCAR

# df_img <- df_mdf %>% mutate_if(is.factor, as.character) %>%  replace(is.character(.), 1)  %>%  replace(is.na(.), 0)
# # df_img <- df_mdf %>%  replace(is.character(.), 1)  %>%  replace(is.na(.), 0)
# df_img [1:10,1:10]
# df_img[is.character(df_img)] <- 1

In the missing data plot above, missing values are illustrated as black segments in the case-by-feature bivariate chart. The hot colormap (17-level) represents the normalized values of the corresponding feature-index pairs, see the mi::image() documentation. Also, test the order, cluster and grayscale options, e.g., image(mdf, x.order = T, clustered = F, grayscale =T).

The histogram plots display the distributions of:

  • The observed data (in blue color),
  • The imputed data (in red color), and
  • The completed values (observed plus imputed, in gray color).
# Next try to impute the missing values.

# Get the Graph Parameters (plotting canvas/margins)
# set to plot the histograms for the 3 imputation chains
# mfcol=c(nr, nc). Subsequent histograms are drawn as nr-by-nc arrays on the graphics device by columns (mfcol), or rows (mfrow)
# oma
# oma=c(bottom, left, top, right) giving the size of the outer margins in lines of text
# mar=c(bottom, left, top, right) gives the number of lines of margin to be specified on the four sides of the plot.
# tcl=length of tick marks as a fraction of the height of a line of text (default=0.5)
par(mfcol=c(5, 5), oma=c(1, 1, 0, 0), mar=c(1, 1, 1, 0), tcl=-0.1, mgp=c(0, 0, 0))

# Note to get verbose output-report, parallel must be OFF: parallel=FALSE, verbose=TRUE
imputations <- mi(sim_data_30pct_missing, n.iter=5, n.chains=3, verbose=TRUE)
hist(imputations)

# Extracts several multiply imputed data.frames from "imputations" object
data.frames <- complete(imputations, 3)

# compare the 3 objects, sim_data, sim_data_30pct_missing, and imputed chain1
# datatable(sim_data, caption = htmltools::tags$caption(
#      style = 'caption-side: bottom; text-align: center;','Table: Initial sim_data'))
df_miss <- sim_data_30pct_missing %>%  mutate_if(is.numeric, round, digits = 2)
datatable(df_miss, caption = htmltools::tags$caption(
     style = 'caption-side: bottom; text-align: center;','Table: Initial sim_data'))
# datatable(sim_data_30pct_missing, caption = htmltools::tags$caption(
#      style = 'caption-side: bottom; text-align: center;', 'Table: sim_data_30pct_missing'))
df_miss30pct <- sim_data_30pct_missing %>%  mutate_if(is.numeric, round, digits = 2)
datatable(df_miss30pct, caption = htmltools::tags$caption(
     style = 'caption-side: bottom; text-align: center;', 
     'Table: sim_data 30% Missing'))
df_chain1 <- data.frames[[1]] %>%  mutate_if(is.numeric, round, digits = 2)
datatable(df_chain1, caption = htmltools::tags$caption(
     style = 'caption-side: bottom; text-align: center;', 'Table: Imputed data (chain 1)'))
# Compare the summary stats for the original data (prior to introducing missing
# values) with missing data and the re-completed data following imputation
# summary(sim_data)
datatable(data.frame(t(as.matrix(unclass(summary(sim_data)))), check.names = FALSE, stringsAsFactors = FALSE), caption = htmltools::tags$caption(
     style = 'caption-side: bottom; text-align: center;', 
     'Table: summary(sim_data)'))
mySummary <- lapply(data.frames, summary)
datatable(data.frame(t(as.matrix(unclass(mySummary$`chain:1`))), check.names = FALSE, stringsAsFactors = FALSE), caption = htmltools::tags$caption(
     style = 'caption-side: bottom; text-align: center;', 
     'Table: Imputed data: summary(chain:1)'))

Let’s check imputation convergence (details provided below).

round(mipply(imputations, mean, to.matrix = TRUE), 3)
##             chain:1 chain:2 chain:3
## y             0.005   0.006   0.008
## x1            0.007  -0.005   0.000
## x2           -0.014  -0.008  -0.023
## x3            1.413   1.405   1.420
## x4            3.001   3.012   3.015
## x5            5.635   5.650   5.437
## x6           -0.010   0.001  -0.004
## x7            0.023   0.014  -0.003
## x8            5.527   5.528   5.538
## x9            0.529   0.540   0.534
## x10           0.014   0.016   0.008
## missing_y     0.300   0.300   0.300
## missing_x1    0.300   0.300   0.300
## missing_x2    0.300   0.300   0.300
## missing_x3    0.300   0.300   0.300
## missing_x4    0.300   0.300   0.300
## missing_x5    0.300   0.300   0.300
## missing_x6    0.300   0.300   0.300
## missing_x7    0.300   0.300   0.300
## missing_x8    0.300   0.300   0.300
## missing_x9    0.300   0.300   0.300
## missing_x10   0.300   0.300   0.300
Rhats(imputations, statistic = "moments") # assess the convergence of MI algorithm
##    mean_y   mean_x1   mean_x2   mean_x3   mean_x4   mean_x5   mean_x6   mean_x7 
## 1.1373913 1.3822138 0.9803310 0.9568204 0.9615887 0.9387548 0.9148436 0.9063065 
##   mean_x8   mean_x9  mean_x10      sd_y     sd_x1     sd_x2     sd_x3     sd_x4 
## 1.1526740 1.0563110 0.9176943 1.0175290 1.0090472 1.5802533 0.9422021 1.0356559 
##     sd_x5     sd_x6     sd_x7     sd_x8     sd_x9    sd_x10 
## 1.0196262 0.9601113 1.0293472 0.9819775 0.9416442 1.0333277
plot(imputations); hist(imputations); image(imputations); summary(imputations)

## $y
## $y$is_missing
## missing
## FALSE  TRUE 
##   700   300 
## 
## $y$imputed
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -1.49285 -0.31246  0.01510  0.02139  0.36033  1.42913 
## 
## $y$observed
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -1.194077 -0.395670  0.001664  0.000000  0.347725  1.423707 
## 
## 
## $x1
## $x1$is_missing
## missing
## FALSE  TRUE 
##   700   300 
## 
## $x1$imputed
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -1.525972 -0.353452 -0.011120  0.002521  0.379387  1.836038 
## 
## $x1$observed
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -0.4805 -0.4805 -0.4805  0.0000  0.4817  1.1365 
## 
## 
## $x2
## $x2$is_missing
## missing
## FALSE  TRUE 
##   700   300 
## 
## $x2$imputed
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -1.55946 -0.42216 -0.05924 -0.05025  0.30196  1.50480 
## 
## $x2$observed
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -0.48420 -0.48420  0.01595  0.00000  0.47667  1.14717 
## 
## 
## $x3
## $x3$crosstab
##    
##     observed imputed
##   0     1215     547
##   1      885     353
## 
## 
## $x4
## $x4$crosstab
##    
##     observed imputed
##   1      408     172
##   2      429     184
##   3      438     176
##   4      411     174
##   5      414     194
## 
## 
## $x5
## $x5$crosstab
##    
##     observed imputed
##   a      186      74
##   b      210      87
##   c      195     102
##   d      231     104
##   e      198      95
##   f      219      72
##   g      210      75
##   h      219     110
##   i      216     115
##   j      216      66
## 
## 
## $x6
## $x6$is_missing
## missing
## FALSE  TRUE 
##   700   300 
## 
## $x6$imputed
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -1.59274 -0.38330 -0.02166 -0.01504  0.35464  1.65916 
## 
## $x6$observed
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -0.81073 -0.41317 -0.01562  0.00000  0.38194  0.77949 
## 
## 
## $x7
## $x7$is_missing
## missing
## FALSE  TRUE 
##   700   300 
## 
## $x7$imputed
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -1.59511 -0.31523  0.04092  0.03773  0.41539  1.54865 
## 
## $x7$observed
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -1.576156 -0.336879  0.003784  0.000000  0.369404  1.773350 
## 
## 
## $x8
## $x8$crosstab
##     
##      observed imputed
##   1       213      81
##   2       222      79
##   3       186      66
##   4       210      93
##   5       228     107
##   6       219      85
##   7       210     105
##   8       207     123
##   9       198      77
##   10      207      84
## 
## 
## $x9
## $x9$is_missing
## missing
## FALSE  TRUE 
##   700   300 
## 
## $x9$imputed
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.004661 0.360787 0.554831 0.548610 0.751360 0.998730 
## 
## $x9$observed
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1010  0.3077  0.5100  0.5286  0.7681  0.9881 
## 
## 
## $x10
## $x10$is_missing
## missing
## FALSE  TRUE 
##   700   300 
## 
## $x10$imputed
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -1.64825 -0.34635  0.03983  0.04177  0.42710  1.63302 
## 
## $x10$observed
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -0.981263 -0.494111 -0.006959  0.000000  0.236617  1.698072

Finally, pool over the \(m = 3\) completed datasets when we fit the “model”. Pool from across the 3 chains - in order to estimate a linear regression model.

model_results <- pool(y ~ x1+x2+x3+x4+x5+x6+x7+x8+x9+x10, data=imputations,  m=3)
display (model_results); summary (model_results)  
## bayesglm(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + 
##     x9 + x10, data = imputations, m = 3)
##             coef.est coef.se
## (Intercept) -0.68     0.59  
## x1           0.93     0.03  
## x2           0.92     0.02  
## x31          0.20     0.23  
## x4.L        -0.15     0.15  
## x4.Q        -0.05     0.21  
## x4.C         0.02     0.16  
## x4^4        -0.03     0.35  
## x5b          0.03     0.22  
## x5c          0.26     0.42  
## x5d          0.41     0.47  
## x5e          0.54     0.27  
## x5f          0.57     0.65  
## x5g          0.04     0.66  
## x5h          0.59     0.45  
## x5i          0.59     0.32  
## x5j         -0.24     0.28  
## x6           0.01     0.03  
## x7           1.04     0.14  
## x82          0.41     0.70  
## x83          0.09     0.45  
## x84          0.37     0.74  
## x85          0.26     0.57  
## x86         -0.14     0.29  
## x87         -0.28     0.48  
## x88          0.30     0.38  
## x89         -0.06     0.35  
## x810        -0.43     1.00  
## x9           1.14     0.41  
## x10          0.11     0.07  
## n = 970, k = 30
## residual deviance = 1878.2, null deviance = 14891.0 (difference = 13012.9)
## overdispersion parameter = 1.9
## residual sd is sqrt(overdispersion) = 1.39
## 
## Call:
## pool(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + 
##     x10, data = imputations, m = 3)
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.67622    0.58950  -1.147   0.3281    
## x1           0.93147    0.03044  30.601 3.65e-06 ***
## x2           0.92049    0.01902  48.391  < 2e-16 ***
## x31          0.20006    0.23223   0.861   0.4570    
## x4.L        -0.15290    0.15309  -0.999   0.3541    
## x4.Q        -0.05366    0.21240  -0.253   0.8156    
## x4.C         0.02152    0.15761   0.137   0.8962    
## x4^4        -0.03186    0.35333  -0.090   0.9354    
## x5b          0.02664    0.22497   0.118   0.9060    
## x5c          0.26188    0.41780   0.627   0.5688    
## x5d          0.41229    0.46559   0.886   0.4408    
## x5e          0.53762    0.26555   2.025   0.0615 .  
## x5f          0.57082    0.65087   0.877   0.4583    
## x5g          0.04444    0.66095   0.067   0.9515    
## x5h          0.58780    0.44970   1.307   0.2785    
## x5i          0.58983    0.32332   1.824   0.1196    
## x5j         -0.24061    0.27827  -0.865   0.4049    
## x6           0.01476    0.02525   0.585   0.5751    
## x7           1.04203    0.14333   7.270   0.0114 *  
## x82          0.40871    0.70257   0.582   0.6126    
## x83          0.08514    0.45155   0.189   0.8616    
## x84          0.36829    0.74270   0.496   0.6641    
## x85          0.25868    0.57355   0.451   0.6880    
## x86         -0.13945    0.28689  -0.486   0.6405    
## x87         -0.27825    0.48116  -0.578   0.6053    
## x88          0.29517    0.38024   0.776   0.4837    
## x89         -0.06421    0.34646  -0.185   0.8606    
## x810        -0.42747    1.00223  -0.427   0.7101    
## x9           1.14356    0.40965   2.792   0.0676 .  
## x10          0.10595    0.06793   1.560   0.2375    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 1.936243)
## 
##     Null deviance: 14891.0  on 999  degrees of freedom
## Residual deviance:  1878.2  on 970  degrees of freedom
## AIC: 3528.3
## 
## Number of Fisher Scoring iterations: 7
# Report the summaries of the imputations
data.frames <- complete(imputations, 3)     # extract the first 3 chains
mySummary <-lapply(data.frames, summary)

datatable(data.frame(t(as.matrix(unclass(mySummary$`chain:1`))), check.names = FALSE, stringsAsFactors = FALSE), caption = htmltools::tags$caption(
     style = 'caption-side: bottom; text-align: center;', 
     'Table: Imputed data: summary(chain:1)'))
datatable(data.frame(t(as.matrix(unclass(mySummary$`chain:2`))), check.names = FALSE, stringsAsFactors = FALSE), caption = htmltools::tags$caption(
     style = 'caption-side: bottom; text-align: center;', 
     'Table: Imputed data: summary(chain:2)'))
datatable(data.frame(t(as.matrix(unclass(mySummary$`chain:3`))), check.names = FALSE, stringsAsFactors = FALSE), caption = htmltools::tags$caption(
     style = 'caption-side: bottom; text-align: center;', 
     'Table: Imputed data: summary(chain:3)'))
coef(summary(model_results))[, 1:2]  # get the model coefficients and their SE's  
##                Estimate Std. Error
## (Intercept) -0.67622094 0.58949544
## x1           0.93146870 0.03043923
## x2           0.92048806 0.01902193
## x31          0.20006237 0.23223298
## x4.L        -0.15289599 0.15308885
## x4.Q        -0.05365643 0.21240082
## x4.C         0.02152036 0.15760678
## x4^4        -0.03185678 0.35333003
## x5b          0.02663712 0.22496747
## x5c          0.26187684 0.41779723
## x5d          0.41229107 0.46559162
## x5e          0.53762308 0.26554726
## x5f          0.57081833 0.65087265
## x5g          0.04443623 0.66095146
## x5h          0.58780289 0.44969556
## x5i          0.58982952 0.32331663
## x5j         -0.24061068 0.27827419
## x6           0.01476263 0.02524864
## x7           1.04202828 0.14333219
## x82          0.40871269 0.70256680
## x83          0.08514471 0.45154512
## x84          0.36829441 0.74269647
## x85          0.25867505 0.57354905
## x86         -0.13944675 0.28689344
## x87         -0.27825103 0.48115701
## x88          0.29516760 0.38024087
## x89         -0.06421226 0.34645818
## x810        -0.42746615 1.00223277
## x9           1.14355929 0.40965449
## x10          0.10594938 0.06792528
library("lattice")
densityplot(y ~ x1 + x2, data=imputations)

# plot_ly(imputations@data$`chain:1`, x=~(x1+x2), y=~density(y))
# To compare the density of observed data and imputed data --
# these should be similar (though not identical) under MAR assumption

Notes:

  • In general, it is recommended to generate multiple imputation chains and then analyze the data (e.g., estimate the model coefficients, obtain inference, compute likelihoods, etc.). Pooling the analytics across all chains accounts for between-chain as well as within-chain variability, Rubin’s rule.
  • When deciding on how many chains to compute, a general rule is to compute \(m\) chains if the rate of incomplete cases in the dataset is about \(m \%\), i.e., 10-chains when 10% of cases are incomplete, White et al.,2011.
  • For categorical features, e.g., binary predictors like \(x_3\), the display() and summary() functions will report coefficient estimates for each (category) level, relative to the base level.

1.13.3 TBI Data Example

Next, we will see an example using the traumatic brain injury (TBI) dataset. More information about the clinical assessment scores (e.g., EGOS, GCS) is available in this publication (DOI: 10.1080/02699050701727460).

# Load the (raw) data from the table into a plain text file "08_EpiBioSData_Incomplete.csv"
TBI_Data <- read.csv("https://umich.instructure.com/files/720782/download?download_frd=1", na.strings=c("", ".", "NA"))    ## 1. read in data
summary(TBI_Data)
##        id             age            sex             mechanism        
##  Min.   : 1.00   Min.   :16.00   Length:46          Length:46         
##  1st Qu.:12.25   1st Qu.:23.00   Class :character   Class :character  
##  Median :23.50   Median :33.00   Mode  :character   Mode  :character  
##  Mean   :23.50   Mean   :36.89                                        
##  3rd Qu.:34.75   3rd Qu.:47.25                                        
##  Max.   :46.00   Max.   :83.00                                        
##                                                                       
##    field.gcs      er.gcs          icu.gcs         worst.gcs       X6m.gose    
##  Min.   : 3   Min.   : 3.000   Min.   : 0.000   Min.   : 0.0   Min.   :2.000  
##  1st Qu.: 3   1st Qu.: 4.000   1st Qu.: 3.000   1st Qu.: 3.0   1st Qu.:3.000  
##  Median : 7   Median : 7.500   Median : 6.000   Median : 3.0   Median :5.000  
##  Mean   : 8   Mean   : 8.182   Mean   : 6.378   Mean   : 5.4   Mean   :4.805  
##  3rd Qu.:12   3rd Qu.:12.250   3rd Qu.: 8.000   3rd Qu.: 7.0   3rd Qu.:6.000  
##  Max.   :15   Max.   :15.000   Max.   :14.000   Max.   :14.0   Max.   :8.000  
##  NA's   :2    NA's   :2        NA's   :1        NA's   :1      NA's   :5      
##    X2013.gose       skull.fx       temp.injury       surgery      
##  Min.   :2.000   Min.   :0.0000   Min.   :0.000   Min.   :0.0000  
##  1st Qu.:5.000   1st Qu.:0.0000   1st Qu.:0.000   1st Qu.:0.0000  
##  Median :7.000   Median :1.0000   Median :1.000   Median :1.0000  
##  Mean   :5.804   Mean   :0.6087   Mean   :0.587   Mean   :0.6304  
##  3rd Qu.:7.000   3rd Qu.:1.0000   3rd Qu.:1.000   3rd Qu.:1.0000  
##  Max.   :8.000   Max.   :1.0000   Max.   :1.000   Max.   :1.0000  
##                                                                   
##    spikes.hr           min.hr           max.hr           acute.sz     
##  Min.   :  1.280   Min.   : 0.000   Min.   :  12.00   Min.   :0.0000  
##  1st Qu.:  5.357   1st Qu.: 0.000   1st Qu.:  35.25   1st Qu.:0.0000  
##  Median : 18.170   Median : 0.000   Median :  97.50   Median :0.0000  
##  Mean   : 52.872   Mean   : 3.571   Mean   : 241.89   Mean   :0.1739  
##  3rd Qu.: 57.227   3rd Qu.: 0.000   3rd Qu.: 312.75   3rd Qu.:0.0000  
##  Max.   :294.000   Max.   :42.000   Max.   :1199.00   Max.   :1.0000  
##  NA's   :18        NA's   :18       NA's   :18                        
##     late.sz          ever.sz     
##  Min.   :0.0000   Min.   :0.000  
##  1st Qu.:0.0000   1st Qu.:0.000  
##  Median :1.0000   Median :1.000  
##  Mean   :0.5652   Mean   :0.587  
##  3rd Qu.:1.0000   3rd Qu.:1.000  
##  Max.   :1.0000   Max.   :1.000  
## 
# Get information matrix of the data
# 2. create an object of class "missing_data.frame" from the data.frame TBI_data
# Convert to a missing_data.frame
# library("betareg"); library("mi")         
mdf <- missing_data.frame(TBI_Data) # warnings about missingness patterns
## NOTE: The following pairs of variables appear to have the same missingness pattern.
##  Please verify whether they are in fact logically distinct variables.
##      [,1]      [,2]       
## [1,] "icu.gcs" "worst.gcs"
datatable(mdf)
# mdf@patterns; 
image(mdf)

# 3. get description of the "family", "imputation_method", "size", "transformation", "type", "link", or "model" of each incomplete variable
# show(mdf)

# 4. change things: mi::change() method changes the family, imputation method,
# size, type, and so forth of a missing variable. It's called 
# before calling mi to affect how the conditional expectation of each 
# missing variable is modeled.

mdf <- change(mdf, y = "spikes.hr", what = "transformation", to = "identity")
# The "to" choices include "identity" = no transformation, "standardize" = standardization, "log" = natural logarithm transformation, "logshift" = log(y + a) transformation, where a is a small constant, or "sqrt" = square-root variable transformation. Changing the transformation will correspondingly change the inverse transformation.
# 5. examine missingness patterns
summary(mdf); hist(mdf); 
##        id             age            sex             mechanism        
##  Min.   : 1.00   Min.   :16.00   Length:46          Length:46         
##  1st Qu.:12.25   1st Qu.:23.00   Class :character   Class :character  
##  Median :23.50   Median :33.00   Mode  :character   Mode  :character  
##  Mean   :23.50   Mean   :36.89                                        
##  3rd Qu.:34.75   3rd Qu.:47.25                                        
##  Max.   :46.00   Max.   :83.00                                        
##                                                                       
##    field.gcs      er.gcs          icu.gcs         worst.gcs       X6m.gose    
##  Min.   : 3   Min.   : 3.000   Min.   : 0.000   Min.   : 0.0   Min.   :2.000  
##  1st Qu.: 3   1st Qu.: 4.000   1st Qu.: 3.000   1st Qu.: 3.0   1st Qu.:3.000  
##  Median : 7   Median : 7.500   Median : 6.000   Median : 3.0   Median :5.000  
##  Mean   : 8   Mean   : 8.182   Mean   : 6.378   Mean   : 5.4   Mean   :4.805  
##  3rd Qu.:12   3rd Qu.:12.250   3rd Qu.: 8.000   3rd Qu.: 7.0   3rd Qu.:6.000  
##  Max.   :15   Max.   :15.000   Max.   :14.000   Max.   :14.0   Max.   :8.000  
##  NA's   :2    NA's   :2        NA's   :1        NA's   :1      NA's   :5      
##    X2013.gose       skull.fx       temp.injury       surgery      
##  Min.   :2.000   Min.   :0.0000   Min.   :0.000   Min.   :0.0000  
##  1st Qu.:5.000   1st Qu.:0.0000   1st Qu.:0.000   1st Qu.:0.0000  
##  Median :7.000   Median :1.0000   Median :1.000   Median :1.0000  
##  Mean   :5.804   Mean   :0.6087   Mean   :0.587   Mean   :0.6304  
##  3rd Qu.:7.000   3rd Qu.:1.0000   3rd Qu.:1.000   3rd Qu.:1.0000  
##  Max.   :8.000   Max.   :1.0000   Max.   :1.000   Max.   :1.0000  
##                                                                   
##    spikes.hr           min.hr           max.hr           acute.sz     
##  Min.   :  1.280   Min.   : 0.000   Min.   :  12.00   Min.   :0.0000  
##  1st Qu.:  5.357   1st Qu.: 0.000   1st Qu.:  35.25   1st Qu.:0.0000  
##  Median : 18.170   Median : 0.000   Median :  97.50   Median :0.0000  
##  Mean   : 52.872   Mean   : 3.571   Mean   : 241.89   Mean   :0.1739  
##  3rd Qu.: 57.227   3rd Qu.: 0.000   3rd Qu.: 312.75   3rd Qu.:0.0000  
##  Max.   :294.000   Max.   :42.000   Max.   :1199.00   Max.   :1.0000  
##  NA's   :18        NA's   :18       NA's   :18                        
##     late.sz          ever.sz     
##  Min.   :0.0000   Min.   :0.000  
##  1st Qu.:0.0000   1st Qu.:0.000  
##  Median :1.0000   Median :1.000  
##  Mean   :0.5652   Mean   :0.587  
##  3rd Qu.:1.0000   3rd Qu.:1.000  
##  Max.   :1.0000   Max.   :1.000  
## 
image(mdf)

# 6. Perform initial imputation
imputations1 <- mi(mdf, n.iter=10, n.chains=5, verbose=TRUE)
hist(imputations1)

# 7. Extracts several multiply imputed data.frames from "imputations" object
data.frames1 <- complete(imputations1, 5)

# 8. Report a list of "summaries" for each element (imputation instance)
mySummary1 <- lapply(data.frames1, summary)

datatable(data.frame(t(as.matrix(unclass(mySummary1$`chain:1`))), check.names = FALSE, stringsAsFactors = FALSE), caption = htmltools::tags$caption(
     style = 'caption-side: bottom; text-align: center;', 
     'Table: Imputed data: summary(chain:1)'),
    extensions = 'Buttons', options = list(dom = 'Bfrtip',
    buttons = c('copy', 'csv', 'excel', 'pdf', 'print')))
datatable(data.frame(t(as.matrix(unclass(mySummary1$`chain:5`))), check.names = FALSE, stringsAsFactors = FALSE), caption = htmltools::tags$caption(
     style = 'caption-side: bottom; text-align: center;', 
     'Table: Imputed data: summary(chain:5)'),
    extensions = 'Buttons', options = list(dom = 'Bfrtip',
    buttons = c('copy', 'csv', 'excel', 'pdf', 'print')))
# 8.a. To cast the imputed numbers as integers (not necessary, but may be useful)
indx <- sapply(data.frames1[[5]], is.numeric)  # get the indices of numeric columns
data.frames1[[5]][indx] <- lapply(data.frames1[[5]][indx], function(x) as.numeric(as.integer(x)))           # cast each value as integer
# data.frames[[5]]$spikes.hr

# 9. Save results out
# write.csv(data.frames1[[5]], "C:\\Users\\IvoD\\Desktop\\TBI_MIData.csv")

# 10. Complete Data analytics functions:
# library("mi")
#lm.mi(); glm.mi(); polr.mi(); bayesglm.mi(); bayespolr.mi(); lmer.mi(); glmer.mi()

# 10.1 Define Linear Regression for multiply imputed dataset - Also see Step (12)
##linear regression for each imputed data set - 5 regression models are fit
fit_lm1 <- glm(ever.sz ~ surgery + worst.gcs + factor(sex) + age, data.frames1$`chain:1`, family = "binomial"); summary(fit_lm1); display(fit_lm1)
## 
## Call:
## glm(formula = ever.sz ~ surgery + worst.gcs + factor(sex) + age, 
##     family = "binomial", data = data.frames1$`chain:1`)
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)
## (Intercept)      0.241688   1.354865   0.178    0.858
## surgery1         0.943172   0.684361   1.378    0.168
## worst.gcs       -0.067326   0.097651  -0.689    0.491
## factor(sex)Male -0.328221   0.842607  -0.390    0.697
## age              0.004515   0.019427   0.232    0.816
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 62.371  on 45  degrees of freedom
## Residual deviance: 60.063  on 41  degrees of freedom
## AIC: 70.063
## 
## Number of Fisher Scoring iterations: 4
## glm(formula = ever.sz ~ surgery + worst.gcs + factor(sex) + age, 
##     family = "binomial", data = data.frames1$`chain:1`)
##                 coef.est coef.se
## (Intercept)      0.24     1.35  
## surgery1         0.94     0.68  
## worst.gcs       -0.07     0.10  
## factor(sex)Male -0.33     0.84  
## age              0.00     0.02  
## ---
##   n = 46, k = 5
##   residual deviance = 60.1, null deviance = 62.4 (difference = 2.3)
# Fit the appropriate model and pool the results (estimates over MI chains)
model_results <- pool(ever.sz ~ surgery + worst.gcs + factor(sex) + age, family = "binomial", data=imputations1,  m=5)
display (model_results); summary (model_results)  
## bayesglm(formula = ever.sz ~ surgery + worst.gcs + factor(sex) + 
##     age, data = imputations1, m = 5, family = "binomial")
##                 coef.est coef.se
## (Intercept)      0.39     1.27  
## surgery1         0.90     0.65  
## worst.gcs       -0.08     0.10  
## factor(sex)Male -0.31     0.76  
## age              0.00     0.02  
## n = 41, k = 5
## residual deviance = 59.6, null deviance = 62.4 (difference = 2.7)
## 
## Call:
## pool(formula = ever.sz ~ surgery + worst.gcs + factor(sex) + 
##     age, data = imputations1, m = 5, family = "binomial")
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)
## (Intercept)      0.388034   1.274033   0.305    0.761
## surgery1         0.899596   0.653402   1.377    0.169
## worst.gcs       -0.079627   0.097349  -0.818    0.413
## factor(sex)Male -0.314043   0.761536  -0.412    0.680
## age              0.002489   0.018189   0.137    0.891
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 62.371  on 45  degrees of freedom
## Residual deviance: 59.624  on 41  degrees of freedom
## AIC: 69.624
## 
## Number of Fisher Scoring iterations: 6.4
# Report the summaries of the imputations
data.frames <- complete(imputations1, 3)    # extract the first 3 chains
mySummary2 <-lapply(data.frames1, summary)

datatable(data.frame(t(as.matrix(unclass(mySummary2$`chain:1`))), check.names = FALSE, stringsAsFactors = FALSE), caption = htmltools::tags$caption(
     style = 'caption-side: bottom; text-align: center;', 
     'Table: Imputed data: summary(chain:1)'))
# 11. Validation: we now verify whether enough iterations were conducted. 
# Validation criteria demands that the mean of each completed variable should
# be similar for each of the k chains (in this case k=5).
# mipply is wrapper for sapply invoked on mi-class objects to compute the col means
round(mipply(imputations1, mean, to.matrix = TRUE), 3)
##                   chain:1 chain:2 chain:3 chain:4 chain:5
## id                 23.500  23.500  23.500  23.500  23.500
## age                 0.000   0.000   0.000   0.000   0.000
## sex                 1.804   1.804   1.804   1.804   1.804
## mechanism           4.261   4.261   4.261   4.261   4.261
## field.gcs          -0.029  -0.028   0.042   0.057   0.004
## er.gcs              0.004   0.045   0.036   0.037   0.064
## icu.gcs             0.039   0.004  -0.053  -0.037  -0.009
## worst.gcs           0.015  -0.001  -0.012  -0.020   0.014
## X6m.gose            0.035  -0.028   0.037   0.081  -0.037
## X2013.gose          0.000   0.000   0.000   0.000   0.000
## skull.fx            1.609   1.609   1.609   1.609   1.609
## temp.injury         1.587   1.587   1.587   1.587   1.587
## surgery             1.630   1.630   1.630   1.630   1.630
## spikes.hr          61.926  69.830  13.245  59.868  50.357
## min.hr              0.036  -0.091  -0.137   0.008  -0.039
## max.hr              0.067   0.055  -0.032  -0.053   0.057
## acute.sz            1.174   1.174   1.174   1.174   1.174
## late.sz             1.565   1.565   1.565   1.565   1.565
## ever.sz             1.587   1.587   1.587   1.587   1.587
## missing_field.gcs   0.043   0.043   0.043   0.043   0.043
## missing_er.gcs      0.043   0.043   0.043   0.043   0.043
## missing_icu.gcs     0.022   0.022   0.022   0.022   0.022
## missing_worst.gcs   0.022   0.022   0.022   0.022   0.022
## missing_X6m.gose    0.109   0.109   0.109   0.109   0.109
## missing_spikes.hr   0.391   0.391   0.391   0.391   0.391
## missing_min.hr      0.391   0.391   0.391   0.391   0.391
## missing_max.hr      0.391   0.391   0.391   0.391   0.391
# Rhat convergence statistics compares the variance between chains to the variance
# within chains (similar to the ANOVA F-test). 
# Rhat Values ~ 1.0 indicate likely convergence, 
# Rhat Values > 1.1 indicate that the chains should be run longer 
# (use large number of iterations)
Rhats(imputations1, statistic = "moments") # assess the convergence of MI algorithm
## mean_field.gcs    mean_er.gcs   mean_icu.gcs mean_worst.gcs  mean_X6m.gose 
##       2.089546       1.523122       2.290793       2.579924       2.826601 
## mean_spikes.hr    mean_min.hr    mean_max.hr   sd_field.gcs      sd_er.gcs 
##       2.005515       1.451672       1.484882       1.349397       1.242987 
##     sd_icu.gcs   sd_worst.gcs    sd_X6m.gose   sd_spikes.hr      sd_min.hr 
##       1.458525       1.343795       1.417582       1.370897       1.294738 
##      sd_max.hr 
##       1.963186
# When convergence is unstable, we can continue the iterations for all chains, e.g.
imputations1 <- mi(imputations1, n.iter=20) # add additional 20 iterations

# To plot the produced mi results, for all missing_variables we can generate
# a histogram of the observed, imputed, and completed data.
# We can compare of the completed data to the fitted values implied by the model
# for the completed data, by plotting binned residuals. 
# hist function works similarly as plot. 
# image function gives a sense of the missingness patterns in the data
plot(imputations1); hist(imputations1); image(imputations1)

mySummary3 <-lapply(data.frames1, summary)

datatable(data.frame(t(as.matrix(unclass(mySummary3$`chain:1`))), check.names = FALSE, stringsAsFactors = FALSE), caption = htmltools::tags$caption(
     style = 'caption-side: bottom; text-align: center;', 
     'Table: Imputed data: summary(chain:1)'))
# 12. Finally, pool over the m = 5 imputed datasets when we fit the "model"
# Pool from across the 4 chains - in order to estimate a linear regression model
# and impact of various predictors

model_results <- pool(ever.sz ~ surgery + worst.gcs + factor(sex) + age, data =  imputations1,  m =  5 ); display (model_results); summary (model_results)
## bayesglm(formula = ever.sz ~ surgery + worst.gcs + factor(sex) + 
##     age, data = imputations1, m = 5)
##                 coef.est coef.se
## (Intercept)      0.43     1.28  
## surgery1         0.92     0.66  
## worst.gcs       -0.09     0.10  
## factor(sex)Male -0.32     0.76  
## age              0.00     0.02  
## n = 41, k = 5
## residual deviance = 59.5, null deviance = 62.4 (difference = 2.9)
## 
## Call:
## pool(formula = ever.sz ~ surgery + worst.gcs + factor(sex) + 
##     age, data = imputations1, m = 5)
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)
## (Intercept)      0.43154    1.27977   0.337    0.736
## surgery1         0.92046    0.65588   1.403    0.160
## worst.gcs       -0.08604    0.09704  -0.887    0.375
## factor(sex)Male -0.32238    0.76271  -0.423    0.673
## age              0.00203    0.01825   0.111    0.911
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 62.371  on 45  degrees of freedom
## Residual deviance: 59.491  on 41  degrees of freedom
## AIC: 69.491
## 
## Number of Fisher Scoring iterations: 6.4
coef(summary(model_results))[, 1:2]  # get the model coefficients and their SE's  
##                    Estimate Std. Error
## (Intercept)      0.43154055 1.27977050
## surgery1         0.92046502 0.65587621
## worst.gcs       -0.08603569 0.09704110
## factor(sex)Male -0.32237656 0.76271187
## age              0.00203028 0.01825169

1.14 Imputation via Expectation-Maximization

Below we present the theory and practice of one specific statistical computing strategy for imputing incomplete datasets.

1.14.1 Types of missing data

  • MCAR: Data which is Missing Completely At Random has nothing systematic about which observations are missing. There is no relationship between missingness and either observed or unobserved covariates.
  • MAR: Missing At Random is weaker than MCAR. The missingness is still random, but solely due to the observed variables. For example, those from a lower socioeconomic status (SES) may be less willing to provide salary information (but we know their SES). The key is that the missingness is not due to the values which are not observed. MCAR implies MAR, but not vice-versa.
  • MNAR: If the data are Missing Not At Random, then the missingness depends on the values of the missing data. Examples include censored data, self-reported data for individuals who are heavier, who are less likely to report their weight, and response-measuring devices that can only measure values above \(0.5\), anything below that is missing.

1.14.2 General Idea of the EM algorithm

Expectation-Maximization (EM) is an iterative process involving two steps - expectation and maximization, which are applied in tandem. EM can be employed to find parameter estimates using maximum likelihood and is specifically useful when the equations determining the relations of the data-parameters cannot be directly solved. For example, a Gaussian mixture modeling assumes that each data point (\(X\)) has a corresponding latent (unobserved) variable or a missing value (\(Y\)), which may be specified as a mixture of coefficients determining the affinity of the data as a linear combination of Gaussian kernels, determined by a set of parameters (\(\theta\)), e.g., means and variance-covariances. Thus, EM estimation relies on:

  • An observed data set \(X\),
  • A set of missing (or latent) values \(Y\),
  • A parameter \(\theta\), which may be a vector of parameters,
  • A likelihood function \(L(\theta | X, Y) =p(X,Y |\theta)\), and
  • The maximum likelihood estimate (MLE) of the unknown parameter(s) \(\theta\) that is computed using the marginal likelihood of the observed data:

\[L(\theta | X) = p(X |\theta) =\int { p(X, Y |\theta)dY}.\]

Most of the time, this equation may not be directly solved, e.g., when \(Y\) is missing.

  • Expectation step (E step): computes the expected value of the log likelihood function, with respect to the conditional distribution of \(Y\) given \(X\) using the parameter estimates at the previous iteration (or at the position of initialization, for the first iteration), \(\theta_t\): \[Q ( \theta | \theta^{(t)} ) = E_{Y | X, \theta^{(t)} }[ log \left ( L(\theta | X , Y ) \right ];\]
  • Maximization step (M step): Determine the parameter, \(\theta\), that maximizes the expectation above, \[\theta^{(t+1)}=\arg\max_{\theta}Q(\theta|\theta^{(t)}).\]

This SOCR EM Activity shows the practical aspects of applying the EM algorithm. Also, in DSPA Chapter 3 we will illustrate the EM method for fitting single distribution models or (linear) mixtures of distributions to data that may represent a blend of heterogeneous observations from multiple different processes.

1.14.3 EM-based imputation

The EM algorithm is an alternative to Newton-Raphson or the method of scoring for computing MLE in cases where there are complications in calculating the MLE. It is applicable for imputing incomplete MAR data, where the missing data mechanism can be ignored and separate parameters may be estimated for each missing feature.

Complete Data: \[Z = \left(\begin{array}{cc} X \\ Y \end{array}\right), ZZ^T = \left(\begin{array}{cc} XX^T & XY^T \\ YX^T & YY^T \end{array}\right),\] where \(X\) is the observed data and \(Y\) is the missing data.

  • E-step: (Expectation) Get the expectations of \(Y\) and \(YY^T\) based on observed data, \(X\).
  • M-step: (Maximization) Maximize the conditional expectation in E-step to estimate the parameters.

Details: If \(o=obs\) and \(m=mis\) stand for observed and missing, the mean vector, \((\mu_{obs}, \mu_{mis})^T\), and the variance-covariance matrix, \(\Sigma^{(t)} = \left(\begin{array}{cc} \Sigma_{oo} & \Sigma_{om} \\ \Sigma_{mo} & \Sigma_{mm} \end{array}\right)\), are represented by:

\[\mu^{(t)} = \left(\begin{array}{cc} \mu_{obs} \\ \mu_{mis} \end{array}\right),\;\;\;\;\; \Sigma^{(t)} = \left(\begin{array}{cc} \Sigma_{oo} & \Sigma_{om} \\ \Sigma_{mo} & \Sigma_{mm} \end{array}\right)\] E-step:

\[E(Z | X) = \left(\begin{array}{cc} X \\ E(Y|X) \end{array}\right),\;\;\;\;\; E(ZZ^T|X) = \left(\begin{array}{cc} XX^T & XE(Y|X)^T \\ E(Y|X)X^T & E(YY^T|X) \end{array}\right).\]

\[E(Y | X) = \mu_{mis} + \Sigma_{mo}\Sigma_{oo}^{-1}(X - \mu_{obs}).\] \[E(YY^T|X) = (\Sigma_{mm}-\Sigma_{mo}\Sigma_{oo}^{-1}\Sigma_{om})+E(Y|X)E(Y|X)^T.\]

M-step: \[\mu^{(t+1)} = \frac{1}{n}\sum_{i=1}^nE(Z|X).\] \[\Sigma^{(t+1)} = \frac{1}{n}\sum_{i=1}^nE(ZZ^T|X) - \mu^{(t+1)}{\mu^{(t+1)}}^T.\]

1.14.4 A simple manual implementation of EM-based imputation

# --- 1. Load Necessary Libraries ---
# These libraries provide functions for data manipulation, plotting, and statistical modeling.
library(ggplot2)       # For data visualization
library(gridExtra)     # For arranging multiple plots
library(MASS)          # For multivariate normal distribution simulation (mvrnorm)
library(knitr)         # For creating dynamic reports
library(dplyr)         # For data manipulation verbs (e.g., mutate_if)
library(DT)            # For creating interactive data tables (datatable)
library(e1071)         # For the rdiscrete function


# --- 2. Simulate Multivariate Normal Data ---
# We will create a dataset that follows a multivariate normal distribution,
# which is a key assumption for this specific EM implementation.

set.seed(202227) # Set a seed for reproducibility of the random data.

# Define the parameters for the multivariate normal distribution.
# Mean vector (mu): 20 features, each with a mean of 2.
mu <- as.matrix(rep(2, 20))
# Covariance matrix (sig): A diagonal matrix, meaning the features are uncorrelated.
# The variance of each feature increases from 1 to 20.
sig <- diag(1:20)

# Generate the primary dataset: 200 observations (cases) for 20 features.
# An additional noise term is added from a standard multivariate normal distribution
# (mean=0, variance=1) to make the data more realistic.
sim_data <- mvrnorm(n = 200, mu, sig) + 
            mvrnorm(n = 200, as.matrix(rep(0, 20)), diag(rep(1, 20)))

# Keep a copy of the original, complete data for potential comparison later.
sim_data.orig <- sim_data


# --- 3. Introduce Missing Values (MCAR) ---
# To test the imputation, we will randomly remove some data points.

# We will introduce approximately 500 missing values.
# The e1071::rdiscrete function samples indices from the flattened data matrix.
# By setting `probs` to a uniform vector, each cell has an equal chance of being selected.
# This creates a "Missing Completely at Random" (MCAR) scenario.
# Note: Since sampling is with replacement, the final count of NAs might be slightly less than 500
# if some indices are selected more than once.
rand.miss <- e1071::rdiscrete(500, probs = rep(1, length(sim_data)), values = 1:length(sim_data))

# Replace the data at the randomly selected indices with NA.
sim_data[rand.miss] <- NA
# Verify the number of missing values introduced.
sum(is.na(sim_data))
## [1] 466
# Convert the matrix to a data frame for easier handling with some packages.
sim_data.df <- data.frame(sim_data)

# Display the data with missing values in an interactive table.
df_mdf <- sim_data.df %>% mutate_if(is.numeric, round, digits = 2)
datatable(df_mdf, caption = htmltools::tags$caption(
  style = 'caption-side: bottom; text-align: center;',
  'Table: Simulated Data with Missing Values (sim_data.df)'),
  extensions = 'Buttons', options = list(dom = 'Bfrtip',
  buttons = c('copy', 'csv', 'excel', 'pdf', 'print')))
# --- 4. Define the EM Imputation Algorithm ---
EM_algorithm <- function(x, tol = 0.001) {
  # --- INITIALIZATION ---
  # Create a boolean matrix to identify the locations of missing values.
  missvals <- is.na(x)
  
  # Initialize matrices for the iterative process.
  new.impute <- x
  old.impute <- x
  
  # Set up iteration counter and convergence flag.
  count.iter <- 1
  reach.tol <- 0
  
  # Initial parameter estimation (theta_0):
  # Calculate the initial covariance matrix and mean vector using only complete cases (listwise deletion).
  # This provides a starting point for the algorithm.
  sigma <- as.matrix(var(na.exclude(x)))
  mean.vec <- as.matrix(apply(na.exclude(x), 2, mean))
  
  # --- ITERATION LOOP (E and M Steps) ---
  while (reach.tol != 1) {
    
    # --- EXPECTATION (E) STEP ---
    # Iterate through each row of the dataset to fill in missing values.
    for (i in 1:nrow(x)) {
      # Identify which columns are missing for the current row `i`.
      pick.miss <- missvals[i, ]
      
      # Only perform imputation if there are missing values in this row.
      if (sum(pick.miss) > 0) {
        
        # Partition the covariance matrix and mean vector into observed and missing parts for row `i`.
        # `!pick.miss` corresponds to the observed part (oo).
        # `pick.miss` corresponds to the missing part (mis).
        
        # Calculate the inverse of the observed portion of the covariance matrix (Sigma_oo^-1).
        # A small tolerance is added for numerical stability.
        inv.S <- solve(sigma[!pick.miss, !pick.miss], tol = 1e-40)
        
        # Estimate the missing values based on the conditional expectation formula:
        # E(Y|X) = mu_mis + Sigma_mo * Sigma_oo^-1 * (X_obs - mu_obs)
        # where Y are the missing values and X are the observed values for that case.
        new.impute[i, pick.miss] <- mean.vec[pick.miss] +            # mu_mis
          sigma[pick.miss, !pick.miss] %*% inv.S %*%                 # Sigma_mo * Sigma_oo^-1
          (t(new.impute[i, !pick.miss]) - mean.vec[!pick.miss])      # (X_obs - mu_obs)
      }
    }
    
    # --- MAXIMIZATION (M) STEP ---
    # After filling in the missing values for this iteration, re-estimate the parameters
    # (mean vector and covariance matrix) using the now-complete dataset.
    # These new parameters will be used in the next E-step.
    
    # Update the covariance matrix (Sigma^(t+1)).
    sigma <- var(new.impute)
    # Update the mean vector (mu^(t+1)).
    mean.vec <- as.matrix(apply(new.impute, 2, mean))
    
    # --- CONVERGENCE CHECK ---
    # Starting from the second iteration, check if the imputed values have stabilized.
    if (count.iter > 1) {
      # Assume tolerance is reached until a value proves otherwise.
      reach.tol <- 1 
      # Compare the newly imputed matrix with the one from the previous iteration.
      # Note: A more efficient check would be `max(abs(old.impute - new.impute), na.rm = TRUE) < tol`.
      # The current loop checks every single cell individually.
      if (max(abs(old.impute - new.impute), na.rm = TRUE) > tol) {
        reach.tol <- 0 # If any change is larger than the tolerance, continue iterating.
      }
    }
    
    # Prepare for the next iteration.
    count.iter <- count.iter + 1
    old.impute <- new.impute
  }
  
  # Return the final imputed dataset once convergence is achieved.
  return(new.impute)
}

# --- 5. Run the Algorithm and Display Results ---
# Execute the EM algorithm on the dataset with missing values.
# A small tolerance is set for faster convergence.
sim_data.imputed <- EM_algorithm(sim_data.df, tol = 0.0001)

# Display the final, imputed data in an interactive table.
df_mdf_imputed <- sim_data.imputed %>% mutate_if(is.numeric, round, digits = 2)
datatable(df_mdf_imputed, caption = htmltools::tags$caption(
  style = 'caption-side: bottom; text-align: center;',
  'Table: EM-Imputed Simulated Data'),
  extensions = 'Buttons', options = list(dom = 'Bfrtip',
  buttons = c('copy', 'csv', 'excel', 'pdf', 'print')))

The above R script demonstrates a practical application of the Expectation-Maximization (EM) algorithm for imputing missing data in a multivariate dataset. The workflow consists of three main parts: data simulation, initialization, and iterative imputation.

  1. Data Simulation and Preparation: The process begins by generating a synthetic dataset of 200 observations and 20 features from a multivariate normal distribution. This is a crucial step, as the specific EM algorithm implemented assumes that the data follows this distribution. After creating the complete dataset, approximately 500 data points are removed completely at random (MCAR) to simulate a realistic scenario where data collection is incomplete.

  2. Initialization: The EM algorithm starts with an initial guess for the statistical parameters that describe the data—specifically, the mean vector (\(\mu\)) and the covariance matrix (\(\Sigma\)). These initial parameters (\(\theta^{(0)}\)) are estimated by performing listwise deletion on the incomplete dataset and calculating the mean and variance-covariance from the remaining complete rows.

Iterative E and M Steps: With the initial parameters established, the algorithm enters an iterative loop that alternates between two steps until the imputed values converge.

  1. Expectation (E) Step: In this step, the algorithm uses the current parameter estimates (\(\mu^{(t)}\) and \(\Sigma^{(t)}\)) to fill in the missing values. For each observation containing NAs, it calculates the conditional expectation of the missing variables, given the values of the observed variables in that same observation. This step produces the most likely estimates for the missing data based on the currently understood relationships (covariances) between all the features.

  2. Maximization (M) Step: After the E-step temporarily completes the dataset, the M-step updates the parameters. It recalculates the mean vector and covariance matrix from this newly filled-in dataset. These updated parameters (\(\mu^{(t+1)}\) and \(\Sigma^{(t+1)}\)) represent a better estimate of the data’s true distribution and will be used in the next E-step.

The two-step EM cycle repeats, with each iteration refining the imputed values and the parameter estimates. The algorithm stops once the changes to the imputed values between successive iterations become negligible (i.e., fall below a predefined tolerance threshold), indicating that a stable solution has been reached. The final output is the dataset with the missing entries replaced by their converged, model-based estimates.

1.14.5 Plotting the complete and imputed data

Smaller points colored in black represent observed data, and the circle-shapes colored in magenta denote the imputed data.

plot.me <- function(index1, index2){
  plot.imputed <- sim_data.imputed[row.names(
    subset(sim_data.df, is.na(sim_data.df[, index1]) | is.na(sim_data.df[, index2]))), ]
  p = ggplot(sim_data.imputed, aes_string( paste0("X",index1)  , paste0("X",index2 ))) + 
  geom_point(alpha = 0.5, size = 0.7)+theme_bw() +
  stat_ellipse(type = "norm", color = "#000099", alpha=0.5) +
  geom_point(data = plot.imputed, aes_string( paste0("X",index1) , paste0("X",(index2))),size = 1.5, color = "Magenta", alpha = 0.8)
} 

gridExtra::grid.arrange( plot.me(1,2), plot.me(5,6), plot.me(13,20), plot.me(18,19), nrow = 2)

index1=1; index2=5
plot.imputed <- sim_data.imputed[row.names(
    subset(sim_data.df, is.na(sim_data.df[, index1]) | is.na(sim_data.df[, index2]))), ]
  p = ggplot(sim_data.imputed, aes_string( paste0("X",index1)  , paste0("X",index2 ))) + 
  geom_point(alpha = 0.5, size = 0.7)+theme_bw() +
  stat_ellipse(type = "norm", color = "#000099", alpha=0.5) +
  geom_point(data = plot.imputed, aes_string( paste0("X",index1) , paste0("X",(index2))),size = 1.5, color = "Magenta", alpha = 0.8)
  
plot_ly(sim_data.imputed, x = ~X1, y = ~X5, type = "scatter", 
        mode = "markers") %>% 
    layout(title='Scatterplot: Improved Water Quality vs. Sanitation Facilities', 
           xaxis = list (title = 'Water Quality'), yaxis = list (title = 'Sanitation'))

1.14.6 Validation of EM-imputation using the R Package Amelia

# knitr::include_graphics("ammelia.png")

1.14.7 Comparison

Let’s use the amelia function to impute the original data sim_data_df and compare the results to the simpler manual EM_algorithm imputation defined above.

# install.packages("Amelia")
library(Amelia)
dim(sim_data.df)
## [1] 200  20
amelia.out <- amelia(sim_data.df, m = 5)
## -- Imputation 1 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
##  21
## 
## -- Imputation 2 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17
## 
## -- Imputation 3 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17
## 
## -- Imputation 4 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12 13 14 15
## 
## -- Imputation 5 --
## 
##   1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
##  21 22 23
amelia.out
## 
## Amelia output with 5 imputed datasets.
## Return code:  1 
## Message:  Normal EM convergence. 
## 
## Chain Lengths:
## --------------
## Imputation 1:  21
## Imputation 2:  17
## Imputation 3:  17
## Imputation 4:  15
## Imputation 5:  23
amelia.imputed.5 <- amelia.out$imputations[[5]]
  • Magenta-color and circle-shape denote manual imputation via EM_algorithm
  • Orange-color and square-shapes denote Amelia imputation
plot.ii2 <- function(index, index2){
  plot.imputed <- sim_data.imputed[row.names(
    subset(sim_data.df, is.na(sim_data.df[, index]) | is.na(sim_data.df[, index2]))), ]
  plot.imputed2 <- amelia.imputed.5[row.names(
  subset(sim_data.df, is.na(sim_data.df[, index]) | is.na(sim_data.df[, index2]))), ]
  p = ggplot(sim_data.imputed, aes_string( paste0("X",index)  , paste0("X",index2 ))) + 
  geom_point(alpha = 0.8, size = 0.7)+theme_bw() +
  stat_ellipse(type = "norm", color = "#000099", alpha=0.5) +
  geom_point(data = plot.imputed, aes_string( paste0("X",index) , paste0("X",(index2))),size = 2.5, color = "Magenta", alpha = 0.9, shape = 16) +  
    geom_point(data = plot.imputed2, aes( X1 , X2),size = 2.5, color = "#FF9933", alpha = 0.8, shape = 18)
  return(p)
} 

plot.ii2(2, 4)

plot.ii2(17, 18)

1.14.8 Density plots

Finally, we can compare the densities of the original, manually-imputed and Amelia-imputed datasets. Remember that in this simulation, we had about \(500\) observations missing out of the \(4,000\) that we synthetically generated.

# plot.ii3 <- function(index){
#   imputed <- sim_data.imputed[is.na(sim_data.df[, index]) , index]
#   imputed.amelia <- amelia.imputed.5[is.na(sim_data.df[, index]) , index]
#   observed <- sim_data.df[!is.na(sim_data.df[, index]) , index]
#   imputed.df <- data.frame(x = c(observed,imputed,imputed.amelia), category = c( rep("obs",length(observed)),rep("simpleImplement",length(imputed)) ,rep("amelia",length(imputed.amelia))   ) )
#   p = ggplot(imputed.df, aes(x=x, y =..density..)) + 
#     geom_density(aes(fill = category),alpha=0.3)+
#     theme_bw()
#   return(p)
#   }
# grid.arrange( plot.ii3(1),plot.ii3(2),plot.ii3(3),plot.ii3(4),plot.ii3(5),
#               plot.ii3(6),plot.ii3(7),plot.ii3(8),plot.ii3(9),plot.ii3(10),
#               nrow = 5)

library(tidyr)

myPlotly <- function(index){
  imputed <- sim_data.imputed[is.na(sim_data.df[, index]) , index]
  imputed.amelia <- amelia.imputed.5[is.na(sim_data.df[, index]) , index]
  observed <- sim_data.df[!is.na(sim_data.df[, index]) , index]
  imputed.df <- data.frame(x = c(observed,imputed,imputed.amelia), 
              category = c( rep("obs",length(observed)),rep("simpleImplement",length(imputed)),
                            rep("amelia",length(imputed.amelia))   ) )
  df_long <- as.data.frame(cbind(index=c(1:length(imputed.df$x)), 
                                 category=imputed.df$category, x=imputed.df$x))
  df_wide <- spread(df_long, category, x)
  
  p = plot_ly() %>%
    add_lines(x = ~density(as.numeric(df_wide$simpleImplement), na.rm = T)$x, 
        y= ~density(as.numeric(df_wide$simpleImplement), na.rm = T)$y, name = "EM", mode = 'lines') %>%
    add_lines(x = density(as.numeric(df_wide$amelia), na.rm = T)$x, 
        y= density(as.numeric(df_wide$amelia), na.rm = T)$y, name = "Amelia", mode = 'lines') %>%
    add_lines(x = ~density(as.numeric(df_wide$obs), na.rm = T)$x, 
        y= ~density(as.numeric(df_wide$obs), na.rm = T)$y, name = "Observed", mode = 'lines') %>%
    layout(title=sprintf("Distributions: Feature X.%d", index),
           xaxis = list(title = 'Measurements'),
           yaxis = list(title ="Densities"),
           legend = list(title="Distributions", orientation = 'h'))
  return(p)
}

# Plot a few features
myPlotly(5)
myPlotly(9)
myPlotly(10)
# grid.arrange( myPlotly(1),myPlotly(2),myPlotly(3),myPlotly(4),myPlotly(5),
#               myPlotly(6),myPlotly(7),myPlotly(8),myPlotly(9),myPlotly(10),
#               nrow = 5)

1.15 Parsing web pages and visualizing tabular HTML data

In this section, we will utilize the Earthquakes dataset on SOCR website. It records information about earthquakes that happened between 1969 and 2007 with magnitudes larger than 5 on the Richter scale. Here is how we parse the data on the source webpage and ingest the information into R:

# install.packages("xml2")
library("XML"); library("xml2")
library("rvest")
wiki_url <- read_html("https://wiki.socr.umich.edu/index.php/SOCR_Data_Dinov_021708_Earthquakes")
html_nodes(wiki_url, "#content")
## {xml_nodeset (1)}
## [1] <div id="content" class="mw-body" role="main">\n\t\t\t<a id="top"></a>\n\ ...
earthquake<- html_table(html_nodes(wiki_url, "table")[[2]])

In this dataset, Magt(magnitude type) may be used as a grouping variable. We will draw a “Longitude vs Latitude” line plot from this dataset. The function we are using is called ggplot() under ggplot2. The input type for this function is mostly data frame. aes() specifies axes.

# library(ggplot2)
# plot4<-ggplot(earthquake, aes(Longitude, Latitude, group=Magt, color=Magt))+ 
#   geom_point(data=earthquake, size=4, mapping=aes(x=Longitude, y=Latitude, shape=Magt))
# plot4  # or plint(plot4)
# https://plotly-r.com/working-with-symbols.html 

glyphication <- function (name) {
  glyph= vector()
  for (i in 1:length(name)){
    glyph[i]="triangle-up"
    if (name[i]=="Md") { glyph[i]="diamond-open" }
    else if (name[i]=="ML") { glyph[i]="circle-open" }
    else if (name[i]=="Mw") { glyph[i]="square-open" }
    else if (name[i]=="Mx") { glyph[i]="x-open" }
  }
  return(glyph)
}
earthquake$glyph <- glyphication(earthquake$Magt)
  
plot_ly(earthquake) %>%
  add_markers(x = ~Longitude, y = ~Latitude, type = "scatter",  color = ~Magt,
        mode = "markers", marker = list(size = ~Depth, color = ~Magt, symbol = ~glyph,
                                        line = list(color = "black",width = 2))) %>% 
    layout(title="California Earthquakes (1969 - 2007)")

The most important line of code has 2 parts. The first part ggplot(earthquake, aes(Longitude, Latitude, group=Magt, color=Magt)) specifies the setting of the plot: dataset, group and color. The second part specifies we are going to draw lines between data points. In later chapters we will frequently use the package ggplot2 and the structure under this great package is always function1+function2.

We can visualize the distribution for different variables using density plots. The following script plots the distribution for Latitude among different Magnitude types, also using the ggplot() function combined with geom_density().

plot5 <- ggplot(earthquake, aes(Latitude, size=1))+geom_density(aes(color=Magt))
plot5

We can also compute and display 2D Kernel Density and 3D Surface Plots. Plotting 2D Kernel Density and 3D Surface plots is very important and useful in multivariate exploratory data analytic.

We will use the plot_ly() function under the plotly package, which takes value from a data frame.

To create a surface plot, we use two vectors: x and y with length m and n respectively. We also need a matrix: z of size \(m\times n\). This z matrix is created from matrix multiplication between x and y.

The kde2d() function is needed for 2D kernel density estimation.

kernal_density <- with(earthquake, MASS::kde2d(Longitude, Latitude, n = 50))

Here z is an estimate of the kernel density function. Then we apply plot_ly to the list kernal_density via the with() function.

library(plotly)
with(kernal_density, plot_ly(x=x, y=y, z=z, type="surface"))

Note that we used the option "surface", however you can experiment with the type option.

Alternatively, one can plot 1D, 2D or 3D plots:

plot_ly(x = ~ earthquake$Longitude)
plot_ly(x = ~ earthquake$Longitude, y = ~earthquake$Latitude)
plot_ly(x = ~ earthquake$Longitude, y = ~earthquake$Latitude, z=~earthquake$Mag)
df3D <- data.frame(x=earthquake$Longitude, y=earthquake$Latitude, z=earthquake$Mag)

# Convert he Long (X, Y, Z) Earthquake format data into a Matrix Format

#  install.packages("Matrix")
library("Matrix")
matrix_EarthQuakes <- with(df3D, sparseMatrix(i = as.numeric(180-x), j=as.numeric(y), x=z, use.last.ij=T, dimnames=list(levels(x), levels(y))))
dim(matrix_EarthQuakes)
## [1] 307  44
# colnames(matrix_EarthQuakes) <- seq(from=earthquake$Longitude[1],
#                                     to=earthquake$Longitude[length(earthquake$Longitude)],
#                                                              length.out=dim(matrix_EarthQuakes)[2])
# rownames(matrix_EarthQuakes) <- seq(from=earthquake$Latitude[1],
#                                     to=earthquake$Latitude[length(earthquake$Latitude)],
#                                                              length.out=dim(matrix_EarthQuakes)[1])
# View(as.matrix(matrix_EarthQuakes))

# view matrix is 2D heatmap: 
library("ggplot2"); library("gplots")

# heatmap.2( as.matrix(matrix_EarthQuakes[280:307, 30:44]), Rowv=FALSE, Colv=FALSE, dendrogram='none', cellnote=as.matrix(matrix_EarthQuakes[280:307, 30:44]), notecol="black", trace='none', key=FALSE, lwid = c(.01, .99), lhei = c(.01, .99), margins = c(5, 15 ))

plot_ly(z = ~as.matrix(matrix_EarthQuakes[280:307, 30:44]), type = "heatmap") %>% hide_colorbar()
# plot_ly(x=~colnames(matrix_EarthQuakes[280:307, 30:44]), 
#         y=~rownames(matrix_EarthQuakes[280:307, 30:44]),
#         z = ~as.matrix(matrix_EarthQuakes[280:307, 30:44]), type = "heatmap") %>%
#   layout(title="California Earthquakes Heatmap", 
#          xaxis=list(title="Longitude"), yaxis=list(title="Latitude")) %>% 
#   hide_colorbar()

# Long -180<x<-170, Lat: 30<y<45, Z: 5<Mag<8
matrix_EarthQuakes <- with(df3D, sparseMatrix(i = as.numeric(180+x), j=as.numeric(y), x=z, use.last.ij=TRUE, dimnames=list(levels(x), levels(y))))
mat1 <- as.matrix(matrix_EarthQuakes)
plot_ly(z = ~mat1, type = "surface")
# To plot the Aggregate (Summed) Magnitudes at all Long/Lat:
matrix_EarthQuakes <- with(df3D, sparseMatrix(i = as.numeric(180+x), j=as.numeric(y), x=z, dimnames=list(levels(x), levels(y))))
mat1 <- as.matrix(matrix_EarthQuakes)
plot_ly(z = ~mat1, type = "surface")
# plot_ly(z = ~mat1[30:60, 20:40], type = "surface")

1.16 Cohort-Rebalancing (for Imbalanced Groups)

Comparing cohorts with imbalanced sample sizes (unbalanced designs) may present hidden biases in the results. Frequently, a cohort-rebalancing protocol is necessary to avoid such unexpected effects. Extremely unequal sample sizes can invalidate various parametric assumptions (e.g., homogeneity of variances). Also, there may be insufficient data representing the patterns belonging to the minority class(es) leading to inadequate capturing of the feature distributions. Although the groups do not have to have equal sizes, a general rule of thumb is that group sizes where one group is more than an order of magnitude larger than the size of another group has the potential for bias.

1.16.1 Example 1: Parkinson’s Diseases Study

This Parkinson’s diseases case-study involves neuroimaging, genetics, clinical, and phenotypic data for over 600 volunteers produced multivariate data for 3 cohorts – HC=Healthy Controls(166) , PD=Parkinson’s (434), SWEDD= subjects without evidence for dopaminergic deficit (61).

# update packages
# update.packages()

# load the data: 06_PPMI_ClassificationValidationData_Short.csv
ppmi_data <-read.csv("https://umich.instructure.com/files/330400/download?download_frd=1", header=TRUE)

table(ppmi_data$ResearchGroup)

# binarize the Dx classes
ppmi_data$ResearchGroup <- ifelse(ppmi_data$ResearchGroup == "Control", "Control", "Patient")
attach(ppmi_data)

head(ppmi_data)

# Model-free analysis, classification
# install.packages("crossval")
# install.packages("ada")
# library("crossval")
library(crossval)
library(ada)
#set up adaboosting prediction function


# Define a new classification result-reporting function
my.ada <- function (train.x, train.y, test.x, test.y, negative, formula){
  ada.fit <- ada(train.x, train.y)
  predict.y <- predict(ada.fit, test.x)
  #count TP, FP, TN, FN, Accuracy, etc.
  out <- confusionMatrix(test.y, predict.y, negative = negative)
 # negative  is the label of a negative "null" sample (default: "control").
  return (out)
}

# balance cases
# SMOTE: Synthetic Minority Oversampling Technique to handle class imbalance in binary classification.
set.seed(1000)
# https://cran.r-project.org/src/contrib/Archive/unbalanced/
# install.packages('mlr', 'FNN', 'RANN', 'unbalanced') to deal with unbalanced group data
library(unbalanced)
ppmi_data$PD <- ifelse(ppmi_data$ResearchGroup=="Control", 1, 0) 
uniqueID <- unique(ppmi_data$FID_IID) 
ppmi_data <- ppmi_data[ppmi_data$VisitID==1, ]
ppmi_data$PD <- factor(ppmi_data$PD)

colnames(ppmi_data)
# ppmi_data.1<-ppmi_data[, c(3:281, 284, 287, 336:340, 341)]
n <- ncol(ppmi_data)
output.1 <- ppmi_data$PD

# remove Default Real Clinical subject classifications! 
ppmi_data$PD <- ifelse(ppmi_data$ResearchGroup=="Control", 1, 0) 
input <- ppmi_data[ , -which(names(ppmi_data) %in% c("ResearchGroup", "PD", "X", "FID_IID"))]
# output <- as.matrix(ppmi_data[ , which(names(ppmi_data) %in% {"PD"})])
output <- as.factor(ppmi_data$PD)
c(dim(input), length(output))

#balance the dataset
data.1<-ubBalance(X= input, Y=output, type="ubSMOTE", percOver=300, percUnder=150, verbose=TRUE)
# percOver = A number that drives the decision of how many extra cases from the minority class are generated (known as over-sampling).
# k = A number indicating the number of nearest neighbors that are used to generate the new examples of the minority class.
# percUnder = A number that drives the decision of how many extra cases from the majority classes are selected for each case generated from the minority class (known as under-sampling)  

balancedData<-cbind(data.1$X, data.1$Y)
table(data.1$Y)

nrow(data.1$X); ncol(data.1$X)
nrow(balancedData); ncol(balancedData)
nrow(input); ncol(input)

colnames(balancedData) <- c(colnames(input), "PD")

# check visually for differences between the distributions of the raw (input) and rebalanced data (for only one variable, in this case)
QQ <- qqplot(input[, 5], balancedData [, 5], plot.it=F)

plot_ly(x=~QQ$x, y = ~QQ$y, type="scatter", mode="markers", showlegend=F) %>%
  add_lines(x=c(0,0.8), y=c(0,0.8), showlegend=F) %>%
  layout(title="QQ-Plot Original vs. Rebalanced Data", xaxis=list(title="original data"),
         yaxis=list(title="Rebalanced data"))

###Check balance
## Wilcoxon test
alpha.0.05 <- 0.05
test.results.bin <- NULL        # binarized/dichotomized p-values
test.results.raw <- NULL        # raw p-values

for (i in 1:(ncol(balancedData)-1)) {
    test.results.raw [i]  <- wilcox.test(input[, i], balancedData [, i])$p.value
    test.results.bin [i] <- ifelse(test.results.raw [i] > alpha.0.05, 1, 0)
    print(c("i=", i, "Wilcoxon-test=", test.results.raw [i]))
}
print(c("Wilcoxon test results: ", test.results.bin))

test.results.corr <- stats::p.adjust(test.results.raw, method = "fdr", n = length(test.results.raw)) 
# where methods are "holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none")
# plot(test.results.raw, test.results.corr)
# zeros (0) are significant independent between-group T-test differences, ones (1) are insignificant

plot_ly(x=~test.results.raw, y = ~test.results.corr, type="scatter", mode="markers", showlegend=F) %>%
  add_lines(x=c(0,1), y=c(0,1), showlegend=F) %>%
  layout(title="Wilcoxon test results - Original vs. Rebalanced Data", xaxis=list(title="Original"),
         yaxis=list(title="Rebalanced"))

# Check the Differences between the rate of significance between the raw and FDR-corrected p-values
test.results.bin <- ifelse(test.results.raw > alpha.0.05, 1, 0)
table(test.results.bin)
test.results.corr.bin <- ifelse(test.results.corr > alpha.0.05, 1, 0)
table(test.results.corr.bin)

Notes

  • SMOTE oversampling of the minority cohort is via generation of synthetic minority samples within the neighborhoods of observed observations. Thus, new minority instances blend observations in the same class and create clusters around each observed minority element.
  • The percOver parameter (perc.over/100) represents the number of new instances generated for each rare instance in the minority sample, when \(perc.over < 100\), a single instance is generated. For example, percOver=300 and percOver=30 would triple (300/100) and leave unchanged (30/100) the size of the minority sample, respectively.
  • The \(k\) parameter represents the number of neighbors to consider as the aggregate pool that the new examples are generated.
  • The percUnder (perc.under/100) represents the number of “normal” (majority class) instances that are randomly selected for each smoted (synthetically generated) observation. For instance, percUnder=300 or percUnder=30 would downsample the majority sample by choosing one-out-of-each-three or all of the majority sample points, respectively.

Continue to Visualization Chapter Part 2 includes exploratory data analytics (EDA), probability distributions, and mixture distribution modeling.

SOCR Resource Visitor number Web Analytics SOCR Email