| SOCR ≫ | DSPA ≫ | Topics ≫ |
Data Science and Predictive Analytics (UMich HS650)
Managing Data in R
SOCR/MIDAS (Ivo Dinov)
September 2017
In this chapter, 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 Saving and Loading R Data Structures
Let’s start by extracting the 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 centimetres). 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")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.
water <- read.csv('https://umich.instructure.com/files/399172/download?download_frd=1', header=T)
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 0colnames(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 0which.max(water$year); ## [1] 913# rowMeans(water[,5:6])
mean(water[,6], trim=0.08, na.rm=T)## [1] 71.63629This code loads CSV files that already include a header line listing the names of the variables. If we don’t have a header in the dataset, we can use the header = FALSE option to fix it. R will assign default names to the column variables of the dataset.
Simulation <- read.csv("https://umich.instructure.com/files/354289/download?download_frd=1", header = FALSE)
Simulation[1:3, ]## V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12
## 1 ID i2 age treat homeless pcs mcs cesd indtot pss_fr drugrisk sexrisk
## 2 1 0 25 0 0 49 7 46 37 0 1 6
## 3 2 18 31 0 0 48 34 17 48 0 0 11
## V13 V14 V15 V16
## 1 satreat female substance racegrp
## 2 0 0 cocaine black
## 3 0 0 alcohol whiteTo save a data frame to CSV files, we could use the write.csv() function. The option file = "a/local/file/path" allow us edit the saved file path.
write.csv(iris, file = "C:/Users/iris.csv")3 Exploring the Structure of Data
We can use the command str() to explore the structure of a dataset.
str(water)## 'data.frame': 3331 obs. of 6 variables:
## $ year : int 1990 1990 1990 1990 1990 1990 1990 1990 1990 1990 ...
## $ region : Factor w/ 6 levels "Africa","Americas",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ country : Factor w/ 192 levels "Afghanistan",..: 3 5 19 23 26 27 30 32 33 37 ...
## $ residence_area : Factor w/ 3 levels "Rural","Total",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ improved_water : num 88 42 49 86 39 67 34 46 37 83 ...
## $ sanitation_facilities: num 77 7 0 22 2 42 27 12 4 11 ...We can see that this World Drinking Water dataset has 3331 observations and 6 variables. The output also give us the class of each variable and first few elements in the variable.
4 Exploring Numeric Variables
Summary statistics for numeric variables in the dataset could be accessed by using the command summary().
summary(water$year)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1990 1995 2005 2002 2010 2012summary(water[c("improved_water", "sanitation_facilities")])## improved_water sanitation_facilities
## Min. : 3.0 Min. : 0.00
## 1st Qu.: 77.0 1st Qu.: 42.00
## Median : 93.0 Median : 81.00
## 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 :135plot(density(water$improved_water,na.rm = T)) # no need to be continuous, we can still get intuition about the variable distribution
The six summary statistics and NA’s (missing data) are in the output.
5 Measuring the 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] 65mean(c(40, 56, 99))## [1] 65median(vec1)## [1] 56median(c(40, 56, 99))## [1] 56# install.packages("psych");
library("psych")
geometric.mean(vec1, na.rm=TRUE)## [1] 60.52866The 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 variable is an example of a bimodal. We also have multimodal that has two or more modes in the data.
Mode is one of the measures for the central tendency. The best way to use it is to comparing the counts of the mode to other values. This help us to judge whether one or several categories dominates all others in the data. After that, we are able to analyze the story behind these common ones.
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.
6 Measuring Spread - 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 quantile/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 quantile/Q3 (
3rd Qu.), representing the \(75^{th}\) percentile, which splits off the lowest 75% of data from the top 25% - Maximum (
Max.), representing hte 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 give us both the minimum and maximum. An combination of range() and diff() could do the trick of getting the actual range value.
range(water$year)## [1] 1990 2012diff(range(water$year))## [1] 22Q1 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$year)## [1] 15summary(water$improved_water)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 3.0 77.0 93.0 84.9 99.0 100.0 32IQR(water$improved_water, na.rm = T)## [1] 22Just 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 = T)## 0% 25% 50% 75% 100%
## 3 77 93 99 100We 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 = T)## 20% 60%
## 71 97When 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 = T)## 0% 20% 40% 60% 80% 100%
## 3 71 89 97 100 100Let’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 interquatile 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. We have a very small minimum that makes the range of first quantile very large. This extreme value impacts the mean less than the median.
7 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")
In the boxplot we have five horizontal lines each represents 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.
8 Visualizing Numeric Variables - histograms
Histogram is another way to show the spread of a numeric variable (See Chapter 3 for additional information). It uses predetermined number of bins as containers for values to divide the original data. The height of the bins indicates frequency.
hist(water$improved_water, main = "Histogram of Percent improved_water", xlab="Percentage")
hist(water$sanitation_facilities, main = "Histogram of Percent sanitation_facilities", xlab = "Percentage")
We could see that the shape of two graphs are somewhat similar. They are both left skewed patterns (\(mean \lt median\)). Other common skew patterns are shown in the following picture.
These plots are generated by R and the code is provided in the appendix.
You can see the density plots of over 80 different probability distributions using the SOCR Java Distribution Calculators or the Distributome HTML5 Distribution Calculators.
9 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.
x <- 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)
Often, but not always, real world processes leand to normally distributed data. A normal distribution would have a higher frequency for middle values and lower frequency for more extreme values. It has a symetric 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.
10 Measuring Spread - variance and standard deviation
Distribution is a great way to characterize data using only a few parameters. For example, normal distribution can be defined by only two parameters center and spread or statistically mean and standard deviation.
The way to get mean value is to divide the sum of the data values by the number of values. So, we have the following formula.
\[Mean(X)=\mu=\frac{1}{n}\sum^{n}_{i=1} x_i\]
The standard deviation is the square root of the variance. Variance is the average sum of square.
\[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 non-normal, we use a new dataset about the demographics of baseball players to illustrate normal distribution properties. The "01_data.txt" in our class file has following variables:
- Name
- Team
- Position
- Height
- Weight
- Age
We check the histogram for approximate normality first.
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")
This plot allows us to visually inspect the normality of the players height and weight. We could also obtain mean and standard deviation of the weight and height variables.
mean(baseball$Weight)## [1] 201.7166mean(baseball$Height)## [1] 73.69729var(baseball$Weight)## [1] 440.9913sd(baseball$Weight)## [1] 20.99979var(baseball$Height)## [1] 5.316798sd(baseball$Height)## [1] 2.305818Larger standard deviation, or variance, suggest the data is more spread out from the mean. Therefore, the weight variable is more spread than the height variable.
Given the firt 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.
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.
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 have 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)
table(water$Year)##
## 1990 1995 2000 2005 2010 2012
## 520 561 570 570 556 554table(water$WHO.region)##
## Africa Americas Eastern Mediterranean
## 797 613 373
## Europe South-East Asia Western Pacific
## 910 191 447table(water$Residence.Area)##
## Rural Total Urban
## 1095 1109 1127Given 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..string.)
prop.table(year_table)##
## 1990 1995 2000 2005 2010 2012
## 0.1561093 0.1684179 0.1711198 0.1711198 0.1669168 0.1663164year_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.612 Exploring Relationships Between Variables
So far the methods and statistics that we have go through are at univariate level. Sometimes we want to examine the relationship between two or multiple variables. For example, did 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.
Visualizing Relationships - scatterplots
Let’s look at 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 dot. If the graph shows a clear pattern rather a group of messy dots or a horizontal line, the two variables may correlated with each other.
In R we can use 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$Population.using.improved,
main= "Scatterplot of Year vs. Improved_water",
xlab= "Year",
ylab= "Percent of Population Using Improved Water")
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, not of the regions had less than 20% of people using improved water sources while there used to be some regions that have such low percentages in the early years.
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 African WHO region first.
water$africa<-water$WHO.region=="Africa"Let’s revisit the table() function to see how many WHO regions are in Africa.
table(water$africa)##
## FALSE TRUE
## 2534 797Now, 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 give us the count that falls into its corresponding category. The Chi-square contribution provide 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.
13 Missing Data
In the previous sections, we simply ignored the incomplete observations in our water dataset (na.rm = TRUE). Is this an apporpriate strategy to handle incmplete 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 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$Population.using.improved.drinking, na.rm = T)
m2<-mean(water$Population.using.improved.sanitation, na.rm = T)
water_imp<-water
for(i in 1:3331){
if(is.na(water_imp$Population.using.improved.drinking[i])){
water_imp$Population.using.improved.drinking[i]=m1
}
if(is.na(water_imp$Population.using.improved.sanitation[i])){
water_imp$Population.using.improved.sanitation=m2
}
}
summary(water_imp)## Year..string. WHO.region..string. Country..string.
## 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.Type..string.
## Rural:1095
## Total:1109
## Urban:1127
##
##
##
##
## Population.using.improved.drinking.water.sources......numeric.
## Min. : 3.0
## 1st Qu.: 77.0
## Median : 93.0
## Mean : 84.9
## 3rd Qu.: 99.0
## Max. :100.0
## NA's :32
## Population.using.improved.sanitation.facilities......numeric.
## Min. : 0.00
## 1st Qu.: 42.00
## Median : 81.00
## Mean : 68.87
## 3rd Qu.: 97.00
## Max. :100.00
## NA's :135
## africa Population.using.improved.sanitation
## Mode :logical Min. :68.87
## FALSE:2534 1st Qu.:68.87
## TRUE :797 Median :68.87
## Mean :68.87
## 3rd Qu.:68.87
## Max. :68.87
##
## Population.using.improved.drinking
## Min. : 3.0
## 1st Qu.: 77.0
## Median : 93.0
## Mean : 84.9
## 3rd Qu.: 99.0
## Max. :100.0
## A more sophisticated way of resolving missing data is to use a model (e.g., linear regression) to predict teh 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", repos = "http://cran.us.r-project.org")
library(mi)## Loading required package: Matrix## Loading required package: stats4## mi (Version 1.0, packaged: 2015-04-16 14:03:10 UTC; goodrich)## mi Copyright (C) 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015 Trustees of Columbia University## This program comes with ABSOLUTELY NO WARRANTY.## This is free software, and you are welcome to redistribute it## under the General Public License version 2 or later.## Execute RShowDoc('COPYING') for details.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..string. WHO.region..string. Country..string.
## 1 1990 Africa Algeria
## 2 1990 Africa Angola
## 3 1990 Africa Benin
## 4 1990 Africa Botswana
## 5 1990 Africa Burkina Faso
## 6 1990 Africa Burundi
## Residence.Area.Type..string.
## 1 Rural
## 2 Rural
## 3 Rural
## 4 Rural
## 5 Rural
## 6 Rural
## Population.using.improved.drinking.water.sources......numeric.
## 1 88
## 2 42
## 3 49
## 4 86
## 5 39
## 6 67
## Population.using.improved.sanitation.facilities......numeric. africa
## 1 77 TRUE
## 2 7 TRUE
## 3 0 TRUE
## 4 22 TRUE
## 5 2 TRUE
## 6 42 TRUE
## missing_Population.using.improved.drinking.water.sources......numeric.
## 1 FALSE
## 2 FALSE
## 3 FALSE
## 4 FALSE
## 5 FALSE
## 6 FALSE
## missing_Population.using.improved.sanitation.facilities......numeric.
## 1 FALSE
## 2 FALSE
## 3 FALSE
## 4 FALSE
## 5 FALSE
## 6 FALSEshow(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
## Year..string. continuous
## WHO.region..string. unordered-categorical
## Country..string. unordered-categorical
## Residence.Area.Type..string. unordered-categorical
## Population.using.improved.drinking.water.sources......numeric. continuous
## Population.using.improved.sanitation.facilities......numeric. continuous
## africa binary
## missing
## Year..string. 0
## WHO.region..string. 0
## Country..string. 0
## Residence.Area.Type..string. 0
## Population.using.improved.drinking.water.sources......numeric. 32
## Population.using.improved.sanitation.facilities......numeric. 135
## africa 0
## method
## Year..string. <NA>
## WHO.region..string. <NA>
## Country..string. <NA>
## Residence.Area.Type..string. <NA>
## Population.using.improved.drinking.water.sources......numeric. ppd
## Population.using.improved.sanitation.facilities......numeric. ppd
## africa <NA>
## model
## Year..string. <NA>
## WHO.region..string. <NA>
## Country..string. <NA>
## Residence.Area.Type..string. <NA>
## Population.using.improved.drinking.water.sources......numeric. linear
## Population.using.improved.sanitation.facilities......numeric. linear
## africa <NA>
##
## family
## Year..string. <NA>
## WHO.region..string. <NA>
## Country..string. <NA>
## Residence.Area.Type..string. <NA>
## Population.using.improved.drinking.water.sources......numeric. gaussian
## Population.using.improved.sanitation.facilities......numeric. gaussian
## africa <NA>
## link
## Year..string. <NA>
## WHO.region..string. <NA>
## Country..string. <NA>
## Residence.Area.Type..string. <NA>
## Population.using.improved.drinking.water.sources......numeric. identity
## Population.using.improved.sanitation.facilities......numeric. identity
## africa <NA>
## transformation
## Year..string. standardize
## WHO.region..string. <NA>
## Country..string. <NA>
## Residence.Area.Type..string. <NA>
## Population.using.improved.drinking.water.sources......numeric. standardize
## Population.using.improved.sanitation.facilities......numeric. standardize
## africa <NA>mdf<-change(mdf, y="Population.using.improved.drinking", what = "imputation_method", to="pmm")
mdf<-change(mdf, y="Population.using.improved.sanitation", what = "imputation_method", to="pmm")- Notes:
- Converting the input
data.frameto amissing_data.frameallows 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. Themissing_data.frameconstructor suggests the most appropriate classes for each missing variable, however, the user often needs to correct, modify or change these meta-data, usingchange().- Use the
change()function to change/correct many meta-data in the constructedmissing_data.frameobject which are incorrect when usingshow(mfd). - To get a sense of the raw data, look at the
summary,image, orhistof the missing_data.frame. - The mi vignettes provide many useful examples of handling missing data.
We can perform the initial imputation. Here we imputed 3 times, which will create 3 different datasets 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.
data.frames <- complete(imputations, 3)
summary(water)## Year..string. WHO.region..string. Country..string.
## 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.Type..string.
## Rural:1095
## Total:1109
## Urban:1127
##
##
##
##
## Population.using.improved.drinking.water.sources......numeric.
## Min. : 3.0
## 1st Qu.: 77.0
## Median : 93.0
## Mean : 84.9
## 3rd Qu.: 99.0
## Max. :100.0
## NA's :32
## Population.using.improved.sanitation.facilities......numeric.
## Min. : 0.00
## 1st Qu.: 42.00
## Median : 81.00
## Mean : 68.87
## 3rd Qu.: 97.00
## Max. :100.00
## NA's :135
## africa
## Mode :logical
## FALSE:2534
## TRUE :797
##
##
##
## summary(data.frames[[1]])## Year..string. WHO.region..string. Country..string.
## 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.Type..string.
## Rural:1095
## Total:1109
## Urban:1127
##
##
##
##
## Population.using.improved.drinking.water.sources......numeric.
## Min. : 3.00
## 1st Qu.: 77.00
## Median : 93.00
## Mean : 84.82
## 3rd Qu.: 99.00
## Max. :100.00
##
## Population.using.improved.sanitation.facilities......numeric.
## Min. : 0.00
## 1st Qu.: 43.00
## Median : 81.00
## Mean : 69.07
## 3rd Qu.: 97.00
## Max. :100.00
##
## africa
## FALSE:2534
## TRUE : 797
##
##
##
##
##
## missing_Population.using.improved.drinking.water.sources......numeric.
## Mode :logical
## FALSE:3299
## TRUE :32
##
##
##
##
## missing_Population.using.improved.sanitation.facilities......numeric.
## Mode :logical
## FALSE:3196
## TRUE :135
##
##
##
## lapply(data.frames, summary)## $`chain:1`
## Year..string. WHO.region..string. Country..string.
## 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.Type..string.
## Rural:1095
## Total:1109
## Urban:1127
##
##
##
##
## Population.using.improved.drinking.water.sources......numeric.
## Min. : 3.00
## 1st Qu.: 77.00
## Median : 93.00
## Mean : 84.82
## 3rd Qu.: 99.00
## Max. :100.00
##
## Population.using.improved.sanitation.facilities......numeric.
## Min. : 0.00
## 1st Qu.: 43.00
## Median : 81.00
## Mean : 69.07
## 3rd Qu.: 97.00
## Max. :100.00
##
## africa
## FALSE:2534
## TRUE : 797
##
##
##
##
##
## missing_Population.using.improved.drinking.water.sources......numeric.
## Mode :logical
## FALSE:3299
## TRUE :32
##
##
##
##
## missing_Population.using.improved.sanitation.facilities......numeric.
## Mode :logical
## FALSE:3196
## TRUE :135
##
##
##
##
##
## $`chain:2`
## Year..string. WHO.region..string. Country..string.
## 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.Type..string.
## Rural:1095
## Total:1109
## Urban:1127
##
##
##
##
## Population.using.improved.drinking.water.sources......numeric.
## Min. : 3.00
## 1st Qu.: 77.00
## Median : 93.00
## Mean : 84.83
## 3rd Qu.: 99.00
## Max. :100.00
##
## Population.using.improved.sanitation.facilities......numeric.
## Min. : 0.00
## 1st Qu.: 43.00
## Median : 81.00
## Mean : 69.11
## 3rd Qu.: 97.00
## Max. :100.00
##
## africa
## FALSE:2534
## TRUE : 797
##
##
##
##
##
## missing_Population.using.improved.drinking.water.sources......numeric.
## Mode :logical
## FALSE:3299
## TRUE :32
##
##
##
##
## missing_Population.using.improved.sanitation.facilities......numeric.
## Mode :logical
## FALSE:3196
## TRUE :135
##
##
##
##
##
## $`chain:3`
## Year..string. WHO.region..string. Country..string.
## 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.Type..string.
## Rural:1095
## Total:1109
## Urban:1127
##
##
##
##
## Population.using.improved.drinking.water.sources......numeric.
## Min. : 3.00
## 1st Qu.: 77.00
## Median : 93.00
## Mean : 84.81
## 3rd Qu.: 99.00
## Max. :100.00
##
## Population.using.improved.sanitation.facilities......numeric.
## Min. : 0.00
## 1st Qu.: 43.00
## Median : 81.00
## Mean : 69.16
## 3rd Qu.: 97.00
## Max. :100.00
##
## africa
## FALSE:2534
## TRUE : 797
##
##
##
##
##
## missing_Population.using.improved.drinking.water.sources......numeric.
## Mode :logical
## FALSE:3299
## TRUE :32
##
##
##
##
## missing_Population.using.improved.sanitation.facilities......numeric.
## Mode :logical
## FALSE:3196
## 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.
13.1 Simulate some real multivariate data
Suppose we would like to generate a synthetic dataset: \[sim\_data=\{y, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10\}.\]
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 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 column has no missing 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 syntheticaly 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)## y x1 x2 x3 x4 x5 x6 x7 x8 x9
## 1 NA NA 0.000000 0 1 h 8 NA 3 NA
## 2 11.449223 NA 5.236938 0 1 i NA NA 10 0.2639489
## 3 -1.188296 0.000000 0.000000 0 5 a 3 -1.1469495 <NA> 0.4753195
## 4 NA NA NA 0 <NA> e 6 1.4810186 10 0.6696932
## 5 4.267916 3.490833 0.000000 0 <NA> <NA> NA 0.9161912 <NA> 0.9578455
## 6 NA 0.000000 4.384732 1 <NA> a NA NA 10 0.6095176
## x10
## 1 1
## 2 2
## 3 NA
## 4 3
## 5 8
## 6 6## y x1 x2 x3
## Min. :-3.846 Min. :0.000 Min. :0.000 Min. :0.0000
## 1st Qu.: 2.410 1st Qu.:0.000 1st Qu.:0.000 1st Qu.:0.0000
## Median : 5.646 Median :0.000 Median :3.068 Median :0.0000
## Mean : 5.560 Mean :2.473 Mean :2.545 Mean :0.4443
## 3rd Qu.: 8.503 3rd Qu.:4.958 3rd Qu.:4.969 3rd Qu.:1.0000
## Max. :16.487 Max. :8.390 Max. :8.421 Max. :1.0000
## NA's :300 NA's :300 NA's :300 NA's :300
## x4 x5 x6 x7 x8
## 1 :138 c : 80 Min. :1.00 Min. :-2.5689 3 : 78
## 2 :129 h : 76 1st Qu.:3.00 1st Qu.:-0.6099 7 : 77
## 3 :147 b : 74 Median :5.00 Median : 0.0202 5 : 75
## 4 :144 a : 73 Mean :4.93 Mean : 0.0435 4 : 73
## 5 :142 j : 72 3rd Qu.:7.00 3rd Qu.: 0.7519 1 : 70
## NA's:300 (Other):325 Max. :9.00 Max. : 3.7157 (Other):327
## NA's :300 NA's :300 NA's :300 NA's :300
## x9 x10
## Min. :0.1001 Min. : 0.000
## 1st Qu.:0.3206 1st Qu.: 2.000
## Median :0.5312 Median : 4.000
## Mean :0.5416 Mean : 3.929
## 3rd Qu.:0.7772 3rd Qu.: 5.000
## Max. :0.9895 Max. :11.000
## NA's :300 NA's :300# 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)## Object of class missing_data.frame with 1000 observations on 11 variables
##
## There are 542 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
## y continuous 300 ppd linear
## x1 continuous 300 ppd linear
## x2 continuous 300 ppd linear
## x3 binary 300 ppd logit
## x4 ordered-categorical 300 ppd ologit
## x5 unordered-categorical 300 ppd mlogit
## x6 continuous 300 ppd linear
## x7 continuous 300 ppd linear
## x8 unordered-categorical 300 ppd mlogit
## x9 proportion 300 ppd betareg
## x10 continuous 300 ppd linear
##
## family link transformation
## y gaussian identity standardize
## x1 gaussian identity standardize
## x2 gaussian identity standardize
## x3 binomial logit <NA>
## x4 multinomial logit <NA>
## x5 multinomial logit <NA>
## x6 gaussian identity standardize
## x7 gaussian identity standardize
## x8 multinomial logit <NA>
## x9 binomial logit identity
## x10 gaussian identity standardize# mdf@patterns # to get the textual missing pattern
image(mdf) # remember the visual pattern of this MCAR
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))
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 summary stats for the original data (prior to introducing missing
# values) with missing data and the re-completed data following imputation
summary(sim_data); summary(sim_data_30pct_missing); summary(data.frames[[1]]);## y x1 x2 x3 x4
## Min. :-3.846 Min. :0.000 Min. :0.000 Min. :0.000 1:200
## 1st Qu.: 2.489 1st Qu.:0.000 1st Qu.:0.000 1st Qu.:0.000 2:200
## Median : 5.549 Median :0.000 Median :2.687 Median :0.000 3:200
## Mean : 5.562 Mean :2.472 Mean :2.516 Mean :0.431 4:200
## 3rd Qu.: 8.325 3rd Qu.:4.996 3rd Qu.:5.007 3rd Qu.:1.000 5:200
## Max. :16.487 Max. :8.390 Max. :8.421 Max. :1.000
##
## x5 x6 x7 x8
## a :100 Min. :1.00 Min. :-3.26215 1 :100
## b :100 1st Qu.:3.00 1st Qu.:-0.62074 2 :100
## c :100 Median :5.00 Median : 0.02879 3 :100
## d :100 Mean :4.99 Mean : 0.03448 4 :100
## e :100 3rd Qu.:7.00 3rd Qu.: 0.72012 5 :100
## f :100 Max. :9.00 Max. : 3.71572 6 :100
## (Other):400 (Other):400
## x9 x10
## Min. :0.1001 Min. : 0.000
## 1st Qu.:0.3157 1st Qu.: 3.000
## Median :0.5263 Median : 4.000
## Mean :0.5362 Mean : 3.989
## 3rd Qu.:0.7603 3rd Qu.: 5.000
## Max. :0.9895 Max. :11.000
## ## y x1 x2 x3
## Min. :-3.846 Min. :0.000 Min. :0.000 Min. :0.0000
## 1st Qu.: 2.410 1st Qu.:0.000 1st Qu.:0.000 1st Qu.:0.0000
## Median : 5.646 Median :0.000 Median :3.068 Median :0.0000
## Mean : 5.560 Mean :2.473 Mean :2.545 Mean :0.4443
## 3rd Qu.: 8.503 3rd Qu.:4.958 3rd Qu.:4.969 3rd Qu.:1.0000
## Max. :16.487 Max. :8.390 Max. :8.421 Max. :1.0000
## NA's :300 NA's :300 NA's :300 NA's :300
## x4 x5 x6 x7 x8
## 1 :138 c : 80 Min. :1.00 Min. :-2.5689 3 : 78
## 2 :129 h : 76 1st Qu.:3.00 1st Qu.:-0.6099 7 : 77
## 3 :147 b : 74 Median :5.00 Median : 0.0202 5 : 75
## 4 :144 a : 73 Mean :4.93 Mean : 0.0435 4 : 73
## 5 :142 j : 72 3rd Qu.:7.00 3rd Qu.: 0.7519 1 : 70
## NA's:300 (Other):325 Max. :9.00 Max. : 3.7157 (Other):327
## NA's :300 NA's :300 NA's :300 NA's :300
## x9 x10
## Min. :0.1001 Min. : 0.000
## 1st Qu.:0.3206 1st Qu.: 2.000
## Median :0.5312 Median : 4.000
## Mean :0.5416 Mean : 3.929
## 3rd Qu.:0.7772 3rd Qu.: 5.000
## Max. :0.9895 Max. :11.000
## NA's :300 NA's :300## y x1 x2 x3 x4
## Min. :-3.846 Min. :-5.566 Min. :-4.686 0:543 1:210
## 1st Qu.: 2.409 1st Qu.: 0.000 1st Qu.: 0.000 1:457 2:196
## Median : 5.496 Median : 2.192 Median : 2.647 3:217
## Mean : 5.481 Mean : 2.464 Mean : 2.532 4:181
## 3rd Qu.: 8.462 3rd Qu.: 4.861 3rd Qu.: 4.815 5:196
## Max. :18.185 Max. :10.301 Max. :10.108
##
## x5 x6 x7 x8
## c :126 Min. :-6.676 Min. :-2.76931 1 :106
## h :105 1st Qu.: 3.000 1st Qu.:-0.60988 10 :106
## a :104 Median : 5.000 Median : 0.03448 3 :105
## d :104 Mean : 4.967 Mean : 0.06266 4 :105
## b :103 3rd Qu.: 7.000 3rd Qu.: 0.75200 7 :105
## j :101 Max. :14.832 Max. : 3.71572 5 :104
## (Other):357 (Other):369
## x9 x10 missing_y missing_x1
## Min. :0.005495 Min. :-1.289 Mode :logical Mode :logical
## 1st Qu.:0.332674 1st Qu.: 2.420 FALSE:700 FALSE:700
## Median :0.538199 Median : 4.000 TRUE :300 TRUE :300
## Mean :0.545506 Mean : 3.895
## 3rd Qu.:0.771012 3rd Qu.: 5.000
## Max. :0.989531 Max. :11.000
##
## missing_x2 missing_x3 missing_x4 missing_x5
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:700 FALSE:700 FALSE:700 FALSE:700
## TRUE :300 TRUE :300 TRUE :300 TRUE :300
##
##
##
##
## missing_x6 missing_x7 missing_x8 missing_x9
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:700 FALSE:700 FALSE:700 FALSE:700
## TRUE :300 TRUE :300 TRUE :300 TRUE :300
##
##
##
##
## missing_x10
## Mode :logical
## FALSE:700
## TRUE :300
##
##
##
## lapply(data.frames, summary)## $`chain:1`
## y x1 x2 x3 x4
## Min. :-3.846 Min. :-5.566 Min. :-4.686 0:543 1:210
## 1st Qu.: 2.409 1st Qu.: 0.000 1st Qu.: 0.000 1:457 2:196
## Median : 5.496 Median : 2.192 Median : 2.647 3:217
## Mean : 5.481 Mean : 2.464 Mean : 2.532 4:181
## 3rd Qu.: 8.462 3rd Qu.: 4.861 3rd Qu.: 4.815 5:196
## Max. :18.185 Max. :10.301 Max. :10.108
##
## x5 x6 x7 x8
## c :126 Min. :-6.676 Min. :-2.76931 1 :106
## h :105 1st Qu.: 3.000 1st Qu.:-0.60988 10 :106
## a :104 Median : 5.000 Median : 0.03448 3 :105
## d :104 Mean : 4.967 Mean : 0.06266 4 :105
## b :103 3rd Qu.: 7.000 3rd Qu.: 0.75200 7 :105
## j :101 Max. :14.832 Max. : 3.71572 5 :104
## (Other):357 (Other):369
## x9 x10 missing_y missing_x1
## Min. :0.005495 Min. :-1.289 Mode :logical Mode :logical
## 1st Qu.:0.332674 1st Qu.: 2.420 FALSE:700 FALSE:700
## Median :0.538199 Median : 4.000 TRUE :300 TRUE :300
## Mean :0.545506 Mean : 3.895
## 3rd Qu.:0.771012 3rd Qu.: 5.000
## Max. :0.989531 Max. :11.000
##
## missing_x2 missing_x3 missing_x4 missing_x5
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:700 FALSE:700 FALSE:700 FALSE:700
## TRUE :300 TRUE :300 TRUE :300 TRUE :300
##
##
##
##
## missing_x6 missing_x7 missing_x8 missing_x9
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:700 FALSE:700 FALSE:700 FALSE:700
## TRUE :300 TRUE :300 TRUE :300 TRUE :300
##
##
##
##
## missing_x10
## Mode :logical
## FALSE:700
## TRUE :300
##
##
##
##
##
## $`chain:2`
## y x1 x2 x3 x4
## Min. :-3.846 Min. :-5.161 Min. :-5.685 0:548 1:204
## 1st Qu.: 2.407 1st Qu.: 0.000 1st Qu.: 0.000 1:452 2:192
## Median : 5.527 Median : 2.383 Median : 2.062 3:206
## Mean : 5.457 Mean : 2.496 Mean : 2.419 4:192
## 3rd Qu.: 8.320 3rd Qu.: 4.891 3rd Qu.: 4.840 5:206
## Max. :16.487 Max. :11.005 Max. : 8.698
##
## x5 x6 x7 x8
## c :119 Min. :-2.956 Min. :-2.56887 1 :111
## b :110 1st Qu.: 3.000 1st Qu.:-0.61088 5 :111
## h :106 Median : 5.000 Median : 0.02879 7 :108
## a :100 Mean : 5.012 Mean : 0.04629 10 :108
## j : 99 3rd Qu.: 7.000 3rd Qu.: 0.76116 3 :103
## d : 98 Max. :11.715 Max. : 3.71572 9 :102
## (Other):368 (Other):357
## x9 x10 missing_y missing_x1
## Min. :0.01524 Min. :-2.241 Mode :logical Mode :logical
## 1st Qu.:0.32805 1st Qu.: 2.239 FALSE:700 FALSE:700
## Median :0.52917 Median : 4.000 TRUE :300 TRUE :300
## Mean :0.53887 Mean : 3.962
## 3rd Qu.:0.76346 3rd Qu.: 5.029
## Max. :0.99758 Max. :11.000
##
## missing_x2 missing_x3 missing_x4 missing_x5
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:700 FALSE:700 FALSE:700 FALSE:700
## TRUE :300 TRUE :300 TRUE :300 TRUE :300
##
##
##
##
## missing_x6 missing_x7 missing_x8 missing_x9
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:700 FALSE:700 FALSE:700 FALSE:700
## TRUE :300 TRUE :300 TRUE :300 TRUE :300
##
##
##
##
## missing_x10
## Mode :logical
## FALSE:700
## TRUE :300
##
##
##
##
##
## $`chain:3`
## y x1 x2 x3 x4
## Min. :-8.447 Min. :-6.476 Min. :-5.852 0:556 1:206
## 1st Qu.: 2.421 1st Qu.: 0.000 1st Qu.: 0.000 1:444 2:198
## Median : 5.502 Median : 2.427 Median : 2.294 3:192
## Mean : 5.448 Mean : 2.515 Mean : 2.425 4:196
## 3rd Qu.: 8.172 3rd Qu.: 4.891 3rd Qu.: 4.798 5:208
## Max. :16.487 Max. :10.438 Max. :10.236
##
## x5 x6 x7 x8
## h :113 Min. :-2.048 Min. :-2.61383 1 :111
## b :109 1st Qu.: 3.000 1st Qu.:-0.57509 10 :110
## c :108 Median : 5.000 Median : 0.07430 3 :109
## j :104 Mean : 4.897 Mean : 0.06363 7 :109
## a :103 3rd Qu.: 7.000 3rd Qu.: 0.75224 5 :108
## d :102 Max. :13.079 Max. : 3.71572 4 :102
## (Other):361 (Other):351
## x9 x10 missing_y missing_x1
## Min. :0.008103 Min. :-1.171 Mode :logical Mode :logical
## 1st Qu.:0.331581 1st Qu.: 2.466 FALSE:700 FALSE:700
## Median :0.543762 Median : 4.000 TRUE :300 TRUE :300
## Mean :0.547881 Mean : 3.933
## 3rd Qu.:0.772536 3rd Qu.: 5.000
## Max. :0.991607 Max. :11.000
##
## missing_x2 missing_x3 missing_x4 missing_x5
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:700 FALSE:700 FALSE:700 FALSE:700
## TRUE :300 TRUE :300 TRUE :300 TRUE :300
##
##
##
##
## missing_x6 missing_x7 missing_x8 missing_x9
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:700 FALSE:700 FALSE:700 FALSE:700
## TRUE :300 TRUE :300 TRUE :300 TRUE :300
##
##
##
##
## missing_x10
## Mode :logical
## FALSE:700
## TRUE :300
##
##
##
## 
Let’s check imputation convergence (details provided below).
round(mipply(imputations, mean, to.matrix = TRUE), 3)## chain:1 chain:2 chain:3
## y -0.010 -0.013 -0.014
## x1 -0.002 0.005 0.008
## x2 -0.002 -0.024 -0.023
## x3 1.457 1.452 1.444
## x4 2.957 3.004 3.002
## x5 5.397 5.390 5.452
## x6 0.007 0.016 -0.006
## x7 0.010 0.001 0.010
## x8 5.446 5.491 5.383
## x9 0.546 0.539 0.548
## x10 -0.009 0.009 0.001
## 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.300Rhats(imputations, statistic = "moments") # assess the convergence of MI algorithm## mean_y mean_x1 mean_x2 mean_x3 mean_x4 mean_x5 mean_x6
## 1.0275098 1.3719292 1.4187538 1.1120566 0.9671518 0.9984926 1.0673339
## mean_x7 mean_x8 mean_x9 mean_x10 sd_y sd_x1 sd_x2
## 1.0235087 1.1472099 1.0590998 1.3056564 1.2606094 0.9871043 1.6797089
## sd_x3 sd_x4 sd_x5 sd_x6 sd_x7 sd_x8 sd_x9
## 1.1093092 0.9032711 0.9416840 1.0689718 1.7036680 0.9389458 0.9210945
## sd_x10
## 1.1690558plot(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.74992 -0.39398 -0.04522 -0.04114 0.28420 1.57721
##
## $y$observed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.17512 -0.39354 0.01069 0.00000 0.36767 1.36503
##
##
## $x1
## $x1$is_missing
## missing
## FALSE TRUE
## 700 300
##
## $x1$imputed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.72547 -0.35089 -0.01430 0.01227 0.36775 1.64527
##
## $x1$observed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.4768 -0.4768 -0.4768 0.0000 0.4793 1.1411
##
##
## $x2
## $x2$is_missing
## missing
## FALSE TRUE
## 700 300
##
## $x2$imputed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.62939 -0.40979 -0.06545 -0.05554 0.31929 1.49239
##
## $x2$observed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.4938 -0.4938 0.1014 0.0000 0.4705 1.1402
##
##
## $x3
## $x3$crosstab
##
## observed imputed
## 0 1167 480
## 1 933 420
##
##
## $x4
## $x4$crosstab
##
## observed imputed
## 1 414 206
## 2 387 199
## 3 441 174
## 4 432 137
## 5 426 184
##
##
## $x5
## $x5$crosstab
##
## observed imputed
## a 219 88
## b 222 100
## c 240 113
## d 204 100
## e 189 75
## f 180 79
## g 198 82
## h 228 96
## i 204 79
## j 216 88
##
##
## $x6
## $x6$is_missing
## missing
## FALSE TRUE
## 700 300
##
## $x6$imputed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.25630 -0.36075 0.01636 0.01882 0.37789 1.92500
##
## $x6$observed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.76403 -0.37521 0.01361 0.00000 0.40243 0.79125
##
##
## $x7
## $x7$is_missing
## missing
## FALSE TRUE
## 700 300
##
## $x7$imputed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.43483 -0.31789 0.02400 0.02384 0.36757 1.66205
##
## $x7$observed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.33258 -0.33329 -0.01189 0.00000 0.36133 1.87322
##
##
## $x8
## $x8$crosstab
##
## observed imputed
## 1 210 118
## 2 204 88
## 3 234 83
## 4 219 81
## 5 225 98
## 6 165 82
## 7 231 91
## 8 192 65
## 9 210 80
## 10 210 114
##
##
## $x9
## $x9$is_missing
## missing
## FALSE TRUE
## 700 300
##
## $x9$imputed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.005495 0.354322 0.557638 0.549790 0.755216 0.997577
##
## $x9$observed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1001 0.3206 0.5312 0.5416 0.7772 0.9895
##
##
## $x10
## $x10$is_missing
## missing
## FALSE TRUE
## 700 300
##
## $x10$imputed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.598702 -0.364685 -0.005922 0.001275 0.347187 1.628924
##
## $x10$observed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.01801 -0.49975 0.01851 0.00000 0.27764 1.83243
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.36 0.31
## x1 0.98 0.02
## x2 0.93 0.03
## x31 -0.02 0.15
## x4.L 0.19 0.32
## x4.Q -0.11 0.14
## x4.C 0.08 0.28
## x4^4 0.14 0.24
## x5b 0.07 0.51
## x5c -0.39 0.32
## x5d -0.29 0.21
## x5e 0.28 0.24
## x5f -0.29 0.24
## x5g -0.44 0.44
## x5h 0.02 0.35
## x5i -0.61 0.79
## x5j 0.46 0.22
## x6 0.02 0.03
## x7 0.94 0.05
## x82 0.83 0.40
## x83 0.73 0.36
## x84 0.10 0.39
## x85 0.92 0.35
## x86 0.55 0.59
## x87 0.62 0.31
## x88 0.65 0.44
## x89 0.39 0.45
## x810 0.29 0.71
## x9 0.94 0.58
## x10 0.01 0.04
## n = 970, k = 30
## residual deviance = 2017.4, null deviance = 15762.3 (difference = 13744.9)
## overdispersion parameter = 2.1
## residual sd is sqrt(overdispersion) = 1.44##
## Call:
## pool(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 +
## x10, data = imputations, m = 3)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.2963 -0.7006 -0.0006 0.6916 3.3562
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.355417 0.306674 -1.159 0.2541
## x1 0.980063 0.019454 50.377 < 2e-16 ***
## x2 0.931090 0.025583 36.395 9.18e-10 ***
## x31 -0.019269 0.151571 -0.127 0.9034
## x4.L 0.193268 0.324358 0.596 0.6016
## x4.Q -0.108087 0.139686 -0.774 0.4583
## x4.C 0.080013 0.279050 0.287 0.7951
## x4^4 0.136430 0.236650 0.577 0.6044
## x5b 0.073158 0.512643 0.143 0.8962
## x5c -0.391046 0.320086 -1.222 0.2741
## x5d -0.294133 0.213543 -1.377 0.1689
## x5e 0.278255 0.238333 1.168 0.2477
## x5f -0.291974 0.242409 -1.204 0.2339
## x5g -0.437671 0.437821 -1.000 0.3840
## x5h 0.024308 0.345927 0.070 0.9469
## x5i -0.610682 0.793366 -0.770 0.5143
## x5j 0.459547 0.220178 2.087 0.0387 *
## x6 0.022907 0.033219 0.690 0.5289
## x7 0.938171 0.053132 17.657 < 2e-16 ***
## x82 0.828925 0.400756 2.068 0.1130
## x83 0.732812 0.364313 2.011 0.1125
## x84 0.096215 0.387074 0.249 0.8166
## x85 0.919834 0.354265 2.596 0.0542 .
## x86 0.552556 0.594507 0.929 0.4304
## x87 0.621744 0.305849 2.033 0.0870 .
## x88 0.651181 0.436424 1.492 0.2223
## x89 0.393601 0.454021 0.867 0.4465
## x810 0.289214 0.713890 0.405 0.7201
## x9 0.942941 0.582727 1.618 0.2256
## x10 0.007564 0.037570 0.201 0.8473
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 2.079832)
##
## Null deviance: 15762.3 on 999 degrees of freedom
## Residual deviance: 2017.4 on 970 degrees of freedom
## AIC: 3599.5
##
## Number of Fisher Scoring iterations: 7# Report the summaries of the imputations
data.frames <- complete(imputations, 3) # extract the first 3 chains
lapply(data.frames, summary)## $`chain:1`
## y x1 x2 x3 x4
## Min. :-3.846 Min. :-5.566 Min. :-4.686 0:543 1:210
## 1st Qu.: 2.409 1st Qu.: 0.000 1st Qu.: 0.000 1:457 2:196
## Median : 5.496 Median : 2.192 Median : 2.647 3:217
## Mean : 5.481 Mean : 2.464 Mean : 2.532 4:181
## 3rd Qu.: 8.462 3rd Qu.: 4.861 3rd Qu.: 4.815 5:196
## Max. :18.185 Max. :10.301 Max. :10.108
##
## x5 x6 x7 x8
## c :126 Min. :-6.676 Min. :-2.76931 1 :106
## h :105 1st Qu.: 3.000 1st Qu.:-0.60988 10 :106
## a :104 Median : 5.000 Median : 0.03448 3 :105
## d :104 Mean : 4.967 Mean : 0.06266 4 :105
## b :103 3rd Qu.: 7.000 3rd Qu.: 0.75200 7 :105
## j :101 Max. :14.832 Max. : 3.71572 5 :104
## (Other):357 (Other):369
## x9 x10 missing_y missing_x1
## Min. :0.005495 Min. :-1.289 Mode :logical Mode :logical
## 1st Qu.:0.332674 1st Qu.: 2.420 FALSE:700 FALSE:700
## Median :0.538199 Median : 4.000 TRUE :300 TRUE :300
## Mean :0.545506 Mean : 3.895
## 3rd Qu.:0.771012 3rd Qu.: 5.000
## Max. :0.989531 Max. :11.000
##
## missing_x2 missing_x3 missing_x4 missing_x5
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:700 FALSE:700 FALSE:700 FALSE:700
## TRUE :300 TRUE :300 TRUE :300 TRUE :300
##
##
##
##
## missing_x6 missing_x7 missing_x8 missing_x9
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:700 FALSE:700 FALSE:700 FALSE:700
## TRUE :300 TRUE :300 TRUE :300 TRUE :300
##
##
##
##
## missing_x10
## Mode :logical
## FALSE:700
## TRUE :300
##
##
##
##
##
## $`chain:2`
## y x1 x2 x3 x4
## Min. :-3.846 Min. :-5.161 Min. :-5.685 0:548 1:204
## 1st Qu.: 2.407 1st Qu.: 0.000 1st Qu.: 0.000 1:452 2:192
## Median : 5.527 Median : 2.383 Median : 2.062 3:206
## Mean : 5.457 Mean : 2.496 Mean : 2.419 4:192
## 3rd Qu.: 8.320 3rd Qu.: 4.891 3rd Qu.: 4.840 5:206
## Max. :16.487 Max. :11.005 Max. : 8.698
##
## x5 x6 x7 x8
## c :119 Min. :-2.956 Min. :-2.56887 1 :111
## b :110 1st Qu.: 3.000 1st Qu.:-0.61088 5 :111
## h :106 Median : 5.000 Median : 0.02879 7 :108
## a :100 Mean : 5.012 Mean : 0.04629 10 :108
## j : 99 3rd Qu.: 7.000 3rd Qu.: 0.76116 3 :103
## d : 98 Max. :11.715 Max. : 3.71572 9 :102
## (Other):368 (Other):357
## x9 x10 missing_y missing_x1
## Min. :0.01524 Min. :-2.241 Mode :logical Mode :logical
## 1st Qu.:0.32805 1st Qu.: 2.239 FALSE:700 FALSE:700
## Median :0.52917 Median : 4.000 TRUE :300 TRUE :300
## Mean :0.53887 Mean : 3.962
## 3rd Qu.:0.76346 3rd Qu.: 5.029
## Max. :0.99758 Max. :11.000
##
## missing_x2 missing_x3 missing_x4 missing_x5
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:700 FALSE:700 FALSE:700 FALSE:700
## TRUE :300 TRUE :300 TRUE :300 TRUE :300
##
##
##
##
## missing_x6 missing_x7 missing_x8 missing_x9
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:700 FALSE:700 FALSE:700 FALSE:700
## TRUE :300 TRUE :300 TRUE :300 TRUE :300
##
##
##
##
## missing_x10
## Mode :logical
## FALSE:700
## TRUE :300
##
##
##
##
##
## $`chain:3`
## y x1 x2 x3 x4
## Min. :-8.447 Min. :-6.476 Min. :-5.852 0:556 1:206
## 1st Qu.: 2.421 1st Qu.: 0.000 1st Qu.: 0.000 1:444 2:198
## Median : 5.502 Median : 2.427 Median : 2.294 3:192
## Mean : 5.448 Mean : 2.515 Mean : 2.425 4:196
## 3rd Qu.: 8.172 3rd Qu.: 4.891 3rd Qu.: 4.798 5:208
## Max. :16.487 Max. :10.438 Max. :10.236
##
## x5 x6 x7 x8
## h :113 Min. :-2.048 Min. :-2.61383 1 :111
## b :109 1st Qu.: 3.000 1st Qu.:-0.57509 10 :110
## c :108 Median : 5.000 Median : 0.07430 3 :109
## j :104 Mean : 4.897 Mean : 0.06363 7 :109
## a :103 3rd Qu.: 7.000 3rd Qu.: 0.75224 5 :108
## d :102 Max. :13.079 Max. : 3.71572 4 :102
## (Other):361 (Other):351
## x9 x10 missing_y missing_x1
## Min. :0.008103 Min. :-1.171 Mode :logical Mode :logical
## 1st Qu.:0.331581 1st Qu.: 2.466 FALSE:700 FALSE:700
## Median :0.543762 Median : 4.000 TRUE :300 TRUE :300
## Mean :0.547881 Mean : 3.933
## 3rd Qu.:0.772536 3rd Qu.: 5.000
## Max. :0.991607 Max. :11.000
##
## missing_x2 missing_x3 missing_x4 missing_x5
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:700 FALSE:700 FALSE:700 FALSE:700
## TRUE :300 TRUE :300 TRUE :300 TRUE :300
##
##
##
##
## missing_x6 missing_x7 missing_x8 missing_x9
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:700 FALSE:700 FALSE:700 FALSE:700
## TRUE :300 TRUE :300 TRUE :300 TRUE :300
##
##
##
##
## missing_x10
## Mode :logical
## FALSE:700
## TRUE :300
##
##
##
## coef(summary(model_results))[, 1:2] # get the model coef's and their SE's ## Estimate Std. Error
## (Intercept) -0.355417206 0.30667422
## x1 0.980063166 0.01945442
## x2 0.931090263 0.02558294
## x31 -0.019269375 0.15157072
## x4.L 0.193268041 0.32435818
## x4.Q -0.108087285 0.13968572
## x4.C 0.080013093 0.27905012
## x4^4 0.136430228 0.23665005
## x5b 0.073158388 0.51264271
## x5c -0.391046074 0.32008602
## x5d -0.294132971 0.21354346
## x5e 0.278255178 0.23833320
## x5f -0.291974092 0.24240915
## x5g -0.437671268 0.43782102
## x5h 0.024307790 0.34592738
## x5i -0.610681553 0.79336593
## x5j 0.459547432 0.22017803
## x6 0.022906845 0.03321945
## x7 0.938171458 0.05313197
## x82 0.828925387 0.40075639
## x83 0.732812041 0.36431326
## x84 0.096214639 0.38707372
## x85 0.919834159 0.35426525
## x86 0.552556161 0.59450709
## x87 0.621743568 0.30584934
## x88 0.651180921 0.43642430
## x89 0.393600695 0.45402091
## x810 0.289214174 0.71389033
## x9 0.942941415 0.58272734
## x10 0.007563922 0.03756981library("lattice")
densityplot(y ~ x1 + x2, data=imputations)
# To compare the density of observed data and imputed data --
# these should be similar (though not identical) under MAR assumption13.2 TBI Data Example
Next, we will see an example using the traumatic brain injury (TBI) dataset.
# 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"))
summary(TBI_Data)## id age sex mechanism
## Min. : 1.00 Min. :16.00 Female: 9 Bike_vs_Auto: 4
## 1st Qu.:12.25 1st Qu.:23.00 Male :37 Blunt : 4
## Median :23.50 Median :33.00 Fall :13
## Mean :23.50 Mean :36.89 GSW : 2
## 3rd Qu.:34.75 3rd Qu.:47.25 MCA : 7
## Max. :46.00 Max. :83.00 MVA :10
## Peds_vs_Auto: 6
## field.gcs er.gcs icu.gcs worst.gcs
## Min. : 3 Min. : 3.000 Min. : 0.000 Min. : 0.0
## 1st Qu.: 3 1st Qu.: 4.000 1st Qu.: 3.000 1st Qu.: 3.0
## Median : 7 Median : 7.500 Median : 6.000 Median : 3.0
## Mean : 8 Mean : 8.182 Mean : 6.378 Mean : 5.4
## 3rd Qu.:12 3rd Qu.:12.250 3rd Qu.: 8.000 3rd Qu.: 7.0
## Max. :15 Max. :15.000 Max. :14.000 Max. :14.0
## NA's :2 NA's :2 NA's :1 NA's :1
## X6m.gose X2013.gose skull.fx temp.injury
## Min. :2.000 Min. :2.000 Min. :0.0000 Min. :0.000
## 1st Qu.:3.000 1st Qu.:5.000 1st Qu.:0.0000 1st Qu.:0.000
## Median :5.000 Median :7.000 Median :1.0000 Median :1.000
## Mean :4.805 Mean :5.804 Mean :0.6087 Mean :0.587
## 3rd Qu.:6.000 3rd Qu.:7.000 3rd Qu.:1.0000 3rd Qu.:1.000
## Max. :8.000 Max. :8.000 Max. :1.0000 Max. :1.000
## NA's :5
## surgery spikes.hr min.hr max.hr
## Min. :0.0000 Min. : 1.280 Min. : 0.000 Min. : 12.00
## 1st Qu.:0.0000 1st Qu.: 5.357 1st Qu.: 0.000 1st Qu.: 35.25
## Median :1.0000 Median : 18.170 Median : 0.000 Median : 97.50
## Mean :0.6304 Mean : 52.872 Mean : 3.571 Mean : 241.89
## 3rd Qu.:1.0000 3rd Qu.: 57.227 3rd Qu.: 0.000 3rd Qu.: 312.75
## Max. :1.0000 Max. :294.000 Max. :42.000 Max. :1199.00
## NA's :18 NA's :18 NA's :18
## acute.sz late.sz ever.sz
## Min. :0.0000 Min. :0.0000 Min. :0.000
## 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.000
## Median :0.0000 Median :1.0000 Median :1.000
## Mean :0.1739 Mean :0.5652 Mean :0.587
## 3rd Qu.:0.0000 3rd Qu.:1.0000 3rd Qu.:1.000
## Max. :1.0000 Max. :1.0000 Max. :1.000
## # Get information matrix of the data
# (1) 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"show(mdf); mdf@patterns; image(mdf)## Object of class missing_data.frame with 46 observations on 19 variables
##
## There are 7 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
## id irrelevant 0 <NA> <NA>
## age continuous 0 <NA> <NA>
## sex binary 0 <NA> <NA>
## mechanism unordered-categorical 0 <NA> <NA>
## field.gcs continuous 2 ppd linear
## er.gcs continuous 2 ppd linear
## icu.gcs continuous 1 ppd linear
## worst.gcs continuous 1 ppd linear
## X6m.gose continuous 5 ppd linear
## X2013.gose continuous 0 <NA> <NA>
## skull.fx binary 0 <NA> <NA>
## temp.injury binary 0 <NA> <NA>
## surgery binary 0 <NA> <NA>
## spikes.hr continuous 18 ppd linear
## min.hr continuous 18 ppd linear
## max.hr continuous 18 ppd linear
## acute.sz binary 0 <NA> <NA>
## late.sz binary 0 <NA> <NA>
## ever.sz binary 0 <NA> <NA>
##
## family link transformation
## id <NA> <NA> <NA>
## age <NA> <NA> standardize
## sex <NA> <NA> <NA>
## mechanism <NA> <NA> <NA>
## field.gcs gaussian identity standardize
## er.gcs gaussian identity standardize
## icu.gcs gaussian identity standardize
## worst.gcs gaussian identity standardize
## X6m.gose gaussian identity standardize
## X2013.gose <NA> <NA> standardize
## skull.fx <NA> <NA> <NA>
## temp.injury <NA> <NA> <NA>
## surgery <NA> <NA> <NA>
## spikes.hr gaussian identity standardize
## min.hr gaussian identity standardize
## max.hr gaussian identity standardize
## acute.sz <NA> <NA> <NA>
## late.sz <NA> <NA> <NA>
## ever.sz <NA> <NA> <NA>## [1] spikes.hr, min.hr, max.hr
## [2] field.gcs
## [3] nothing
## [4] nothing
## [5] nothing
## [6] nothing
## [7] spikes.hr, min.hr, max.hr
## [8] nothing
## [9] nothing
## [10] nothing
## [11] nothing
## [12] nothing
## [13] spikes.hr, min.hr, max.hr
## [14] nothing
## [15] spikes.hr, min.hr, max.hr
## [16] spikes.hr, min.hr, max.hr
## [17] nothing
## [18] spikes.hr, min.hr, max.hr
## [19] spikes.hr, min.hr, max.hr
## [20] spikes.hr, min.hr, max.hr
## [21] X6m.gose, spikes.hr, min.hr, max.hr
## [22] nothing
## [23] spikes.hr, min.hr, max.hr
## [24] spikes.hr, min.hr, max.hr
## [25] spikes.hr, min.hr, max.hr
## [26] spikes.hr, min.hr, max.hr
## [27] spikes.hr, min.hr, max.hr
## [28] X6m.gose
## [29] spikes.hr, min.hr, max.hr
## [30] nothing
## [31] X6m.gose, spikes.hr, min.hr, max.hr
## [32] spikes.hr, min.hr, max.hr
## [33] nothing
## [34] nothing
## [35] nothing
## [36] nothing
## [37] field.gcs, er.gcs, icu.gcs, worst.gcs, X6m.gose
## [38] er.gcs
## [39] nothing
## [40] nothing
## [41] nothing
## [42] spikes.hr, min.hr, max.hr
## [43] nothing
## [44] nothing
## [45] nothing
## [46] X6m.gose
## 7 Levels: nothing field.gcs X6m.gose er.gcs ... field.gcs, er.gcs, icu.gcs, worst.gcs, X6m.gose
# (2) 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")# (3) examine missingness patterns
summary(mdf); hist(mdf); ## id age sex mechanism
## Min. : 1.00 Min. :16.00 Female: 9 Bike_vs_Auto: 4
## 1st Qu.:12.25 1st Qu.:23.00 Male :37 Blunt : 4
## Median :23.50 Median :33.00 Fall :13
## Mean :23.50 Mean :36.89 GSW : 2
## 3rd Qu.:34.75 3rd Qu.:47.25 MCA : 7
## Max. :46.00 Max. :83.00 MVA :10
## Peds_vs_Auto: 6
## field.gcs er.gcs icu.gcs worst.gcs
## Min. : 3 Min. : 3.000 Min. : 0.000 Min. : 0.0
## 1st Qu.: 3 1st Qu.: 4.000 1st Qu.: 3.000 1st Qu.: 3.0
## Median : 7 Median : 7.500 Median : 6.000 Median : 3.0
## Mean : 8 Mean : 8.182 Mean : 6.378 Mean : 5.4
## 3rd Qu.:12 3rd Qu.:12.250 3rd Qu.: 8.000 3rd Qu.: 7.0
## Max. :15 Max. :15.000 Max. :14.000 Max. :14.0
## NA's :2 NA's :2 NA's :1 NA's :1
## X6m.gose X2013.gose skull.fx temp.injury
## Min. :2.000 Min. :2.000 Min. :0.0000 Min. :0.000
## 1st Qu.:3.000 1st Qu.:5.000 1st Qu.:0.0000 1st Qu.:0.000
## Median :5.000 Median :7.000 Median :1.0000 Median :1.000
## Mean :4.805 Mean :5.804 Mean :0.6087 Mean :0.587
## 3rd Qu.:6.000 3rd Qu.:7.000 3rd Qu.:1.0000 3rd Qu.:1.000
## Max. :8.000 Max. :8.000 Max. :1.0000 Max. :1.000
## NA's :5
## surgery spikes.hr min.hr max.hr
## Min. :0.0000 Min. : 1.280 Min. : 0.000 Min. : 12.00
## 1st Qu.:0.0000 1st Qu.: 5.357 1st Qu.: 0.000 1st Qu.: 35.25
## Median :1.0000 Median : 18.170 Median : 0.000 Median : 97.50
## Mean :0.6304 Mean : 52.872 Mean : 3.571 Mean : 241.89
## 3rd Qu.:1.0000 3rd Qu.: 57.227 3rd Qu.: 0.000 3rd Qu.: 312.75
## Max. :1.0000 Max. :294.000 Max. :42.000 Max. :1199.00
## NA's :18 NA's :18 NA's :18
## acute.sz late.sz ever.sz
## Min. :0.0000 Min. :0.0000 Min. :0.000
## 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.000
## Median :0.0000 Median :1.0000 Median :1.000
## Mean :0.1739 Mean :0.5652 Mean :0.587
## 3rd Qu.:0.0000 3rd Qu.:1.0000 3rd Qu.:1.000
## Max. :1.0000 Max. :1.0000 Max. :1.000
## image(mdf)
# (4) Perform initial imputation
imputations <- mi(mdf, n.iter=10, n.chains=5, verbose=TRUE)
hist(imputations)
# (5) Extracts several multiply imputed data.frames from "imputations" object
data.frames <- complete(imputations, 5)
# (6) Report a list of "summaries" for each element (imputation instance)
lapply(data.frames, summary)## $`chain:1`
## id age sex mechanism
## Min. : 1.00 Min. :16.00 Female: 9 Bike_vs_Auto: 4
## 1st Qu.:12.25 1st Qu.:23.00 Male :37 Blunt : 4
## Median :23.50 Median :33.00 Fall :13
## Mean :23.50 Mean :36.89 GSW : 2
## 3rd Qu.:34.75 3rd Qu.:47.25 MCA : 7
## Max. :46.00 Max. :83.00 MVA :10
## Peds_vs_Auto: 6
## field.gcs er.gcs icu.gcs worst.gcs
## Min. :-5.113 Min. :-1.043 Min. :-0.429 Min. : 0.000
## 1st Qu.: 3.000 1st Qu.: 4.000 1st Qu.: 3.000 1st Qu.: 3.000
## Median : 7.000 Median : 7.500 Median : 6.000 Median : 3.000
## Mean : 7.807 Mean : 8.167 Mean : 6.230 Mean : 5.329
## 3rd Qu.:12.000 3rd Qu.:12.750 3rd Qu.: 7.750 3rd Qu.: 7.000
## Max. :15.000 Max. :16.710 Max. :14.000 Max. :14.000
##
## X6m.gose X2013.gose skull.fx temp.injury surgery
## Min. :2.000 Min. :2.000 0:18 0:19 0:17
## 1st Qu.:3.000 1st Qu.:5.000 1:28 1:27 1:29
## Median :5.000 Median :7.000
## Mean :4.971 Mean :5.804
## 3rd Qu.:6.423 3rd Qu.:7.000
## Max. :8.639 Max. :8.000
##
## spikes.hr min.hr max.hr acute.sz late.sz
## Min. :-170.001 Min. :-19.18 Min. :-1320.74 0:38 0:20
## 1st Qu.: 2.225 1st Qu.: 0.00 1st Qu.: 16.12 1: 8 1:26
## Median : 18.170 Median : 0.00 Median : 76.00
## Mean : 41.218 Mean : 2.65 Mean : 101.13
## 3rd Qu.: 73.225 3rd Qu.: 3.75 3rd Qu.: 276.68
## Max. : 294.000 Max. : 42.00 Max. : 1199.00
##
## ever.sz missing_field.gcs missing_er.gcs missing_icu.gcs
## 0:19 Mode :logical Mode :logical Mode :logical
## 1:27 FALSE:44 FALSE:44 FALSE:45
## TRUE :2 TRUE :2 TRUE :1
##
##
##
##
## missing_worst.gcs missing_X6m.gose missing_spikes.hr missing_min.hr
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:45 FALSE:41 FALSE:28 FALSE:28
## TRUE :1 TRUE :5 TRUE :18 TRUE :18
##
##
##
##
## missing_max.hr
## Mode :logical
## FALSE:28
## TRUE :18
##
##
##
##
##
## $`chain:2`
## id age sex mechanism
## Min. : 1.00 Min. :16.00 Female: 9 Bike_vs_Auto: 4
## 1st Qu.:12.25 1st Qu.:23.00 Male :37 Blunt : 4
## Median :23.50 Median :33.00 Fall :13
## Mean :23.50 Mean :36.89 GSW : 2
## 3rd Qu.:34.75 3rd Qu.:47.25 MCA : 7
## Max. :46.00 Max. :83.00 MVA :10
## Peds_vs_Auto: 6
## field.gcs er.gcs icu.gcs worst.gcs
## Min. : 1.588 Min. :-4.273 Min. : 0.000 Min. :-3.496
## 1st Qu.: 3.000 1st Qu.: 4.000 1st Qu.: 3.000 1st Qu.: 3.000
## Median : 7.000 Median : 7.500 Median : 6.000 Median : 3.000
## Mean : 8.151 Mean : 8.005 Mean : 6.583 Mean : 5.207
## 3rd Qu.:12.000 3rd Qu.:12.376 3rd Qu.: 8.000 3rd Qu.: 7.000
## Max. :21.344 Max. :15.000 Max. :15.825 Max. :14.000
##
## X6m.gose X2013.gose skull.fx temp.injury surgery
## Min. :2.000 Min. :2.000 0:18 0:19 0:17
## 1st Qu.:3.000 1st Qu.:5.000 1:28 1:27 1:29
## Median :5.000 Median :7.000
## Mean :4.975 Mean :5.804
## 3rd Qu.:6.608 3rd Qu.:7.000
## Max. :8.000 Max. :8.000
##
## spikes.hr min.hr max.hr acute.sz late.sz
## Min. :-100.965 Min. :-21.251 Min. : -76.25 0:38 0:20
## 1st Qu.: 4.577 1st Qu.: 0.000 1st Qu.: 37.50 1: 8 1:26
## Median : 26.317 Median : 0.000 Median : 197.46
## Mean : 57.065 Mean : 3.915 Mean : 271.03
## 3rd Qu.: 75.442 3rd Qu.: 6.142 3rd Qu.: 383.24
## Max. : 303.904 Max. : 42.000 Max. :1199.00
##
## ever.sz missing_field.gcs missing_er.gcs missing_icu.gcs
## 0:19 Mode :logical Mode :logical Mode :logical
## 1:27 FALSE:44 FALSE:44 FALSE:45
## TRUE :2 TRUE :2 TRUE :1
##
##
##
##
## missing_worst.gcs missing_X6m.gose missing_spikes.hr missing_min.hr
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:45 FALSE:41 FALSE:28 FALSE:28
## TRUE :1 TRUE :5 TRUE :18 TRUE :18
##
##
##
##
## missing_max.hr
## Mode :logical
## FALSE:28
## TRUE :18
##
##
##
##
##
## $`chain:3`
## id age sex mechanism
## Min. : 1.00 Min. :16.00 Female: 9 Bike_vs_Auto: 4
## 1st Qu.:12.25 1st Qu.:23.00 Male :37 Blunt : 4
## Median :23.50 Median :33.00 Fall :13
## Mean :23.50 Mean :36.89 GSW : 2
## 3rd Qu.:34.75 3rd Qu.:47.25 MCA : 7
## Max. :46.00 Max. :83.00 MVA :10
## Peds_vs_Auto: 6
## field.gcs er.gcs icu.gcs worst.gcs
## Min. : 3.000 Min. :-12.584 Min. :-7.138 Min. :-2.715
## 1st Qu.: 3.250 1st Qu.: 3.250 1st Qu.: 3.000 1st Qu.: 3.000
## Median : 7.000 Median : 7.000 Median : 6.000 Median : 3.000
## Mean : 8.237 Mean : 7.601 Mean : 6.084 Mean : 5.224
## 3rd Qu.:12.000 3rd Qu.: 12.000 3rd Qu.: 7.750 3rd Qu.: 7.000
## Max. :19.527 Max. : 15.000 Max. :14.000 Max. :14.000
##
## X6m.gose X2013.gose skull.fx temp.injury surgery
## Min. :-0.2918 Min. :2.000 0:18 0:19 0:17
## 1st Qu.: 3.0000 1st Qu.:5.000 1:28 1:27 1:29
## Median : 5.0000 Median :7.000
## Mean : 4.4636 Mean :5.804
## 3rd Qu.: 6.0000 3rd Qu.:7.000
## Max. : 8.0000 Max. :8.000
##
## spikes.hr min.hr max.hr acute.sz late.sz
## Min. : 1.280 Min. :-13.939 Min. :-1083.22 0:38 0:20
## 1st Qu.: 7.209 1st Qu.: 0.000 1st Qu.: 28.75 1: 8 1:26
## Median : 46.665 Median : 0.000 Median : 128.41
## Mean : 83.302 Mean : 4.905 Mean : 253.89
## 3rd Qu.:129.118 3rd Qu.: 7.862 3rd Qu.: 432.71
## Max. :391.360 Max. : 42.000 Max. : 1548.28
##
## ever.sz missing_field.gcs missing_er.gcs missing_icu.gcs
## 0:19 Mode :logical Mode :logical Mode :logical
## 1:27 FALSE:44 FALSE:44 FALSE:45
## TRUE :2 TRUE :2 TRUE :1
##
##
##
##
## missing_worst.gcs missing_X6m.gose missing_spikes.hr missing_min.hr
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:45 FALSE:41 FALSE:28 FALSE:28
## TRUE :1 TRUE :5 TRUE :18 TRUE :18
##
##
##
##
## missing_max.hr
## Mode :logical
## FALSE:28
## TRUE :18
##
##
##
##
##
## $`chain:4`
## id age sex mechanism
## Min. : 1.00 Min. :16.00 Female: 9 Bike_vs_Auto: 4
## 1st Qu.:12.25 1st Qu.:23.00 Male :37 Blunt : 4
## Median :23.50 Median :33.00 Fall :13
## Mean :23.50 Mean :36.89 GSW : 2
## 3rd Qu.:34.75 3rd Qu.:47.25 MCA : 7
## Max. :46.00 Max. :83.00 MVA :10
## Peds_vs_Auto: 6
## field.gcs er.gcs icu.gcs worst.gcs
## Min. :-2.922 Min. :-10.221 Min. : 0.000 Min. : 0.000
## 1st Qu.: 3.000 1st Qu.: 4.000 1st Qu.: 3.000 1st Qu.: 3.000
## Median : 6.500 Median : 7.500 Median : 6.000 Median : 3.000
## Mean : 7.677 Mean : 7.994 Mean : 6.555 Mean : 5.463
## 3rd Qu.:12.000 3rd Qu.: 12.750 3rd Qu.: 8.000 3rd Qu.: 7.750
## Max. :15.000 Max. : 17.962 Max. :14.544 Max. :14.000
##
## X6m.gose X2013.gose skull.fx temp.injury surgery
## Min. :2.000 Min. :2.000 0:18 0:19 0:17
## 1st Qu.:3.000 1st Qu.:5.000 1:28 1:27 1:29
## Median :5.000 Median :7.000
## Mean :4.983 Mean :5.804
## 3rd Qu.:6.579 3rd Qu.:7.000
## Max. :8.914 Max. :8.000
##
## spikes.hr min.hr max.hr acute.sz late.sz
## Min. :-44.726 Min. :-20.896 Min. :-177.6 0:38 0:20
## 1st Qu.: 6.393 1st Qu.: 0.000 1st Qu.: 31.5 1: 8 1:26
## Median : 51.160 Median : 0.000 Median : 150.0
## Mean : 76.172 Mean : 2.828 Mean : 268.3
## 3rd Qu.:128.062 3rd Qu.: 3.970 3rd Qu.: 454.2
## Max. :294.000 Max. : 42.000 Max. :1199.0
##
## ever.sz missing_field.gcs missing_er.gcs missing_icu.gcs
## 0:19 Mode :logical Mode :logical Mode :logical
## 1:27 FALSE:44 FALSE:44 FALSE:45
## TRUE :2 TRUE :2 TRUE :1
##
##
##
##
## missing_worst.gcs missing_X6m.gose missing_spikes.hr missing_min.hr
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:45 FALSE:41 FALSE:28 FALSE:28
## TRUE :1 TRUE :5 TRUE :18 TRUE :18
##
##
##
##
## missing_max.hr
## Mode :logical
## FALSE:28
## TRUE :18
##
##
##
##
##
## $`chain:5`
## id age sex mechanism
## Min. : 1.00 Min. :16.00 Female: 9 Bike_vs_Auto: 4
## 1st Qu.:12.25 1st Qu.:23.00 Male :37 Blunt : 4
## Median :23.50 Median :33.00 Fall :13
## Mean :23.50 Mean :36.89 GSW : 2
## 3rd Qu.:34.75 3rd Qu.:47.25 MCA : 7
## Max. :46.00 Max. :83.00 MVA :10
## Peds_vs_Auto: 6
## field.gcs er.gcs icu.gcs worst.gcs
## Min. : 3.000 Min. : 3.000 Min. : 0.000 Min. : 0.000
## 1st Qu.: 3.250 1st Qu.: 4.250 1st Qu.: 3.000 1st Qu.: 3.000
## Median : 7.000 Median : 7.500 Median : 6.000 Median : 3.000
## Mean : 7.989 Mean : 8.434 Mean : 6.276 Mean : 5.361
## 3rd Qu.:12.000 3rd Qu.:12.750 3rd Qu.: 7.750 3rd Qu.: 7.000
## Max. :15.000 Max. :20.978 Max. :14.000 Max. :14.000
##
## X6m.gose X2013.gose skull.fx temp.injury surgery
## Min. :-0.3156 Min. :2.000 0:18 0:19 0:17
## 1st Qu.: 3.0000 1st Qu.:5.000 1:28 1:27 1:29
## Median : 5.0000 Median :7.000
## Mean : 4.5839 Mean :5.804
## 3rd Qu.: 6.0000 3rd Qu.:7.000
## Max. : 8.0000 Max. :8.000
##
## spikes.hr min.hr max.hr acute.sz late.sz
## Min. :-128.214 Min. :-16.903 Min. :-150.41 0:38 0:20
## 1st Qu.: 5.518 1st Qu.: 0.000 1st Qu.: 42.25 1: 8 1:26
## Median : 44.740 Median : 0.000 Median : 226.75
## Mean : 74.843 Mean : 3.293 Mean : 293.07
## 3rd Qu.: 124.832 3rd Qu.: 6.250 3rd Qu.: 420.23
## Max. : 451.717 Max. : 42.000 Max. :1302.91
##
## ever.sz missing_field.gcs missing_er.gcs missing_icu.gcs
## 0:19 Mode :logical Mode :logical Mode :logical
## 1:27 FALSE:44 FALSE:44 FALSE:45
## TRUE :2 TRUE :2 TRUE :1
##
##
##
##
## missing_worst.gcs missing_X6m.gose missing_spikes.hr missing_min.hr
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:45 FALSE:41 FALSE:28 FALSE:28
## TRUE :1 TRUE :5 TRUE :18 TRUE :18
##
##
##
##
## missing_max.hr
## Mode :logical
## FALSE:28
## TRUE :18
##
##
##
## # (6.a) To cast the imputed numbers as integers (not necessary, but may be useful)
indx <- sapply(data.frames[[5]], is.numeric) # get the indices of numeric columns
data.frames[[5]][indx] <- lapply(data.frames[[5]][indx], function(x) as.numeric(as.integer(x))) # cast each value as integer
# data.frames[[5]]$spikes.hr
# (7) Save results out
write.csv(data.frames[[5]], "C:\\Users\\Dinov\\Desktop\\TBI_MIData.csv")
# (8) Complete Data analytics functions:
# library("mi")
#lm.mi(); glm.mi(); polr.mi(); bayesglm.mi(); bayespolr.mi(); lmer.mi(); glmer.mi()
# (8.1) Define Linear Regression for multiply imputed dataset - Also see Step (10)
##linear regression for each imputed data set - 5 regression models are fit
fit_lm1 <- glm(ever.sz ~ surgery + worst.gcs + factor(sex) + age, data.frames$`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.frames$`chain:1`)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7497 -1.1362 0.7954 1.0197 1.4494
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.5931082 1.4204271 0.418 0.676
## surgery1 1.1165195 0.7162158 1.559 0.119
## worst.gcs -0.1182699 0.1052842 -1.123 0.261
## factor(sex)Male -0.4009911 0.8530059 -0.470 0.638
## age 0.0008816 0.0197598 0.045 0.964
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 62.371 on 45 degrees of freedom
## Residual deviance: 59.249 on 41 degrees of freedom
## AIC: 69.249
##
## Number of Fisher Scoring iterations: 4## glm(formula = ever.sz ~ surgery + worst.gcs + factor(sex) + age,
## family = "binomial", data = data.frames$`chain:1`)
## coef.est coef.se
## (Intercept) 0.59 1.42
## surgery1 1.12 0.72
## worst.gcs -0.12 0.11
## factor(sex)Male -0.40 0.85
## age 0.00 0.02
## ---
## n = 46, k = 5
## residual deviance = 59.2, null deviance = 62.4 (difference = 3.1)# 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=imputations, m=5)
display (model_results); summary (model_results) ## bayesglm(formula = ever.sz ~ surgery + worst.gcs + factor(sex) +
## age, data = imputations, m = 5, family = "binomial")
## coef.est coef.se
## (Intercept) 0.49 1.34
## surgery1 0.95 0.66
## worst.gcs -0.09 0.10
## factor(sex)Male -0.34 0.77
## age 0.00 0.02
## n = 41, k = 5
## residual deviance = 59.3, null deviance = 62.4 (difference = 3.1)##
## Call:
## pool(formula = ever.sz ~ surgery + worst.gcs + factor(sex) +
## age, data = imputations, m = 5, family = "binomial")
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6847 -1.1710 0.8342 1.0231 1.3889
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.492162 1.341430 0.367 0.714
## surgery1 0.953403 0.660858 1.443 0.149
## worst.gcs -0.094873 0.097244 -0.976 0.329
## factor(sex)Male -0.336636 0.769951 -0.437 0.662
## age 0.001334 0.019605 0.068 0.946
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 62.371 on 45 degrees of freedom
## Residual deviance: 59.263 on 41 degrees of freedom
## AIC: 69.263
##
## Number of Fisher Scoring iterations: 6.8# Report the summaries of the imputations
data.frames <- complete(imputations, 3) # extract the first 3 chains
lapply(data.frames, summary)## $`chain:1`
## id age sex mechanism
## Min. : 1.00 Min. :16.00 Female: 9 Bike_vs_Auto: 4
## 1st Qu.:12.25 1st Qu.:23.00 Male :37 Blunt : 4
## Median :23.50 Median :33.00 Fall :13
## Mean :23.50 Mean :36.89 GSW : 2
## 3rd Qu.:34.75 3rd Qu.:47.25 MCA : 7
## Max. :46.00 Max. :83.00 MVA :10
## Peds_vs_Auto: 6
## field.gcs er.gcs icu.gcs worst.gcs
## Min. :-5.113 Min. :-1.043 Min. :-0.429 Min. : 0.000
## 1st Qu.: 3.000 1st Qu.: 4.000 1st Qu.: 3.000 1st Qu.: 3.000
## Median : 7.000 Median : 7.500 Median : 6.000 Median : 3.000
## Mean : 7.807 Mean : 8.167 Mean : 6.230 Mean : 5.329
## 3rd Qu.:12.000 3rd Qu.:12.750 3rd Qu.: 7.750 3rd Qu.: 7.000
## Max. :15.000 Max. :16.710 Max. :14.000 Max. :14.000
##
## X6m.gose X2013.gose skull.fx temp.injury surgery
## Min. :2.000 Min. :2.000 0:18 0:19 0:17
## 1st Qu.:3.000 1st Qu.:5.000 1:28 1:27 1:29
## Median :5.000 Median :7.000
## Mean :4.971 Mean :5.804
## 3rd Qu.:6.423 3rd Qu.:7.000
## Max. :8.639 Max. :8.000
##
## spikes.hr min.hr max.hr acute.sz late.sz
## Min. :-170.001 Min. :-19.18 Min. :-1320.74 0:38 0:20
## 1st Qu.: 2.225 1st Qu.: 0.00 1st Qu.: 16.12 1: 8 1:26
## Median : 18.170 Median : 0.00 Median : 76.00
## Mean : 41.218 Mean : 2.65 Mean : 101.13
## 3rd Qu.: 73.225 3rd Qu.: 3.75 3rd Qu.: 276.68
## Max. : 294.000 Max. : 42.00 Max. : 1199.00
##
## ever.sz missing_field.gcs missing_er.gcs missing_icu.gcs
## 0:19 Mode :logical Mode :logical Mode :logical
## 1:27 FALSE:44 FALSE:44 FALSE:45
## TRUE :2 TRUE :2 TRUE :1
##
##
##
##
## missing_worst.gcs missing_X6m.gose missing_spikes.hr missing_min.hr
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:45 FALSE:41 FALSE:28 FALSE:28
## TRUE :1 TRUE :5 TRUE :18 TRUE :18
##
##
##
##
## missing_max.hr
## Mode :logical
## FALSE:28
## TRUE :18
##
##
##
##
##
## $`chain:2`
## id age sex mechanism
## Min. : 1.00 Min. :16.00 Female: 9 Bike_vs_Auto: 4
## 1st Qu.:12.25 1st Qu.:23.00 Male :37 Blunt : 4
## Median :23.50 Median :33.00 Fall :13
## Mean :23.50 Mean :36.89 GSW : 2
## 3rd Qu.:34.75 3rd Qu.:47.25 MCA : 7
## Max. :46.00 Max. :83.00 MVA :10
## Peds_vs_Auto: 6
## field.gcs er.gcs icu.gcs worst.gcs
## Min. : 1.588 Min. :-4.273 Min. : 0.000 Min. :-3.496
## 1st Qu.: 3.000 1st Qu.: 4.000 1st Qu.: 3.000 1st Qu.: 3.000
## Median : 7.000 Median : 7.500 Median : 6.000 Median : 3.000
## Mean : 8.151 Mean : 8.005 Mean : 6.583 Mean : 5.207
## 3rd Qu.:12.000 3rd Qu.:12.376 3rd Qu.: 8.000 3rd Qu.: 7.000
## Max. :21.344 Max. :15.000 Max. :15.825 Max. :14.000
##
## X6m.gose X2013.gose skull.fx temp.injury surgery
## Min. :2.000 Min. :2.000 0:18 0:19 0:17
## 1st Qu.:3.000 1st Qu.:5.000 1:28 1:27 1:29
## Median :5.000 Median :7.000
## Mean :4.975 Mean :5.804
## 3rd Qu.:6.608 3rd Qu.:7.000
## Max. :8.000 Max. :8.000
##
## spikes.hr min.hr max.hr acute.sz late.sz
## Min. :-100.965 Min. :-21.251 Min. : -76.25 0:38 0:20
## 1st Qu.: 4.577 1st Qu.: 0.000 1st Qu.: 37.50 1: 8 1:26
## Median : 26.317 Median : 0.000 Median : 197.46
## Mean : 57.065 Mean : 3.915 Mean : 271.03
## 3rd Qu.: 75.442 3rd Qu.: 6.142 3rd Qu.: 383.24
## Max. : 303.904 Max. : 42.000 Max. :1199.00
##
## ever.sz missing_field.gcs missing_er.gcs missing_icu.gcs
## 0:19 Mode :logical Mode :logical Mode :logical
## 1:27 FALSE:44 FALSE:44 FALSE:45
## TRUE :2 TRUE :2 TRUE :1
##
##
##
##
## missing_worst.gcs missing_X6m.gose missing_spikes.hr missing_min.hr
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:45 FALSE:41 FALSE:28 FALSE:28
## TRUE :1 TRUE :5 TRUE :18 TRUE :18
##
##
##
##
## missing_max.hr
## Mode :logical
## FALSE:28
## TRUE :18
##
##
##
##
##
## $`chain:3`
## id age sex mechanism
## Min. : 1.00 Min. :16.00 Female: 9 Bike_vs_Auto: 4
## 1st Qu.:12.25 1st Qu.:23.00 Male :37 Blunt : 4
## Median :23.50 Median :33.00 Fall :13
## Mean :23.50 Mean :36.89 GSW : 2
## 3rd Qu.:34.75 3rd Qu.:47.25 MCA : 7
## Max. :46.00 Max. :83.00 MVA :10
## Peds_vs_Auto: 6
## field.gcs er.gcs icu.gcs worst.gcs
## Min. : 3.000 Min. :-12.584 Min. :-7.138 Min. :-2.715
## 1st Qu.: 3.250 1st Qu.: 3.250 1st Qu.: 3.000 1st Qu.: 3.000
## Median : 7.000 Median : 7.000 Median : 6.000 Median : 3.000
## Mean : 8.237 Mean : 7.601 Mean : 6.084 Mean : 5.224
## 3rd Qu.:12.000 3rd Qu.: 12.000 3rd Qu.: 7.750 3rd Qu.: 7.000
## Max. :19.527 Max. : 15.000 Max. :14.000 Max. :14.000
##
## X6m.gose X2013.gose skull.fx temp.injury surgery
## Min. :-0.2918 Min. :2.000 0:18 0:19 0:17
## 1st Qu.: 3.0000 1st Qu.:5.000 1:28 1:27 1:29
## Median : 5.0000 Median :7.000
## Mean : 4.4636 Mean :5.804
## 3rd Qu.: 6.0000 3rd Qu.:7.000
## Max. : 8.0000 Max. :8.000
##
## spikes.hr min.hr max.hr acute.sz late.sz
## Min. : 1.280 Min. :-13.939 Min. :-1083.22 0:38 0:20
## 1st Qu.: 7.209 1st Qu.: 0.000 1st Qu.: 28.75 1: 8 1:26
## Median : 46.665 Median : 0.000 Median : 128.41
## Mean : 83.302 Mean : 4.905 Mean : 253.89
## 3rd Qu.:129.118 3rd Qu.: 7.862 3rd Qu.: 432.71
## Max. :391.360 Max. : 42.000 Max. : 1548.28
##
## ever.sz missing_field.gcs missing_er.gcs missing_icu.gcs
## 0:19 Mode :logical Mode :logical Mode :logical
## 1:27 FALSE:44 FALSE:44 FALSE:45
## TRUE :2 TRUE :2 TRUE :1
##
##
##
##
## missing_worst.gcs missing_X6m.gose missing_spikes.hr missing_min.hr
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:45 FALSE:41 FALSE:28 FALSE:28
## TRUE :1 TRUE :5 TRUE :18 TRUE :18
##
##
##
##
## missing_max.hr
## Mode :logical
## FALSE:28
## TRUE :18
##
##
##
## # (9) 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(imputations, 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.021 0.017 0.026 -0.036 -0.001
## er.gcs -0.002 -0.021 -0.070 -0.022 0.030
## icu.gcs -0.023 0.032 -0.045 0.027 -0.016
## worst.gcs -0.011 -0.029 -0.026 0.009 -0.006
## X6m.gose 0.048 0.049 -0.098 0.051 -0.063
## 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 41.218 57.065 83.302 76.172 74.843
## min.hr -0.047 0.018 0.068 -0.038 -0.014
## max.hr -0.232 0.048 0.020 0.043 0.084
## 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
# across chains.
# 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(imputations, statistic = "moments") # assess the convergence of MI algorithm## mean_field.gcs mean_er.gcs mean_icu.gcs mean_worst.gcs mean_X6m.gose
## 2.084123 2.023806 1.729935 1.545901 2.403296
## mean_spikes.hr mean_min.hr mean_max.hr sd_field.gcs sd_er.gcs
## 2.099167 1.531441 1.909896 1.401151 1.323421
## sd_icu.gcs sd_worst.gcs sd_X6m.gose sd_spikes.hr sd_min.hr
## 1.121314 1.231602 1.070467 1.726487 1.170868
## sd_max.hr
## 2.973138# When convergence is unstable, we can continue the iterations for all chains, e.g.
imputations <- mi(imputations, 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(imputations); hist(imputations); image(imputations); summary(imputations)
## $id
## $id$is_missing
## [1] "all values observed"
##
## $id$observed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.00 12.25 23.50 23.50 34.75 46.00
##
##
## $age
## $age$is_missing
## [1] "all values observed"
##
## $age$observed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.6045 -0.4019 -0.1126 0.0000 0.2997 1.3341
##
##
## $sex
## $sex$is_missing
## [1] "all values observed"
##
## $sex$observed
##
## 1 2
## 9 37
##
##
## $mechanism
## $mechanism$is_missing
## [1] "all values observed"
##
## $mechanism$observed
##
## 1 2 3 4 5 6 7
## 4 4 13 2 7 10 6
##
##
## $field.gcs
## $field.gcs$is_missing
## missing
## FALSE TRUE
## 44 2
##
## $field.gcs$imputed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.0804 -1.3838 -0.8009 -0.4953 -0.5242 2.8615
##
## $field.gcs$observed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.5545 -0.5545 -0.1109 0.0000 0.4436 0.7763
##
##
## $er.gcs
## $er.gcs$is_missing
## missing
## FALSE TRUE
## 44 2
##
## $er.gcs$imputed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.8863 -0.5104 0.7002 0.5768 1.8371 3.8846
##
## $er.gcs$observed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.62098 -0.50114 -0.08171 0.00000 0.48752 0.81708
##
##
## $icu.gcs
## $icu.gcs$is_missing
## missing
## FALSE TRUE
## 45 1
##
## $icu.gcs$imputed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.4552 -1.1662 -0.8078 -0.7051 -0.5337 0.4377
##
## $icu.gcs$observed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.98350 -0.52088 -0.05826 0.00000 0.25016 1.17540
##
##
## $worst.gcs
## $worst.gcs$is_missing
## missing
## FALSE TRUE
## 45 1
##
## $worst.gcs$imputed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.2106 -0.9268 -0.3532 -0.4798 -0.1576 0.2492
##
## $worst.gcs$observed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.8068 -0.3586 -0.3586 0.0000 0.2390 1.2849
##
##
## $X6m.gose
## $X6m.gose$is_missing
## missing
## FALSE TRUE
## 41 5
##
## $X6m.gose$imputed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.4935 -0.8170 -0.5343 -0.4788 -0.1878 1.4437
##
## $X6m.gose$observed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.80159 -0.51581 0.05576 0.00000 0.34155 0.91312
##
##
## $X2013.gose
## $X2013.gose$is_missing
## [1] "all values observed"
##
## $X2013.gose$observed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.9987 -0.2112 0.3139 0.0000 0.3139 0.5764
##
##
## $skull.fx
## $skull.fx$is_missing
## [1] "all values observed"
##
## $skull.fx$observed
##
## 1 2
## 18 28
##
##
## $temp.injury
## $temp.injury$is_missing
## [1] "all values observed"
##
## $temp.injury$observed
##
## 1 2
## 19 27
##
##
## $surgery
## $surgery$is_missing
## [1] "all values observed"
##
## $surgery$observed
##
## 1 2
## 17 29
##
##
## $spikes.hr
## $spikes.hr$is_missing
## missing
## FALSE TRUE
## 28 18
##
## $spikes.hr$imputed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -225.70 22.51 127.26 129.91 230.00 513.79
##
## $spikes.hr$observed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.280 5.357 18.170 52.872 57.227 294.000
##
##
## $min.hr
## $min.hr$is_missing
## missing
## FALSE TRUE
## 28 18
##
## $min.hr$imputed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.7672 -0.3709 0.1133 0.1950 0.8848 2.5339
##
## $min.hr$observed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.1828 -0.1828 -0.1828 0.0000 -0.1828 1.9668
##
##
## $max.hr
## $max.hr$is_missing
## missing
## FALSE TRUE
## 28 18
##
## $max.hr$imputed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.09247 0.06149 0.48296 0.46771 0.92447 2.18651
##
## $max.hr$observed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.3783 -0.3400 -0.2376 0.0000 0.1166 1.5750
##
##
## $acute.sz
## $acute.sz$is_missing
## [1] "all values observed"
##
## $acute.sz$observed
##
## 1 2
## 38 8
##
##
## $late.sz
## $late.sz$is_missing
## [1] "all values observed"
##
## $late.sz$observed
##
## 1 2
## 20 26
##
##
## $ever.sz
## $ever.sz$is_missing
## [1] "all values observed"
##
## $ever.sz$observed
##
## 1 2
## 19 27# (10) 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 = imputations, m = 5 ); display (model_results); summary (model_results)## bayesglm(formula = ever.sz ~ surgery + worst.gcs + factor(sex) +
## age, data = imputations, m = 5)
## coef.est coef.se
## (Intercept) 0.47 1.33
## surgery1 0.95 0.66
## worst.gcs -0.09 0.10
## factor(sex)Male -0.33 0.77
## age 0.00 0.02
## n = 41, k = 5
## residual deviance = 59.3, null deviance = 62.4 (difference = 3.0)##
## Call:
## pool(formula = ever.sz ~ surgery + worst.gcs + factor(sex) +
## age, data = imputations, m = 5)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.684 -1.173 0.835 1.020 1.387
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.473429 1.334510 0.355 0.723
## surgery1 0.946044 0.658712 1.436 0.151
## worst.gcs -0.093440 0.096987 -0.963 0.335
## factor(sex)Male -0.330296 0.768624 -0.430 0.667
## age 0.001662 0.019502 0.085 0.932
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 62.371 on 45 degrees of freedom
## Residual deviance: 59.331 on 41 degrees of freedom
## AIC: 69.331
##
## Number of Fisher Scoring iterations: 6.6coef(summary(model_results))[, 1:2] # get the model coef's and their SE's ## Estimate Std. Error
## (Intercept) 0.473428525 1.33450965
## surgery1 0.946044017 0.65871196
## worst.gcs -0.093439549 0.09698740
## factor(sex)Male -0.330295948 0.76862392
## age 0.001661906 0.01950154# Report the summaries of the imputations
data.frames <- complete(imputations, 3) # extract the first 3 chains
lapply(data.frames, summary) # report summaries## [[1]]
## id age sex mechanism
## Min. : 1.00 Min. :16.00 Female: 9 Bike_vs_Auto: 4
## 1st Qu.:12.25 1st Qu.:23.00 Male :37 Blunt : 4
## Median :23.50 Median :33.00 Fall :13
## Mean :23.50 Mean :36.89 GSW : 2
## 3rd Qu.:34.75 3rd Qu.:47.25 MCA : 7
## Max. :46.00 Max. :83.00 MVA :10
## Peds_vs_Auto: 6
## field.gcs er.gcs icu.gcs worst.gcs
## Min. : 3.000 Min. : 3.000 Min. :-3.059 Min. :-0.8032
## 1st Qu.: 3.116 1st Qu.: 4.250 1st Qu.: 3.000 1st Qu.: 3.0000
## Median : 7.000 Median : 8.000 Median : 6.000 Median : 3.0000
## Mean : 8.462 Mean : 8.607 Mean : 6.173 Mean : 5.2651
## 3rd Qu.:12.000 3rd Qu.:12.750 3rd Qu.: 7.750 3rd Qu.: 7.0000
## Max. :33.801 Max. :25.117 Max. :14.000 Max. :14.0000
##
## X6m.gose X2013.gose skull.fx temp.injury surgery
## Min. :-0.4212 Min. :2.000 0:18 0:19 0:17
## 1st Qu.: 3.0000 1st Qu.:5.000 1:28 1:27 1:29
## Median : 5.0000 Median :7.000
## Mean : 4.5295 Mean :5.804
## 3rd Qu.: 6.0000 3rd Qu.:7.000
## Max. : 8.0000 Max. :8.000
##
## spikes.hr min.hr max.hr acute.sz late.sz
## Min. :-76.990 Min. :-16.321 Min. :-295.8 0:38 0:20
## 1st Qu.: 5.518 1st Qu.: 0.000 1st Qu.: 47.0 1: 8 1:26
## Median : 45.290 Median : 0.000 Median : 210.3
## Mean :102.316 Mean : 6.737 Mean : 316.9
## 3rd Qu.:195.119 3rd Qu.: 14.636 3rd Qu.: 561.7
## Max. :429.659 Max. : 42.000 Max. :1199.0
##
## ever.sz missing_field.gcs missing_er.gcs missing_icu.gcs
## 0:19 Mode :logical Mode :logical Mode :logical
## 1:27 FALSE:44 FALSE:44 FALSE:45
## TRUE :2 TRUE :2 TRUE :1
##
##
##
##
## missing_worst.gcs missing_X6m.gose missing_spikes.hr missing_min.hr
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:45 FALSE:41 FALSE:28 FALSE:28
## TRUE :1 TRUE :5 TRUE :18 TRUE :18
##
##
##
##
## missing_max.hr
## Mode :logical
## FALSE:28
## TRUE :18
##
##
##
##
##
## [[2]]
## id age sex mechanism
## Min. : 1.00 Min. :16.00 Female: 9 Bike_vs_Auto: 4
## 1st Qu.:12.25 1st Qu.:23.00 Male :37 Blunt : 4
## Median :23.50 Median :33.00 Fall :13
## Mean :23.50 Mean :36.89 GSW : 2
## 3rd Qu.:34.75 3rd Qu.:47.25 MCA : 7
## Max. :46.00 Max. :83.00 MVA :10
## Peds_vs_Auto: 6
## field.gcs er.gcs icu.gcs worst.gcs
## Min. :-10.759 Min. :-15.903 Min. : 0.000 Min. : 0.000
## 1st Qu.: 3.000 1st Qu.: 4.000 1st Qu.: 3.000 1st Qu.: 3.000
## Median : 7.000 Median : 7.500 Median : 6.000 Median : 3.000
## Mean : 7.857 Mean : 8.363 Mean : 6.439 Mean : 5.377
## 3rd Qu.: 12.000 3rd Qu.: 12.750 3rd Qu.: 8.000 3rd Qu.: 7.000
## Max. : 20.189 Max. : 40.597 Max. :14.000 Max. :14.000
##
## X6m.gose X2013.gose skull.fx temp.injury surgery
## Min. :2.000 Min. :2.000 0:18 0:19 0:17
## 1st Qu.:3.000 1st Qu.:5.000 1:28 1:27 1:29
## Median :5.000 Median :7.000
## Mean :4.663 Mean :5.804
## 3rd Qu.:6.000 3rd Qu.:7.000
## Max. :8.000 Max. :8.000
##
## spikes.hr min.hr max.hr acute.sz late.sz
## Min. :-133.833 Min. :-50.4955 Min. :-1029.70 0:38 0:20
## 1st Qu.: 5.293 1st Qu.: 0.0000 1st Qu.: 42.25 1: 8 1:26
## Median : 39.101 Median : 0.0000 Median : 255.00
## Mean : 73.543 Mean : -0.8856 Mean : 370.56
## 3rd Qu.: 126.787 3rd Qu.: 4.3046 3rd Qu.: 731.93
## Max. : 416.147 Max. : 43.5564 Max. : 1570.63
##
## ever.sz missing_field.gcs missing_er.gcs missing_icu.gcs
## 0:19 Mode :logical Mode :logical Mode :logical
## 1:27 FALSE:44 FALSE:44 FALSE:45
## TRUE :2 TRUE :2 TRUE :1
##
##
##
##
## missing_worst.gcs missing_X6m.gose missing_spikes.hr missing_min.hr
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:45 FALSE:41 FALSE:28 FALSE:28
## TRUE :1 TRUE :5 TRUE :18 TRUE :18
##
##
##
##
## missing_max.hr
## Mode :logical
## FALSE:28
## TRUE :18
##
##
##
##
##
## [[3]]
## id age sex mechanism
## Min. : 1.00 Min. :16.00 Female: 9 Bike_vs_Auto: 4
## 1st Qu.:12.25 1st Qu.:23.00 Male :37 Blunt : 4
## Median :23.50 Median :33.00 Fall :13
## Mean :23.50 Mean :36.89 GSW : 2
## 3rd Qu.:34.75 3rd Qu.:47.25 MCA : 7
## Max. :46.00 Max. :83.00 MVA :10
## Peds_vs_Auto: 6
## field.gcs er.gcs icu.gcs worst.gcs
## Min. :-5.02 Min. : 3.000 Min. : 0.000 Min. : 0.000
## 1st Qu.: 3.00 1st Qu.: 4.000 1st Qu.: 3.000 1st Qu.: 3.000
## Median : 6.50 Median : 7.500 Median : 6.000 Median : 3.000
## Mean : 7.57 Mean : 8.459 Mean : 6.264 Mean : 5.436
## 3rd Qu.:12.00 3rd Qu.:12.750 3rd Qu.: 7.750 3rd Qu.: 7.051
## Max. :15.00 Max. :25.241 Max. :14.000 Max. :14.000
##
## X6m.gose X2013.gose skull.fx temp.injury surgery
## Min. :-0.2791 Min. :2.000 0:18 0:19 0:17
## 1st Qu.: 3.0000 1st Qu.:5.000 1:28 1:27 1:29
## Median : 5.0000 Median :7.000
## Mean : 4.5770 Mean :5.804
## 3rd Qu.: 6.0000 3rd Qu.:7.000
## Max. : 8.0000 Max. :8.000
##
## spikes.hr min.hr max.hr acute.sz late.sz
## Min. :-211.798 Min. :-19.621 Min. :-616.49 0:38 0:20
## 1st Qu.: 5.965 1st Qu.: 0.000 1st Qu.: 42.25 1: 8 1:26
## Median : 46.665 Median : 0.000 Median : 183.56
## Mean : 80.427 Mean : 4.789 Mean : 322.36
## 3rd Qu.: 151.828 3rd Qu.: 5.533 3rd Qu.: 636.82
## Max. : 458.114 Max. : 53.081 Max. :1497.96
##
## ever.sz missing_field.gcs missing_er.gcs missing_icu.gcs
## 0:19 Mode :logical Mode :logical Mode :logical
## 1:27 FALSE:44 FALSE:44 FALSE:45
## TRUE :2 TRUE :2 TRUE :1
##
##
##
##
## missing_worst.gcs missing_X6m.gose missing_spikes.hr missing_min.hr
## Mode :logical Mode :logical Mode :logical Mode :logical
## FALSE:45 FALSE:41 FALSE:28 FALSE:28
## TRUE :1 TRUE :5 TRUE :18 TRUE :18
##
##
##
##
## missing_max.hr
## Mode :logical
## FALSE:28
## TRUE :18
##
##
##
## 
13.3 Imputation via Expectation-Maximization
Below we present the theory and practice of one specific statistical computing strategy for imputing incomplete datasets.
13.3.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 device that can only measure values above \(0.5\), anything below that is missing.
13.3.2 General Idea of 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:
- A 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 Z 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 parameters, \(\theta\), that maximizes the expectation above, \[\theta^{(t+1)}=\arg\max_{\theta}Q(\theta|\theta^{(t)}).\]
13.3.3 EM-based imputation
The EM algorithm is an alternative to Newton-Raphson or the method of scoring for computing MLE in cases where the complications in calculating the MLE are due to incomplete observation and data are MAR, missing at random, with separate parameters for observation and the missing data mechanism, so the missing data mechanism can be ignored.
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.
- 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.\]
13.3.4 A simple manual implementation of EM-based imputation
# simulate 20 (feature) vectors of 200 (cases) Normal Distributed random values (\mu, \Sigma)
# You can choose multiple distribution for testing
# sim_data <- replicate(20, rpois(50, 10))
set.seed(202227)
mu <- as.matrix(rep(2,20) )
sig <- diag(c(1:20) )
# Add a noise item. The noise is $ \epsilon ~ MVN(as.matrix(rep(0,20)), diag(rep(1,20)))$
sim_data <- mvrnorm(n = 200, mu, sig) +
mvrnorm(n=200, as.matrix(rep(0,20)), diag( rep(1,20) ))
# save these in the "original" object
sim_data.orig <- sim_data
# introduce 500 random missing indices (in the total of 4000=200*20)
# discrete distribution where the probability of the elements of values is proportional to probs,
# which are normalized to add up to 1.
rand.miss <- e1071::rdiscrete(500, probs = rep(1,length(sim_data)), values = seq(1, length(sim_data)))
sim_data[rand.miss] <- NA
sum(is.na(sim_data)) # check now many missing (NA) are there < 500## [1] 459# cast the data into a data.frame object and report 15*10 elements
sim_data.df <- data.frame(sim_data)
kable( sim_data.df[1:15, 1:10], caption = "The first 15 rows and first 10 columns of the simulation data")| X1 | X2 | X3 | X4 | X5 | X6 | X7 | X8 | X9 | X10 |
|---|---|---|---|---|---|---|---|---|---|
| 0.1710559 | 4.8500284 | 1.7282025 | NA | 4.5511409 | 2.1699719 | 4.7118571 | 2.6972905 | NA | -0.9250101 |
| 1.7721898 | -1.6955535 | 0.4234295 | 6.2619764 | -3.1339307 | -1.8043889 | 1.2932913 | 0.9906410 | 1.0790366 | 1.1638161 |
| 2.2859237 | 4.2713494 | 3.8900670 | 3.5104578 | 3.0025607 | 3.9058619 | 3.0548309 | NA | 2.7208548 | 1.3601231 |
| 1.4906345 | 6.9253703 | 1.3387561 | -2.6177675 | -1.4074433 | -2.0861790 | -0.8847125 | 5.6485105 | 1.3348689 | 0.6512754 |
| -0.0305062 | 2.3069337 | NA | 1.2474398 | -0.2091862 | 1.1726298 | 1.4622491 | 0.6492390 | 0.1458781 | 7.2195026 |
| 1.9411674 | 3.6606019 | 6.3004943 | 0.3543817 | 1.0875460 | 1.3077307 | 0.0269935 | -0.0112352 | 2.9724355 | 2.4616466 |
| 0.5209321 | 4.2225983 | -0.6328796 | NA | 5.7103015 | NA | 4.5803136 | 0.2535158 | -0.9758147 | NA |
| 3.3111727 | 0.1245427 | 0.9768425 | 4.7764851 | 4.9810860 | -0.0438375 | -0.3213523 | 4.3612903 | 2.7483497 | 3.7744105 |
| 2.1097813 | 0.2427543 | 1.1718823 | 3.6791032 | 3.7784795 | 5.5533794 | 1.8635898 | 0.1565858 | 4.6743012 | -1.0232277 |
| 1.6203382 | NA | 1.3625543 | 1.4006816 | -1.3467049 | 6.1636769 | 3.1703070 | -3.4308300 | 0.2195447 | -2.2921107 |
| NA | 0.5601779 | NA | 2.1652769 | 0.1060700 | 2.3763668 | -3.7629196 | -0.9718406 | 1.2196606 | 0.7038253 |
| 1.8570634 | 1.6714614 | 0.2000255 | 2.3308968 | 0.4593543 | 0.1716613 | -0.5795159 | -1.3327330 | NA | 1.2607465 |
| 2.5416914 | 1.0388186 | 4.7421120 | 2.3110187 | NA | -1.5291668 | 2.4574501 | NA | 2.5364297 | NA |
| 0.9393925 | 1.8668513 | 3.4702117 | 1.6797620 | 0.7797323 | 2.5548284 | -3.5498813 | 3.6595296 | 3.2124694 | 5.3120133 |
| -0.8334973 | NA | 1.1949414 | -0.6573305 | 1.7484169 | -1.7743495 | 1.7168340 | 2.1065878 | 2.9478528 | -0.5562349 |
# Define the EM imputation method
EM_algorithm <- function(x, tol = 0.001) {
# identify the missing data entries (Boolean indices)
missvals <- is.na(x)
# instantialize the EM-iteration
new.impute <- x
old.impute <- x
count.iter <- 1
reach.tol <- 0
# compute \Sigma on complete data
sigma <- as.matrix(var(na.exclude(x)))
# compute the vector of feature (column) means
mean.vec <- as.matrix(apply(na.exclude(x), 2, mean))
while (reach.tol != 1) {
for (i in 1:nrow(x)) {
pick.miss <- (c(missvals[i, ]))
if (sum(pick.miss) != 0) {
# compute invese-Sigma_completeData, variance-covariance matrix
inv.S <- solve(sigma[!pick.miss, !pick.miss], tol = 1e-40)
# Expectation Step
# $$E(Y|X)=\mu_{mis}+\Sigma_{mo}\Sigma_{oo}^{-1}(X-\mu_{obs})$$
new.impute[i, pick.miss] <- mean.vec[pick.miss] +
sigma[pick.miss,!pick.miss] %*% inv.S %*%
(t(new.impute[i, !pick.miss]) - t(t(mean.vec[!pick.miss])))
}
}
# Maximization Step
# Recompute the complete \Sigma the complete vector of feature (column) means
#$$\Sigma^{(t+1)} = \frac{1}{n}\sum_{i=1}^nE(ZZ^T|X) - \mu^{(t+1)}{\mu^{(t+1)}}^T$$
sigma <- var((new.impute))
#$$\mu^{(t+1)} = \frac{1}{n}\sum_{i=1}^nE(Z|X)$$
mean.vec <- as.matrix(apply(new.impute, 2, mean))
# Inspect for convergence tolerance, start with the 2nd iteration
if (count.iter > 1) {
for (l in 1:nrow(new.impute)) {
for (m in 1:ncol(new.impute)) {
if (abs((old.impute[l, m] - new.impute[l, m])) > tol) {
reach.tol < -0
} else {
reach.tol <- 1
}
}
}
}
count.iter <- count.iter + 1
old.impute <- new.impute
}
# return the imputation output of the current iteration that passed the tolerance level
return(new.impute)
}
sim_data.imputed <- EM_algorithm(sim_data.df, tol=0.0001)13.3.5 Plotting complete and imputed data
Smaller black colored points represent observed data, and magenta-color and circle-shapes denote the imputated 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)
13.3.6 Validation of EM-imputation using the Amelia R Package
13.3.6.1 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)## Loading required package: Rcpp## ##
## ## Amelia II: Multiple Imputation
## ## (Version 1.7.4, built: 2015-12-05)
## ## Copyright (C) 2005-2017 James Honaker, Gary King and Matthew Blackwell
## ## Refer to http://gking.harvard.edu/amelia/ for more information
## ##dim(sim_data.df)## [1] 200 20amelia.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
##
##
## -- Imputation 2 --
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
##
## -- Imputation 3 --
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
##
## -- Imputation 4 --
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 21 22 23 24
##
## -- Imputation 5 --
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17amelia.out##
## Amelia output with 5 imputed datasets.
## Return code: 1
## Message: Normal EM convergence.
##
## Chain Lengths:
## --------------
## Imputation 1: 20
## Imputation 2: 15
## Imputation 3: 16
## Imputation 4: 24
## Imputation 5: 17amelia.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)
13.3.6.2 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 hte \(4,000\) that we synthetically generated.
14 Parsing webpages and visualizing tabular HTML data
In this section, we will utilize the Earthquakes dataset on SOCR website. It records information about earthquakes happened between 1969 and 2007 with magnitudes larger than 5 on the Richter’s 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")##
## Attaching package: 'rvest'## The following object is masked from 'package:XML':
##
## xmlwiki_url <- read_html("http://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-primary" role="main">\n\t<a id="top ...earthquake<- html_table(html_nodes(wiki_url, "table")[[2]])In this dataset, Magt(magnitude type) may be used as 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_line()
plot4 # or plint(plot4)
We can see the most important line of code was made up with 2 parts. The first part ggplot(earthquake, aes(Longiture, 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 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 chunk of codes plots the distribution for Latitude among different Magnitude types. Also, it is using ggplot() function but 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 plot_ly() function under 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. However, we need to register in plotly website to make our own plots. Still, there are some built-in datasets that can be used to demonstrate this type of graph.
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 with() function.
library(plotly)##
## Attaching package: 'plotly'## The following object is masked from 'package:MASS':
##
## select## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutwith(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)## No trace type specified:
## Based on info supplied, a 'histogram' trace seems appropriate.
## Read more about this trace type -> https://plot.ly/r/reference/#histogramplot_ly(x = ~ earthquake$Longitude, y = ~earthquake$Latitude)## No trace type specified:
## Based on info supplied, a 'scatter' trace seems appropriate.
## Read more about this trace type -> https://plot.ly/r/reference/#scatter## No scatter mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plot.ly/r/reference/#scatter-modeplot_ly(x = ~ earthquake$Longitude, y = ~earthquake$Latitude, z=~earthquake$Mag)## No trace type specified:
## Based on info supplied, a 'scatter3d' trace seems appropriate.
## Read more about this trace type -> https://plot.ly/r/reference/#scatter3d## No scatter3d mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plot.ly/r/reference/#scatter-modedf3D <- 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 44View(as.matrix(matrix_EarthQuakes))
# view matrix ias 2D heatmap:
library("ggplot2"); library("gplots")##
## Attaching package: 'gplots'## The following object is masked from 'package:stats':
##
## lowessheatmap.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 ))
# 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")15 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 groups is more than an order of magnitude larger than the size of another group has the potential for bias.
Example 1: Parkinson’s Diseases Study involving 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")
require(crossval)
require(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 misbalance in binary classification.
set.seed(1000)
# install.packages("unbalanced") to deal with unbalanced group data
require(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), dim(output))
#balance the dataset
data.1<-ubBalance(X= input, Y=output, type="ubSMOTE", percOver=300, percUnder=150, verbose=TRUE)## Error in ubBalance(X = input, Y = output, type = "ubSMOTE", percOver = 300, : could not find function "ubBalance"# 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)## Error in cbind(data.1$X, data.1$Y): object 'data.1' not foundtable(data.1$Y)## Error in table(data.1$Y): object 'data.1' not foundnrow(data.1$X); ncol(data.1$X)## Error in nrow(data.1$X): object 'data.1' not found## Error in ncol(data.1$X): object 'data.1' not foundnrow(balancedData); ncol(balancedData)## Error in nrow(balancedData): object 'balancedData' not found## Error in ncol(balancedData): object 'balancedData' not foundnrow(input); ncol(input)
colnames(balancedData) <- c(colnames(input), "PD")## Error in colnames(balancedData) <- c(colnames(input), "PD"): object 'balancedData' not found# check visually for differences between the distributions of the raw (input) and rebalanced data (for only one variable, in this case)
qqplot(input[, 5], balancedData [, 5])## Error in sort(y): object 'balancedData' not found###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]))
}## Error in ncol(balancedData): object 'balancedData' not foundprint(c("Wilcoxon test results: ", test.results.bin))
test.results.corr <- stats::p.adjust(test.results.raw, method = "fdr", n = length(test.results.raw)) ## Error: n >= lp is not TRUE# where methods are "holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none")
plot(test.results.raw, test.results.corr)## Error in plot(test.results.raw, test.results.corr): object 'test.results.corr' not found# zeros (0) are significant independent between-group T-test differences, ones (1) are insignificant
# 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)## Error in ifelse(test.results.corr > alpha.0.05, 1, 0): object 'test.results.corr' not foundtable(test.results.corr.bin)## Error in table(test.results.corr.bin): object 'test.results.corr.bin' not found16 Appendix
16.1 Importing Data from SQL Databases
We can also import SQL databases in to R. First, we need to install and load the RODBC(R Open Database Connectivity) package.
install.packages("RODBC", repos = "http://cran.us.r-project.org")
library(RODBC)Then, we could open a connection to the SQL server database with Data Source Name (DSN), via Microsoft Access. More details are provided here and here.
16.2 R codes for graphs shown in this chapter
#Right Skewed
N <- 10000
x <- rnbinom(N, 10, .5)
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)
#No Skew
N <- 10000
x <- rnorm(N, 0, 1)
hist(x, probability=T,
col='lightblue', xlab=' ', ylab=' ', axes=F,
main='No Skew')
lines(density(x, bw=0.4), col='red', lwd=3)
#Uniform density
x<-runif(1000, 1, 50)
hist(x, col='lightblue', main="Uniform Distribution", probability = T, xlab="", ylab="Density", axes=F)
abline(h=0.02, col='red', lwd=3)
#68-95-99.7 rule
x <- rnorm(N, 0, 1)
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)
