| SOCR ≫ | DSPA ≫ | DSPA2 Topics ≫ |
DSPA2: Data Science and Predictive Analytics (UMich HS650)
Basic Visualization and Exploratory Data Analytics
Basic Visualization and Exploratory Data Analytics
SOCR/MIDAS (Ivo Dinov)
SOCR/MIDAS (Ivo Dinov)
June 2023
In this chapter, we will present a number of complementary strategies for data wrangling, harmonization, manipulation, aggregation, visualization, and graphical exploration. Specifically, we will discuss alternative methods for loading and saving computable data objects, importing and exporting different data structures, measuring sample statistics for quantitative variables, plotting sample histograms and model distribution functions, and scraping data from websites. In addition, we will cover exploratory data analytical (EDA) techniques, handling of incomplete (missing) data, and cohort-rebalancing of imbalanced groups.
1 Data Handling
In this section, we will discuss strategies to import
data and export results. Also, we are going to learn the
basic tricks we need to know about processing different types of data.
Specifically, we will illustrate common R data structures
and strategies for loading (ingesting) and saving (regurgitating) data.
In addition, we will (1) present some basic statistics, e.g., for
measuring central tendency (mean, median, mode) or dispersion (variance,
quartiles, range), (2) explore simple plots, (3) demonstrate the uniform
and normal distributions, (4) contrast numerical and categorical types
of variables, (5) present strategies for handling incomplete (missing)
data, and (6) show the need for cohort-rebalancing when comparing
imbalanced groups of subjects, cases or units.
1.1 Saving and Loading
R Data Structures
Let’s start by extracting Edgar Anderson’s Iris Data from the package
datasets. The iris
dataset quantifies morphologic shape variations of 50 Iris flowers
of three related genera - Iris setosa, Iris virginica
and Iris versicolor. Four shape features were measured from
each sample - length and the width of the sepals and petals (in
centimeters). These data were used by Ronald Fisher in
his 1936
linear discriminant analysis paper.
## [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.
1.2 Importing and Saving Data from CSV Files
Importing the data from
"CaseStudy07_WorldDrinkingWater_Data.csv" from these
case-studies and saving it into the R dataset named
“water”. The variables in the dataset are as follows:
- Time: Years (1990, 1995, 2000, 2005, 2010, 2012)
- Demographic: Country (across the world)
- Residence Area Type: Urban, rural, or total
- WHO Region
- Population using improved drinking-water sources: The percentage of the population using an improved drinking water source.
- Population using improved sanitation facilities: The percentage of the population using an improved sanitation facility.
Generally, the separator of a CSV file is comma. By default, we have
optionsep=", " in the command read.csv().
Also, we can use colnames() to rename the column variables.
Let’s use CaseStudy07_WorldDrinkingWater_Data.csv
from out Canvas
Data Archive as an example. This code loads CSV files that already
include a header line containing the names of the variables. If we don’t
have a header in the dataset, we can use the header = FALSE
option to read the first row in the file as data. In such cases,
R will assign default names to the column variables of the
dataset.
water <- read.csv(
'https://umich.instructure.com/files/399172/download?download_frd=1',
header=TRUE, fileEncoding = "UTF-8") #, fileEncoding = "UTF-8")
water[1:3, ]## Year..string. WHO.region..string. Country..string.
## 1 1990 Africa Algeria
## 2 1990 Africa Angola
## 3 1990 Africa Benin
## Residence.Area.Type..string.
## 1 Rural
## 2 Rural
## 3 Rural
## Population.using.improved.drinking.water.sources......numeric.
## 1 88
## 2 42
## 3 49
## Population.using.improved.sanitation.facilities......numeric.
## 1 77
## 2 7
## 3 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 0## [1] 1## [1] 20.4To save a data frame to CSV files, we could use the
write.csv() function. The option
file = "a/local/file/path" allows us to specify the output
file name and location.
1.3 Importing Data from ZIP and SAV Files
This example demonstrates data import from a compressed (ZIP) SPSS (SAV) file. In this case, we utilize DSPA Case-Study 25: National Ambulatory Medical Care Survey (NAMCS).
# install.packages("foreign")
library("foreign")
pathToZip <- tempfile()
download.file("https://umich.instructure.com/files/8111611/download?download_frd=1", pathToZip, mode = "wb")
dataset <- read.spss(unzip(pathToZip, files = "namcs2015-spss.sav", list = F, overwrite = TRUE), to.data.frame=TRUE)
dim(dataset)## [1] 28332 10961.4 Exploring the Structure of Data
We can use the command str() and describe()
to explore the structure of a dataset (in this case the
CaseStudy07_WorldDrinkingWater_Data).
## 'data.frame': 11 obs. of 6 variables:
## $ year : int 1990 1990 1990 1990 1990 1990 1990 1990 1990 1990 ...
## $ region : chr "Africa" "Africa" "Africa" "Africa" ...
## $ country : chr "Algeria" "Angola" "Benin" "Botswana" ...
## $ residence_area : chr "Rural" "Rural" "Rural" "Rural" ...
## $ improved_water : int 88 42 49 86 39 67 34 46 37 83 ...
## $ sanitation_facilities: int 77 7 0 22 2 42 27 12 4 11 ...## water
##
## 6 Variables 11 Observations
## --------------------------------------------------------------------------------
## year
## n missing distinct Info Mean Gmd
## 11 0 1 0 1990 0
##
## Value 1990
## Frequency 11
## Proportion 1
## --------------------------------------------------------------------------------
## region
## n missing distinct value
## 11 0 1 Africa
##
## Value Africa
## Frequency 11
## Proportion 1
## --------------------------------------------------------------------------------
## country
## n missing distinct
## 11 0 11
##
## lowest : Algeria Angola Benin Botswana Burkina Faso
## highest: C Cameroon Central African Republic Chad Comoros
## --------------------------------------------------------------------------------
## residence_area
## n missing distinct value
## 10 1 1 Rural
##
## Value Rural
## Frequency 10
## Proportion 1
## --------------------------------------------------------------------------------
## improved_water
## n missing distinct Info Mean Gmd .05 .10
## 10 1 10 1 57.1 25.04 35.35 36.70
## .25 .50 .75 .90 .95
## 39.75 47.50 79.00 86.20 87.10
##
## Value 34.00 36.70 38.86 41.56 45.88 48.58 66.94 82.60 85.84 88.00
## Frequency 1 1 1 1 1 1 1 1 1 1
## Proportion 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
##
## For the frequency table, variable is rounded to the nearest 0.54
## --------------------------------------------------------------------------------
## sanitation_facilities
## n missing distinct Info Mean Gmd .05 .10
## 10 1 10 1 20.4 25.2 0.90 1.80
## .25 .50 .75 .90 .95
## 4.75 11.50 25.75 45.50 61.25
##
## Value 0.00 1.54 3.85 6.93 10.78 11.55 21.56 26.95 41.58 77.00
## Frequency 1 1 1 1 1 1 1 1 1 1
## Proportion 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
##
## For the frequency table, variable is rounded to the nearest 0.77
## --------------------------------------------------------------------------------We can see that this World Drinking Water dataset has
3331 observations and 6 variables. The output also includes the class of
each variable and first few elements in the variable. The dimension of
the other dataset (Case-Study 25: National Ambulatory Medical Care
Survey) is much larger, \(28,332\times
1,096\).
1.5 Exploring Numeric Variables
Summary statistics for numeric variables in the dataset could be
accessed by using the command summary().
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1990 1990 1990 1990 1990 1990## improved_water sanitation_facilities
## Min. :34.00 Min. : 0.00
## 1st Qu.:39.75 1st Qu.: 4.75
## Median :47.50 Median :11.50
## Mean :57.10 Mean :20.40
## 3rd Qu.:79.00 3rd Qu.:25.75
## Max. :88.00 Max. :77.00
## NA's :1 NA's :1# plot(density(water$improved_water,na.rm = T)) # no need to be continuous, we can still get intuition about the variable distribution
fit <- density(as.numeric(water$improved_water),na.rm = T)
plot_ly(x = fit$x, y = fit$y, type = "scatter", mode = "lines",
fill = "tozeroy", name = "Density") %>%
layout(title='Density of (%) Improved Water Quality',
xaxis = list (title = 'Percent'), yaxis = list (title = 'Density'))The six summary statistics and NA’s (missing data) are
reported in the output.
1.6 Measuring Central Tendency - mean, median, mode
Mean and median are two frequent
measurements of the central tendency. Mean is “the sum of all values
divided by the number of values”. Median is the number in the middle of
an ordered list of values. In R, mean() and
median() functions can provide us with these two
measurements.
## [1] 65## [1] 65## [1] 56## [1] 56## [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 represent an example of a multimodal variable that has two, or more, modes.
Mode is one of the measures for the central tendency. The best way to use it is to compare the mode to other values in the data. This helps us determine whether one or several categories dominate all others in the data. In numeric datasets, we could think mode as the highest bin in the histogram, since it is unlikely to have many repeated measurements for continuous variables. In this way, we can also examine if the numeric data is multimodal.
More information about measures of centrality is available here.
1.7 Measuring Spread - variance, quartiles and the five-number summary
The five-number summary describes the spread of a dataset. They are:
- Minimum (
Min.), representing the smallest value in the data - First quartile/Q1 (
1st Qu.), representing the \(25^{th}\) percentile, which splits off the lowest 25% of data from the highest 75% - Median/Q2 (
Median), representing the \(50^{th}\) percentile, which splits off the lowest 50% of data from the top 50% - Third quartile/Q3 (
3rd Qu.), representing the \(75^{th}\) percentile, which splits off the lowest 75% of data from the top 25% - Maximum (
Max.), representing the largest value in the data.
Min and Max can be obtained by using
min() and max() respectively.
The difference between maximum and minimum is known as range. In R,
range() function gives us both the minimum and maximum. A
combination of range() and diff() could do the
trick of getting the actual range value. To avoid problems with missing
values, we will ignore them using the option
na.rm=TRUE.
## [1] 34 88## [1] 54Q1 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.
## [1] 39.25## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 34.00 39.75 47.50 57.10 79.00 88.00 1Just 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.
## 0% 25% 50% 75% 100%
## 34.00 39.75 47.50 79.00 88.00We can also calculate specific percentiles in the data. For example, if we want the 20th and 60th percentiles, we can do the following.
## 20% 60%
## 38.6 56.2When we include the seq() function, generating
percentiles of evenly-spaced values is available.
## 0% 20% 40% 60% 80% 100%
## 34.0 38.6 44.4 56.2 83.6 88.0Let’s re-examine the five-number summary for the
improved_water variable. When we ignore the
NA’s, the difference between minimum and Q1 is 74 while the
difference between Q3 and maximum is only 1. The interquartile range is
22%. Combining these facts, the first quarter is more widely spread than
the middle 50 percent of values. The last quarter is the most condensed
one that has only two percentages 99% and 100%. Also, we can notice that
the mean is smaller than the median. The mean is more sensitive to the
extreme values than the median. Having some of very small values may
spread out the first quartile, skew the distribution to the left and
make the mean less than the median.
Distribution models offer a way to characterize data using only a few parameters. For example, the normal distribution can be defined by only two parameters - center and spread, statistically speaking, mean and standard deviation.
The mean value is obtained by arithmetic averaging of all data points.
\[mean(X)=\mu=\frac{1}{n}\sum_{i=1}^{n} x_i\]
The standard deviation is the square root of the variance. And the variance is the average sum of square deviation from the mean.
\[Var(X)=\sigma^2=\frac{1}{n-1}\sum^{n}_{i=1} (x_i-\mu)^2\] \[StdDev(X)=\sigma=\sqrt{Var(X)}\]
Since the water dataset is not close to normal, in this example, we will use MLB baseball players dataset to illustrate normal distribution properties. The MLB dataset (01a_data.txt) in our class file data archive has following variables - Name, Team, Position, Height, Weight, and Age.
We can use histograms to visually assess approximate normality of baseball players’ Height and Weight.
baseball<-read.table("https://umich.instructure.com/files/330381/download?download_frd=1", header=T)
# hist(baseball$Weight, main = "Histogram for Baseball Player's Weight", xlab="weight")
# hist(baseball$Height, main = "Histogram for Baseball Player's Height", xlab="height")
x <- rnorm(10000, mean=mean(baseball$Weight, na.rm=T), sd=sd(baseball$Weight, na.rm=T))
fit <- density(x, bw=10)
plot_ly(x=~baseball$Weight, type = "histogram", name = "Weight Histogram", histnorm = "probability") %>%
add_trace(x =~fit$x, y =~5*fit$y, type = "scatter", mode = "lines", opacity=0.1,
fill = "tozeroy", name = "Normal Density") %>%
layout(title='Baseball Weight Histogram & Model Normal Distribution',
xaxis = list(title = "Weight"), yaxis = list(title = "relative frequency/density"),
legend = list(orientation = 'h')) x <- rnorm(10000, mean=mean(baseball$Height, na.rm=T), sd=sd(baseball$Height, na.rm=T))
fit <- density(x, bw=1)
plot_ly(x=~baseball$Height, type = "histogram", name = "Height Histogram", histnorm = "probability") %>%
add_trace(x=~fit$x, y=~fit$y, type = "scatter", mode = "lines", opacity=0.1,
fill = "tozeroy", name = "Normal Density") %>%
layout(title='Baseball Height Histogram & Model Normal Distribution',
xaxis = list(title = "Height"), yaxis = list(title = "relative frequency/density"),
legend = list(orientation = 'h')) TWe could also report the mean and standard deviation of the weight and height variables.
## [1] 201.7166## [1] 73.69729## [1] 440.9913## [1] 20.99979## [1] 5.316798## [1] 2.305818Larger standard deviation, or variance, suggests the data is more spread out from the mean. Therefore, for MLB players, weights appear to be more spread than heights.
Given the first two moments (mean and standard deviation), we can
easily estimate how extreme a specific value is. Assuming we have a
normal distribution, the values follow a \(68-95-99.7\) rule. This means 68% of the
data lies within the interval \([\mu-\sigma,
\mu+\sigma]\);95% of the data lies within the interval \([\mu-2*\sigma, \mu+2*\sigma]\) and 99.7% of
the data lies within the interval \([\mu-3*\sigma, \mu+3*\sigma]\). The
following graph plotted by R illustrates the \(68-95-99.7\) rule.
# hist(x, probability=T,
# col='lightblue', xlab=' ', ylab=' ', axes = F,
# main='68-95-99.7 Rule')
# lines(density(x, bw=0.4), col='red', lwd=3)
# axis(1, at=c(-3, -2, -1, 0, 1, 2, 3), labels = expression(mu-3*sigma, mu-2*sigma, mu-sigma, mu, mu+sigma, mu+2*sigma, mu+3*sigma))
# abline(v=-1, lwd=3, lty=2)
# abline(v=1, lwd=3, lty=2)
# abline(v=-2, lwd=3, lty=2)
# abline(v=2, lwd=3, lty=2)
# abline(v=-3, lwd=3, lty=2)
# abline(v=3, lwd=3, lty=2)
# text(0, 0.2, "68%")
# segments(-1, 0.2, -0.3, 0.2, col = 'red', lwd=2)
# segments(1, 0.2, 0.3, 0.2, col = 'red', lwd=2)
# text(0, 0.15, "95%")
# segments(-2, 0.15, -0.3, 0.15, col = 'red', lwd=2)
# segments(2, 0.15, 0.3, 0.15, col = 'red', lwd=2)
# text(0, 0.1, "99.7%")
# segments(-3, 0.1, -0.3, 0.1, col = 'red', lwd=2)
# segments(3, 0.1, 0.3, 0.1, col = 'red', lwd=2)
N<- 1000
norm <- rnorm(N, 0, 1)
# hist(x, probability=T,
# col='lightblue', xlab=' ', ylab=' ', axes=F,
# main='Normal Distribution')
# lines(density(x, bw=0.4), col='red', lwd=3)
normDensity <- density(norm, bw=0.5)
dens <- data.frame(x = normDensity$x, y = normDensity$y)
miny <- 0
maxy <- max(dens$y)
xLabels <- c("μ-3σ","μ-2σ", "μ-σ", "μ", "μ+σ", "μ+2σ", "μ+3σ")
labelColors <- c("green", "red", "orange", "black", "orange", "red", "green")
xLocation <- c(-3, -2, -1, 0, 1, 2, 3)
yLocation <- 0.2
data <- data.frame(xLabels, xLocation, yLocation)
plot_ly(dens) %>%
add_histogram(x = norm, name="Normal Histogram") %>%
add_lines(data = dens, x = ~x, y = ~y+0.05, yaxis = "y2",
line = list(width = 3), name="N(0,1)") %>%
add_annotations(x = ~xLocation, y = ~yLocation, type = 'scatter', ax = 20, ay = 20,
mode = 'text', text = ~xLabels, textposition = 'middle right',
textfont = list(color = labelColors, size = 16)) %>%
add_segments(x=-3, xend=-3, y=0, yend=100, name="99.7%", line=list(dash="dash", color="green")) %>%
add_segments(x=-2, xend=-2, y=0, yend=90, name="95%", line=list(dash="dash", color="red")) %>%
add_segments(x=-1, xend=-1, y=0, yend=80, name="68%", line=list(dash="dash", color="orange")) %>%
add_segments(x=1, xend=1, y=0, yend=80, name="68%", line = list(dash = "dash", color="orange")) %>%
add_segments(x=2, xend=2, y=0, yend=90, name="95%", line=list(dash="dash", color="red")) %>%
add_segments(x=3, xend=3, y=0, yend=100, name="99.7%", line=list(dash="dash", color="green")) %>%
add_segments(x=-3, xend=3, y=100, yend=100, name="99.7%", line=list(dash="dash", color="green")) %>%
add_segments(x=-2, xend=2, y=90, yend=90, name="95%", line=list(dash="dash", color="red")) %>%
add_segments(x=-1, xend=1, y=80, yend=80, name="68%", line=list(dash="dash", color="orange")) %>%
layout(bargap=0.1, xaxis=list(name=""), yaxis=list(title="density/frequency"),
yaxis2 = list(overlaying = "y", side = "right", # title="prob",
range = c(miny, maxy+0.1), showgrid = F, zeroline = F),
legend = list(orientation = 'h'), title="Normal 68-95-99.7% Rule")Applying the 68-95-99.7 rule to our baseball weight variable, we know that 68% of our players weighted between 180.7168 pounds and 222.7164 pounds; 95% of the players weighted between 159.7170 pounds and 243.7162 pounds; And 99.7% of the players weighted between 138.7172 pounds and 264.7160 pounds.
1.8 Visualizing Numeric Variables - boxplots
We can visualize the five-number summary by a boxplot
(box-and-whiskers plot). With the boxplot() function we can
manage the title (main="") and labels for x
(xlab="") and y (ylab="") axis.
# boxplot(water$improved_water, main="Boxplot for Percent improved_water", ylab="Percentage")
plot_ly(y = ~water$improved_water, type = "box", name="improved water qual") %>%
add_trace(y = ~water$sanitation_facilities, name ="sanitation") %>%
layout(title='Boxplots of Improved Water Quality and Sanitation Facilities',
yaxis = list (title = 'Percent'))In the boxplot we have five horizontal lines each representing the corresponding value in the five-number summary. The box in the middle represents the middle 50 percent of values. The bold line in the box is the median. Mean value is not illustrated on the graph.
Boxplots only allow the two ends to extend to a minimum or maximum of 1.5 times the IQR. Therefore, any value that falls outside of the \(3\times IQR\) range will be represented as circles or dots. They are considered outliers. We can see that there are a lot of outliers with small values on the low ends of the graph.
1.9 Visualizing Numeric Variables - histograms
Histograms offer another way to show the distribution spread of numeric variables. They require a specification of a number of bins, value containers, to divide and stratify the original data. The heights of the bins indicate the observed frequencies within each bin.
# hist(water$improved_water, main = "Histogram of Percent improved_water", xlab="Percentage")
# hist(water$sanitation_facilities, main = "Histogram of Percent sanitation_facilities", xlab = "Percentage")
plot_ly(x = ~water$improved_water, type = "histogram", name="improved_water") %>%
add_trace(x = ~water$sanitation_facilities, type = "histogram", name="sanitation_facilities") %>%
layout(bargap=0.1, title='Histograms', legend = list(orientation = 'h'),
xaxis = list(title = 'Percent'), yaxis = list (title = 'Frequency'))We could see that the shape of two graphs are somewhat similar. They both appear to have left skewed patterns (\(mean \lt median\)). Other common skew patterns are shown in the following graph.
N <- 10000
x <- rnbinom(N, 5, 0.1)
# hist(x,
# xlim=c(min(x), max(x)), probability=T, nclass=max(x)-min(x)+1,
# col='lightblue', xlab=' ', ylab=' ', axes=F,
# main='Right Skewed')
# lines(density(x, bw=1), col='red', lwd=3)
fit <- density(x)
plot_ly(x = x, type = "histogram", name = "Data Histogram") %>%
add_trace(x = fit$x, y = fit$y, type = "scatter", mode = "lines", opacity=0.3,
fill = "tozeroy", yaxis = "y2", name = "Density (rnbinom(N, 5, 0.1))") %>%
layout(title='Right Skewed Process', yaxis2 = list(overlaying = "y", side = "right"),
legend = list(orientation = 'h'))N <- 10000
x <- rnorm(N, 15, 3.7)
# hist(x,
# xlim=c(min(x), max(x)), probability=T, nclass=max(x)-min(x)+1,
# col='lightblue', xlab=' ', ylab=' ', axes=F,
# main='Right Skewed')
# lines(density(x, bw=1), col='red', lwd=3)
fit <- density(x)
plot_ly(x = x, type = "histogram", name = "Data Histogram") %>%
add_trace(x = fit$x, y = fit$y, type = "scatter", mode = "lines", opacity=0.3,
fill = "tozeroy", yaxis = "y2", name = "Density (rnorm(N, 15, 3.7))") %>%
layout(title='Symmetric Process', yaxis2 = list(overlaying = "y", side = "right"),
legend = list(orientation = 'h'))# N <- 10000
# xNu <- extraDistr::rlaplace(N, mu = 0, sigma = 0.4)
# yNu <- density(xNu, bw=0.2)
# xMu <- extraDistr::rlaplace(N, mu = 0, sigma = 0.5)
# yMu <- density(xMu, bw=0.2)
# # correct second Laplace Density (mu) to ensure absolute continuity, nu << mu
# yMu$y <- 2*yMu$y
# plot_ly(x = x, type = "histogram", name = "Data Histogram") %>%
# add_trace(x = yNu$x, y = yNu$y, type = "scatter", mode = "lines", opacity=0.3,
# fill = "tozeroy", yaxis = "y2", name = "nu, Laplace(N,0,0.4) Density") %>%
# add_trace(x = yMu$x, y = yMu$y, type="scatter", mode="lines", opacity=0.3,
# fill = "tozeroy", yaxis = "y2", name = "mu, Laplace(N,0,0.5) Density") %>%
# layout(title="Absolutely Continuous Laplace Distributions, nu<<mu",
# yaxis2 = list(overlaying = "y", side = "right"),
# xaxis = list(range = list(-pi, pi)),
# legend = list(orientation = 'h'))
# integrate(approxfun(yNu), -pi, pi)
# integrate(approxfun(yMu), -pi, pi)You can learn more about Probability Distributions in the SOCR EBook and see the density plots of over 80 different probability distributions using the SOCR Java Distribution Calculators or the Distributome HTML5 Distribution Calculators.
For each probability distribution defined in R, there
are four functions that provide the density (e.g., dnorm),
the cumulative probability (e.g., pnorm), the inverse
cumulative distribution (quantile) function (e.g., qnorm),
and the random sampling (simulation) function (e.g.,
rnorm). The plots below show the standard normal
density, cumulative probability and the quantile functions. As the
density is very small outside of the interval \((-4,4)\), the plots are restricted to this
domain.
z<-seq(-4, 4, 0.1) # points from -4 to 4 in 0.1 steps
q<-seq(0.001, 0.999, 0.001) # probability quantile values from 0.1% to 99.9% in 0.1% steps
dStandardNormal <- data.frame(Z=z, Density=dnorm(z, mean=0, sd=1), Distribution=pnorm(z, mean=0, sd=1))
qStandardNormal <- data.frame(Q=q, Quantile=qnorm(q, mean=0, sd=1))
head(dStandardNormal)## Z Density Distribution
## 1 -4.0 0.0001338302 3.167124e-05
## 2 -3.9 0.0001986555 4.809634e-05
## 3 -3.8 0.0002919469 7.234804e-05
## 4 -3.7 0.0004247803 1.077997e-04
## 5 -3.6 0.0006119019 1.591086e-04
## 6 -3.5 0.0008726827 2.326291e-04# plot(z, dStandardNormal$Density, main="Normal Density Curve", type = "l", xlab = "critical values", ylab="density", lwd=4, col="blue")
# polygon(z, dStandardNormal$Density, col="red", border="blue")
# plot(z, dStandardNormal$Distribution, main="Normal Distribution", type = "l", xlab = "critical values", ylab="Cumulative Distribution", lwd=4, col="blue")
# plot(q, qStandardNormal$Quantile, main="Normal Quantile Function (Inverse CDF)", type = "l", xlab = "p-values", ylab="Critical Values", lwd=4, col="blue")
plot_ly(x = z, y= dStandardNormal$Density, name = "Normal Density Curve",
mode = 'lines') %>%
layout(title='Normal Density Curve',
xaxis = list(title = 'critical values'),
yaxis = list(title ="Density"),
legend = list(orientation = 'h'))1.10 Understanding Numeric Data - uniform and normal distributions
If the data follows a uniform distribution, then all values are equally likely to occur. The histogram for a uniformly distributed data would have equal heights for each bin like the following graph.
Often, but not always, real world processes may appear as normally distributed data. A normal distribution would have a higher frequency for middle values and lower frequency for more extreme values. It has a symmetric and bell-curved shape just like the following diagram generated by R. Many parametric-based statistical approaches assume normality of the data. In cases where this parametric assumption is violated, variable transformations or distribution-free tests may be more appropriate.
N<- 1000
norm <- rnorm(N, 0, 1)
# hist(x, probability=T,
# col='lightblue', xlab=' ', ylab=' ', axes=F,
# main='Normal Distribution')
# lines(density(x, bw=0.4), col='red', lwd=3)
normDensity <- density(norm, bw=0.5)
dens <- data.frame(x = normDensity$x, y = normDensity$y)
miny <- 0
maxy <- max(dens$y)
plot_ly(dens) %>%
add_histogram(x = norm, name="Normal Histogram") %>%
add_lines(data = dens, x = ~x, y = ~y, yaxis = "y2",
line = list(width = 3), name="N(0,1)") %>%
layout(bargap=0.1, yaxis2 = list(overlaying = "y", side = "right",
range = c(miny, maxy), showgrid = F, zeroline = F),
legend = list(orientation = 'h'), title="Normal(0,1)")1.11 Exploring Categorical Variables
Back to our water dataset, we can treat the year variable as categorical rather than a numeric variable. Since the year variable only has six distinctive values, it is rational to treat it as a categorical variable where each value is a category that could apply to multiple WHO regions. Moreover, region and residence area variables are also categorical.
Different from numeric variables, the categorical variables are
better examined by tables rather than summary statistics. One-way table
represents a single categorical variable. It gives us the counts of
different categories. table() function can create one-way
tables for our water dataset:
water <- read.csv('https://umich.instructure.com/files/399172/download?download_frd=1', header=T, stringsAsFactors=FALSE, fileEncoding="latin1")
colnames(water)<-c("year", "region", "country", "residence_area", "improved_water", "sanitation_facilities")
table(water$year)##
## 1990 1995 2000 2005 2010 2012
## 520 561 570 570 556 554##
## Africa Americas Eastern Mediterranean
## 797 613 373
## Europe South-East Asia Western Pacific
## 910 191 447##
## 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.
##
## 1990 1995 2000 2005 2010 2012
## 0.1561093 0.1684179 0.1711198 0.1711198 0.1669168 0.1663164##
## 1990 1995 2000 2005 2010 2012
## 15.6 16.8 17.1 17.1 16.7 16.61.12 Exploring Relationships Between Variables
So far, the methods and statistics that we have gone through are at univariate level. Sometimes we want to examine the relationship between two or multiple variables. For example, does the percentage of population that uses improved drinking-water sources increase over time? To address these problems we need to look at bivariate or multivariate relationships.
1.12.1 Visualizing Relationships - scatterplots
Let’s look at the bivariate case first. A scatterplot is a good way to visualize bivariate relationships. We have x axis and y axis each representing one of the variables. Each observation is illustrated on the graph by a glyph, e.g., a solid point. If the graph shows a clear pattern, rather than a random scatter of points or a horizontal line, the two variables may be correlated with each other.
In R we can use the plot() function to
create scatterplots. We have to define the variables for x-axis and
y-axis. The labels in the graph are editable.
# plot.window(c(400,1000), c(500,1000))
# plot(x=water$year, y=water$improved_water,
# main= "Scatterplot of Year vs. Improved_water",
# xlab= "Year",
# ylab= "Percent of Population Using Improved Water")
plot_ly(x = ~water$sanitation_facilities, y = ~water$improved_water, type = "scatter",
mode = "markers") %>%
layout(title='Scatterplot: Improved Water Quality vs. Sanitation Facilities',
xaxis = list (title = 'Water Quality'), yaxis = list (title = 'Sanitation'))We can see from the scatterplot that there is an increasing pattern. In later years, the percentages are more centered around one hundred. Especially, in 2012, none of the regions had less than 20% of people using improved water sources while there used to be some regions that had such low percentages in the early years.
1.12.2 Examining Relationships - two-way cross-tabulations
Scatterplot is a useful tool to examine the relationship between two variables where at least one of them is numeric. When both variables are nominal, two-way cross-tabulation would be a better choice (also named as crosstab or contingency table).
The function CrossTable() is available in R
under the package gmodels. Let’s install it first.
We are interested in investigating the relationship between WHO region and residence area type in the water study. We might want to know if there is a difference in terms of residence area type between the African WHO region and all other WHO regions.
To address this problem we need to create an indicator variable for the African WHO region first.
Let’s revisit the table() function to see how many WHO
regions are in Africa.
##
## FALSE TRUE
## 2534 797Now, let’s create a two-way cross-tabulation using
CrossTable().
##
##
## Cell Contents
## |-------------------------|
## | N |
## | Chi-square contribution |
## | N / Row Total |
## | N / Col Total |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 3331
##
##
## | water$africa
## water$residence_area | FALSE | TRUE | Row Total |
## ---------------------|-----------|-----------|-----------|
## Rural | 828 | 267 | 1095 |
## | 0.030 | 0.096 | |
## | 0.756 | 0.244 | 0.329 |
## | 0.327 | 0.335 | |
## | 0.249 | 0.080 | |
## ---------------------|-----------|-----------|-----------|
## Total | 845 | 264 | 1109 |
## | 0.002 | 0.007 | |
## | 0.762 | 0.238 | 0.333 |
## | 0.333 | 0.331 | |
## | 0.254 | 0.079 | |
## ---------------------|-----------|-----------|-----------|
## Urban | 861 | 266 | 1127 |
## | 0.016 | 0.050 | |
## | 0.764 | 0.236 | 0.338 |
## | 0.340 | 0.334 | |
## | 0.258 | 0.080 | |
## ---------------------|-----------|-----------|-----------|
## Column Total | 2534 | 797 | 3331 |
## | 0.761 | 0.239 | |
## ---------------------|-----------|-----------|-----------|
##
## Each cell in the table contains five numbers. The first one N gives us the count that falls into its corresponding category. The Chi-square contribution provides us information about the cell’s contribution in the Pearson’s Chi-squared test for independence between two variables. This number measures the probability that the differences in cell counts are due to chance alone.
The number of most interest is the N/ Col Total or the
counts over column total. In this case, these numbers represent the
distribution for residence area type among African regions and the
regions in the rest of the world. We can see the numbers are very close
between African and non-African regions for each type of residence area.
Therefore, we can conclude that African WHO regions do not have a
difference in terms of residence area types compared to the rest of the
world.
1.13 Missing Data
In the previous sections, we simply ignored the missing observations
in our water dataset (na.rm = TRUE). Is this an appropriate
strategy to handle incomplete data? Could the missingness pattern of
those incomplete observations be important? It is possible that the
arrangement of the missing observations may reflect an important factor
that was not accounted for in our statistics or our models.
Missing Completely at Random (MCAR) is an assumption about the probability of missingness being equal for all cases; Missing at Random (MAR) assumes the probability of missingness has a known but random mechanism (e.g., different rates for different groups); Missing not at Random (MNAR) suggest a missingness mechanism linked to the values of predictors and/or response, e.g., some participants may drop out of a drug trial when they have side-effects.
There are a number of strategies to impute missing data. The expectation maximization (EM) algorithm provides one example for handling missing data. The SOCR EM tutorial, activity, and documentations provides the theory, applications and practice for effective (multidimensional) EM parameter estimation.
The simplest way to handle incomplete data is to substitute each missing value with its (feature or column) average. When the missingness proportion is small, the effect of substituting the means for the missing values will have little effect on the mean, variance, or other important statistics of the data. Also, this will preserve those non-missing values of the same observation or row.
m1<-mean(water$improved_water, na.rm = T)
m2<-mean(water$sanitation_facilities, na.rm = T)
water_imp<-water
for(i in 1:3331){
if(is.na(water_imp$improved_water[i])){
water_imp$improved_water[i] <- m1
}
if(is.na(water_imp$sanitation_facilities[i])){
water_imp$sanitation_facilities[i] <- m2
}
}
summary(water_imp)## year region country residence_area
## Min. :1990 Length:3331 Length:3331 Length:3331
## 1st Qu.:1995 Class :character Class :character Class :character
## Median :2005 Mode :character Mode :character Mode :character
## Mean :2002
## 3rd Qu.:2010
## Max. :2012
## improved_water sanitation_facilities africa
## Min. : 3.0 Min. : 0.00 Mode :logical
## 1st Qu.: 77.0 1st Qu.: 44.00 FALSE:2534
## Median : 93.0 Median : 79.00 TRUE :797
## Mean : 84.9 Mean : 68.87
## 3rd Qu.: 99.0 3rd Qu.: 97.00
## Max. :100.0 Max. :100.00A more sophisticated way of resolving missing data is to use a model
(e.g., linear regression) to predict the missing feature and impute its
missing values. This is called the
predictive mean matching approach. This method is good for
data with multivariate normality. However, a disadvantage of it is that
it can only predict one value at a time, which is very time consuming.
Also, the multivariate normality assumption might not be satisfied and
there may be important multivariate relations that are not accounted
for. We are using the mi package for the predictive mean
matching procedure.
Let’s install the mi package first.
Then we need to get the missing information matrix. We are using the
imputation method pmm(predictive mean matching approach)
for both missing variables.
## year region country residence_area improved_water sanitation_facilities
## 1 1990 Africa Algeria Rural 88 77
## 2 1990 Africa Angola Rural 42 7
## 3 1990 Africa Benin Rural 49 0
## 4 1990 Africa Botswana Rural 86 22
## 5 1990 Africa Burkina Faso Rural 39 2
## 6 1990 Africa Burundi Rural 67 42
## africa missing_improved_water missing_sanitation_facilities
## 1 TRUE FALSE FALSE
## 2 TRUE FALSE FALSE
## 3 TRUE FALSE FALSE
## 4 TRUE FALSE FALSE
## 5 TRUE FALSE FALSE
## 6 TRUE FALSE FALSE## Object of class missing_data.frame with 3331 observations on 7 variables
##
## There are 3 missing data patterns
##
## Append '@patterns' to this missing_data.frame to access the corresponding pattern for every observation or perhaps use table()
##
## type missing method model
## year continuous 0 <NA> <NA>
## region unordered-categorical 0 <NA> <NA>
## country unordered-categorical 0 <NA> <NA>
## residence_area unordered-categorical 0 <NA> <NA>
## improved_water continuous 32 ppd linear
## sanitation_facilities continuous 135 ppd linear
## africa binary 0 <NA> <NA>
##
## family link transformation
## year <NA> <NA> standardize
## region <NA> <NA> <NA>
## country <NA> <NA> <NA>
## residence_area <NA> <NA> <NA>
## improved_water gaussian identity standardize
## sanitation_facilities gaussian identity standardize
## africa <NA> <NA> <NA>mdf<-change(mdf, y="improved_water", what = "imputation_method", to="pmm")
mdf<-change(mdf, y="sanitation_facilities", what = "imputation_method", to="pmm")Notes:
Converting the input
data.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 three times, which will create three different (complete) datasets, three chains, with slightly different imputed values.
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.
## year region country residence_area
## Min. :1990 Length:3331 Length:3331 Length:3331
## 1st Qu.:1995 Class :character Class :character Class :character
## Median :2005 Mode :character Mode :character Mode :character
## Mean :2002
## 3rd Qu.:2010
## Max. :2012
##
## improved_water sanitation_facilities africa
## Min. : 3.0 Min. : 0.00 Mode :logical
## 1st Qu.: 77.0 1st Qu.: 42.00 FALSE:2534
## Median : 93.0 Median : 81.00 TRUE :797
## Mean : 84.9 Mean : 68.87
## 3rd Qu.: 99.0 3rd Qu.: 97.00
## Max. :100.0 Max. :100.00
## NA's :32 NA's :135## year region country
## Min. :1990 Africa :797 Albania : 18
## 1st Qu.:1995 Americas :613 Algeria : 18
## Median :2005 Eastern Mediterranean:373 Andorra : 18
## Mean :2002 Europe :910 Angola : 18
## 3rd Qu.:2010 South-East Asia :191 Antigua and Barbuda: 18
## Max. :2012 Western Pacific :447 Argentina : 18
## (Other) :3223
## residence_area improved_water sanitation_facilities africa
## Rural:1095 Min. : 3.00 Min. : 0.00 FALSE:2534
## Total:1109 1st Qu.: 77.00 1st Qu.: 43.00 TRUE : 797
## Urban:1127 Median : 93.00 Median : 81.00
## Mean : 84.82 Mean : 69.29
## 3rd Qu.: 99.00 3rd Qu.: 97.00
## Max. :100.00 Max. :100.00
##
## missing_improved_water missing_sanitation_facilities
## Mode :logical Mode :logical
## FALSE:3299 FALSE:3196
## TRUE :32 TRUE :135
##
##
##
## mySummary <- lapply(data.frames, summary)
mySummary$`chain:1` # report just the summary of the first chain.## year region country
## Min. :1990 Africa :797 Albania : 18
## 1st Qu.:1995 Americas :613 Algeria : 18
## Median :2005 Eastern Mediterranean:373 Andorra : 18
## Mean :2002 Europe :910 Angola : 18
## 3rd Qu.:2010 South-East Asia :191 Antigua and Barbuda: 18
## Max. :2012 Western Pacific :447 Argentina : 18
## (Other) :3223
## residence_area improved_water sanitation_facilities africa
## Rural:1095 Min. : 3.00 Min. : 0.00 FALSE:2534
## Total:1109 1st Qu.: 77.00 1st Qu.: 43.00 TRUE : 797
## Urban:1127 Median : 93.00 Median : 81.00
## Mean : 84.82 Mean : 69.29
## 3rd Qu.: 99.00 3rd Qu.: 97.00
## Max. :100.00 Max. :100.00
##
## missing_improved_water missing_sanitation_facilities
## Mode :logical Mode :logical
## FALSE:3299 FALSE:3196
## TRUE :32 TRUE :135
##
##
##
## This is just a brief introduction for handling incomplete datasets. In later chapters, we will discuss more about missing data with different imputation methods and how to evaluate the complete imputed results.
1.13.1 Simulate some real multivariate data
Suppose we would like to generate a synthetic dataset: \[sim\_data=\{y, x_1, x_2, x_3, x_4, x_5, x_6, x_7, x_8, x_9, x_{10}\}.\]
Then, we can introduce a method that takes a dataset and a desired
proportion of missingness and wipes out the same proportion of the data,
i.e., introduces random patterns of missingness. Note that there are
already R functions that automate the introduction of
missingness, e.g., missForest::prodNA(), however writing
such a method from scratch is also useful.
set.seed(123)
# create MCAR missing-data generator
create.missing <- function (data, pct.mis = 10)
{
n <- nrow(data)
J <- ncol(data)
if (length(pct.mis) == 1) {
if(pct.mis>= 0 & pct.mis <=100) {
n.mis <- rep((n * (pct.mis/100)), J)
}
else {
warning("Percent missing values should be an integer between 0 and 100! Exiting"); break
}
}
else {
if (length(pct.mis) < J)
stop("The length of the missing-vector is not equal to the number of columns in the data! Exiting!")
n.mis <- n * (pct.mis/100)
}
for (i in 1:ncol(data)) {
if (n.mis[i] == 0) { # if the column has no missing values, do nothing
data[, i] <- data[, i]
}
else {
data[sample(1:n, n.mis[i], replace = FALSE), i] <- NA
# For each given column (i), sample the row indices (1:n),
# a number of indices to replace as "missing", n.mis[i], "NA",
# without replacement
}
}
return(as.data.frame(data))
}Next, let’s synthetically generate (simulate) \(1,000\) cases including all 11 features in the data (\(\{y, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10\}\)).
n <- 1000; u1 <- rbinom(n, 1, .5); v1 <- log(rnorm(n, 5, 1)); x1 <- u1*exp(v1)
u2 <- rbinom(n, 1, .5); v2 <- log(rnorm(n, 5, 1)); x2 <- u2*exp(v2)
x3 <- rbinom(n, 1, prob=0.45); x4 <- ordered(rep(seq(1, 5), n)[sample(1:n, n)])
x5 <- rep(letters[1:10], n)[sample(1:n, n)]; x6 <- trunc(runif(n, 1, 10))
x7 <- rnorm(n); x8 <- factor(rep(seq(1, 10), n)[sample(1:n, n)])
x9 <- runif(n, 0.1, .99); x10 <- rpois(n, 4)
y <- x1 + x2 + x7 + x9 + rnorm(n)
# package the simulated data as a data frame object
sim_data <- cbind.data.frame(y, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10)
# randomly create missing values
sim_data_30pct_missing <- create.missing(sim_data, pct.mis=30);
# head(sim_data_30pct_missing); summary(sim_data_30pct_missing)
# install.packages("DT")
library("DT")
library(dplyr)
df_raw <- sim_data %>% mutate_if(is.numeric, round, digits = 2)
datatable(df_raw)# install.packages("mi")
# install.packages("betareg")
library("betareg"); library("mi")
# get show the missing information matrix
mdf <- missing_data.frame(sim_data_30pct_missing)
# show(mdf)
df_mdf <- as.data.frame(mdf) %>% mutate_if(is.numeric, round, digits = 2)
datatable(df_mdf)# mdf@patterns # to get the textual missing pattern
image(mdf) # remember the visual pattern of this MCAR
# df_img <- df_mdf %>% mutate_if(is.factor, as.character) %>% replace(is.character(.), 1) %>% replace(is.na(.), 0)
# # df_img <- df_mdf %>% replace(is.character(.), 1) %>% replace(is.na(.), 0)
# df_img [1:10,1:10]
# df_img[is.character(df_img)] <- 1In the missing data plot above, missing values are illustrated as
black segments in the case-by-feature bivariate chart. The
hot colormap (17-level) represents the normalized
values of the corresponding feature-index pairs, see the mi::image()
documentation. Also, test the order,
cluster and grayscale options, e.g.,
image(mdf, x.order = T, clustered = F, grayscale =T).
The histogram plots display the distributions of:
- The observed data (in blue color),
- The imputed data (in red color), and
- The completed values (observed plus imputed, in gray color).
# Next try to impute the missing values.
# Get the Graph Parameters (plotting canvas/margins)
# set to plot the histograms for the 3 imputation chains
# mfcol=c(nr, nc). Subsequent histograms are drawn as nr-by-nc arrays on the graphics device by columns (mfcol), or rows (mfrow)
# oma
# oma=c(bottom, left, top, right) giving the size of the outer margins in lines of text
# mar=c(bottom, left, top, right) gives the number of lines of margin to be specified on the four sides of the plot.
# tcl=length of tick marks as a fraction of the height of a line of text (default=0.5)
par(mfcol=c(5, 5), oma=c(1, 1, 0, 0), mar=c(1, 1, 1, 0), tcl=-0.1, mgp=c(0, 0, 0))
# Note to get verbose output-report, parallel must be OFF: parallel=FALSE, verbose=TRUE
imputations <- mi(sim_data_30pct_missing, n.iter=5, n.chains=3, verbose=TRUE)
hist(imputations)
# Extracts several multiply imputed data.frames from "imputations" object
data.frames <- complete(imputations, 3)
# compare the 3 objects, sim_data, sim_data_30pct_missing, and imputed chain1
# datatable(sim_data, caption = htmltools::tags$caption(
# style = 'caption-side: bottom; text-align: center;','Table: Initial sim_data'))
df_miss <- sim_data_30pct_missing %>% mutate_if(is.numeric, round, digits = 2)
datatable(df_miss, caption = htmltools::tags$caption(
style = 'caption-side: bottom; text-align: center;','Table: Initial sim_data'))# datatable(sim_data_30pct_missing, caption = htmltools::tags$caption(
# style = 'caption-side: bottom; text-align: center;', 'Table: sim_data_30pct_missing'))
df_miss30pct <- sim_data_30pct_missing %>% mutate_if(is.numeric, round, digits = 2)
datatable(df_miss30pct, caption = htmltools::tags$caption(
style = 'caption-side: bottom; text-align: center;',
'Table: sim_data 30% Missing'))df_chain1 <- data.frames[[1]] %>% mutate_if(is.numeric, round, digits = 2)
datatable(df_chain1, caption = htmltools::tags$caption(
style = 'caption-side: bottom; text-align: center;', 'Table: Imputed data (chain 1)'))# Compare the summary stats for the original data (prior to introducing missing
# values) with missing data and the re-completed data following imputation
# summary(sim_data)
datatable(data.frame(t(as.matrix(unclass(summary(sim_data)))), check.names = FALSE, stringsAsFactors = FALSE), caption = htmltools::tags$caption(
style = 'caption-side: bottom; text-align: center;',
'Table: summary(sim_data)'))mySummary <- lapply(data.frames, summary)
datatable(data.frame(t(as.matrix(unclass(mySummary$`chain:1`))), check.names = FALSE, stringsAsFactors = FALSE), caption = htmltools::tags$caption(
style = 'caption-side: bottom; text-align: center;',
'Table: Imputed data: summary(chain:1)'))
Let’s check imputation convergence (details provided below).
## chain:1 chain:2 chain:3
## y 0.006 0.023 0.013
## x1 0.001 0.016 0.017
## x2 -0.020 -0.011 -0.004
## x3 1.421 1.428 1.407
## x4 2.990 2.966 2.956
## x5 5.593 5.464 5.674
## x6 -0.004 0.008 0.011
## x7 -0.007 0.019 -0.004
## x8 5.576 5.362 5.479
## x9 0.534 0.526 0.532
## x10 0.016 0.001 -0.023
## missing_y 0.300 0.300 0.300
## missing_x1 0.300 0.300 0.300
## missing_x2 0.300 0.300 0.300
## missing_x3 0.300 0.300 0.300
## missing_x4 0.300 0.300 0.300
## missing_x5 0.300 0.300 0.300
## missing_x6 0.300 0.300 0.300
## missing_x7 0.300 0.300 0.300
## missing_x8 0.300 0.300 0.300
## missing_x9 0.300 0.300 0.300
## missing_x10 0.300 0.300 0.300## mean_y mean_x1 mean_x2 mean_x3 mean_x4 mean_x5 mean_x6 mean_x7
## 1.0080355 1.2331703 1.0922942 1.3672522 0.9897120 0.9217467 0.9322015 1.1089271
## mean_x8 mean_x9 mean_x10 sd_y sd_x1 sd_x2 sd_x3 sd_x4
## 0.9176461 1.1675012 1.9580126 2.7245301 1.2391157 1.7068228 1.3820948 0.9336540
## sd_x5 sd_x6 sd_x7 sd_x8 sd_x9 sd_x10
## 1.0719535 0.9253061 0.8979584 1.0798253 1.0073414 1.0774781
## $y
## $y$is_missing
## missing
## FALSE TRUE
## 700 300
##
## $y$imputed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.28172 -0.30749 0.03040 0.04626 0.37343 1.34069
##
## $y$observed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.194077 -0.395670 0.001664 0.000000 0.347725 1.423707
##
##
## $x1
## $x1$is_missing
## missing
## FALSE TRUE
## 700 300
##
## $x1$imputed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.31408 -0.31263 0.04437 0.03690 0.35678 1.55641
##
## $x1$observed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.4805 -0.4805 -0.4805 0.0000 0.4817 1.1365
##
##
## $x2
## $x2$is_missing
## missing
## FALSE TRUE
## 700 300
##
## $x2$imputed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.72775 -0.38002 -0.05710 -0.03903 0.29907 1.49582
##
## $x2$observed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.48420 -0.48420 0.01595 0.00000 0.47667 1.14717
##
##
## $x3
## $x3$crosstab
##
## observed imputed
## 0 1215 529
## 1 885 371
##
##
## $x4
## $x4$crosstab
##
## observed imputed
## 1 408 217
## 2 429 161
## 3 438 180
## 4 411 171
## 5 414 171
##
##
## $x5
## $x5$crosstab
##
## observed imputed
## a 186 95
## b 210 88
## c 195 85
## d 231 93
## e 198 84
## f 219 96
## g 210 72
## h 219 100
## i 216 82
## j 216 105
##
##
## $x6
## $x6$is_missing
## missing
## FALSE TRUE
## 700 300
##
## $x6$imputed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.865324 -0.331768 -0.003448 0.016464 0.381098 1.644210
##
## $x6$observed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.81073 -0.41317 -0.01562 0.00000 0.38194 0.77949
##
##
## $x7
## $x7$is_missing
## missing
## FALSE TRUE
## 700 300
##
## $x7$imputed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.769935 -0.334713 0.018176 0.008462 0.352540 1.643058
##
## $x7$observed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.576156 -0.336879 0.003784 0.000000 0.369404 1.773350
##
##
## $x8
## $x8$crosstab
##
## observed imputed
## 1 213 103
## 2 222 83
## 3 186 83
## 4 210 85
## 5 228 106
## 6 219 80
## 7 210 97
## 8 207 79
## 9 198 89
## 10 207 95
##
##
## $x9
## $x9$is_missing
## missing
## FALSE TRUE
## 700 300
##
## $x9$imputed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.01106 0.34674 0.53456 0.53610 0.73667 0.99484
##
## $x9$observed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1010 0.3077 0.5100 0.5286 0.7681 0.9881
##
##
## $x10
## $x10$is_missing
## missing
## FALSE TRUE
## 700 300
##
## $x10$imputed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.758461 -0.350683 -0.026282 -0.005796 0.358105 1.829094
##
## $x10$observed
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.981263 -0.494111 -0.006959 0.000000 0.236617 1.698072
Finally, pool over the \(m = 3\) completed datasets when we fit the “model”. Pool from across the 3 chains - in order to estimate a linear regression model.
model_results <- pool(y ~ x1+x2+x3+x4+x5+x6+x7+x8+x9+x10, data=imputations, m=3)
display (model_results); summary (model_results) ## bayesglm(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 +
## x9 + x10, data = imputations, m = 3)
## coef.est coef.se
## (Intercept) -0.20 0.35
## x1 0.94 0.03
## x2 0.95 0.04
## x31 0.03 0.21
## x4.L 0.05 0.31
## x4.Q -0.04 0.12
## x4.C -0.16 0.12
## x4^4 0.13 0.20
## x5b -0.10 0.81
## x5c -0.07 0.70
## x5d -0.27 0.53
## x5e 0.53 0.55
## x5f 0.57 0.69
## x5g 0.05 1.16
## x5h 0.46 0.89
## x5i 0.50 0.83
## x5j 0.00 0.83
## x6 0.02 0.04
## x7 0.90 0.10
## x82 0.23 0.75
## x83 -0.10 0.39
## x84 0.30 0.72
## x85 -0.03 0.26
## x86 -0.41 0.57
## x87 -0.24 0.67
## x88 -0.39 0.47
## x89 -0.54 0.41
## x810 -0.56 0.85
## x9 0.99 0.21
## x10 0.09 0.07
## n = 970, k = 30
## residual deviance = 2056.5, null deviance = 15061.5 (difference = 13005.0)
## overdispersion parameter = 2.1
## residual sd is sqrt(overdispersion) = 1.46##
## Call:
## pool(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 +
## x10, data = imputations, m = 3)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.195106 0.349525 -0.558 0.58317
## x1 0.939325 0.033975 27.648 1.22e-05 ***
## x2 0.948332 0.038885 24.388 8.25e-05 ***
## x31 0.025527 0.213319 0.120 0.91198
## x4.L 0.048785 0.313158 0.156 0.88787
## x4.Q -0.041400 0.124436 -0.333 0.74255
## x4.C -0.159832 0.123819 -1.291 0.20794
## x4^4 0.129699 0.199864 0.649 0.55412
## x5b -0.101297 0.805897 -0.126 0.91032
## x5c -0.068001 0.698835 -0.097 0.93001
## x5d -0.272700 0.533810 -0.511 0.64724
## x5e 0.532085 0.551007 0.966 0.41108
## x5f 0.568015 0.687619 0.826 0.48305
## x5g 0.049315 1.159327 0.043 0.96989
## x5h 0.455083 0.894766 0.509 0.65824
## x5i 0.498684 0.834413 0.598 0.60532
## x5j -0.003877 0.832570 -0.005 0.99667
## x6 0.023434 0.035234 0.665 0.54387
## x7 0.898158 0.103438 8.683 0.00303 **
## x82 0.226402 0.746398 0.303 0.78716
## x83 -0.099904 0.388335 -0.257 0.80909
## x84 0.298364 0.722217 0.413 0.71472
## x85 -0.027635 0.262176 -0.105 0.91776
## x86 -0.413341 0.574078 -0.720 0.53038
## x87 -0.236774 0.665402 -0.356 0.75091
## x88 -0.393474 0.465155 -0.846 0.45674
## x89 -0.544436 0.407180 -1.337 0.25746
## x810 -0.561647 0.854912 -0.657 0.57393
## x9 0.993793 0.212026 4.687 2.63e-05 ***
## x10 0.085931 0.065402 1.314 0.29504
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 2.120145)
##
## Null deviance: 15061.5 on 999 degrees of freedom
## Residual deviance: 2056.5 on 970 degrees of freedom
## AIC: 3620.2
##
## Number of Fisher Scoring iterations: 7# Report the summaries of the imputations
data.frames <- complete(imputations, 3) # extract the first 3 chains
mySummary <-lapply(data.frames, summary)
datatable(data.frame(t(as.matrix(unclass(mySummary$`chain:1`))), check.names = FALSE, stringsAsFactors = FALSE), caption = htmltools::tags$caption(
style = 'caption-side: bottom; text-align: center;',
'Table: Imputed data: summary(chain:1)'))datatable(data.frame(t(as.matrix(unclass(mySummary$`chain:2`))), check.names = FALSE, stringsAsFactors = FALSE), caption = htmltools::tags$caption(
style = 'caption-side: bottom; text-align: center;',
'Table: Imputed data: summary(chain:2)'))datatable(data.frame(t(as.matrix(unclass(mySummary$`chain:3`))), check.names = FALSE, stringsAsFactors = FALSE), caption = htmltools::tags$caption(
style = 'caption-side: bottom; text-align: center;',
'Table: Imputed data: summary(chain:3)'))## Estimate Std. Error
## (Intercept) -0.195105544 0.34952452
## x1 0.939325088 0.03397494
## x2 0.948331866 0.03888537
## x31 0.025526807 0.21331910
## x4.L 0.048784767 0.31315836
## x4.Q -0.041399675 0.12443592
## x4.C -0.159832067 0.12381871
## x4^4 0.129699198 0.19986373
## x5b -0.101297131 0.80589708
## x5c -0.068001494 0.69883523
## x5d -0.272700159 0.53381011
## x5e 0.532085340 0.55100726
## x5f 0.568015316 0.68761887
## x5g 0.049314779 1.15932662
## x5h 0.455083260 0.89476618
## x5i 0.498683570 0.83441336
## x5j -0.003876698 0.83256971
## x6 0.023434393 0.03523444
## x7 0.898157946 0.10343788
## x82 0.226402068 0.74639772
## x83 -0.099904413 0.38833549
## x84 0.298363578 0.72221733
## x85 -0.027634790 0.26217633
## x86 -0.413341308 0.57407764
## x87 -0.236773958 0.66540156
## x88 -0.393474003 0.46515519
## x89 -0.544435930 0.40717954
## x810 -0.561646848 0.85491165
## x9 0.993793037 0.21202608
## x10 0.085930536 0.06540160
# plot_ly(imputations@data$`chain:1`, x=~(x1+x2), y=~density(y))
# To compare the density of observed data and imputed data --
# these should be similar (though not identical) under MAR assumptionNotes:
- In general, it is recommended to generate multiple imputation chains and then analyze the data (e.g., estimate the model coefficients, obtain inference, compute likelihoods, etc.). Pooling the analytics across all chains accounts for between-chain as well as within-chain variability, Rubin’s rule.
- When deciding on how many chains to compute, a general rule is to compute \(m\) chains if the rate of incomplete cases in the dataset is about \(m \%\), i.e., 10-chains when 10% of cases are incomplete, White et al.,2011.
- For categorical features, e.g., binary predictors like \(x_3\), the display() and summary() functions will report coefficient estimates for each (category) level, relative to the base level.
1.13.2 TBI Data Example
Next, we will see an example using the traumatic brain injury (TBI) dataset. More information about the clinical assessment scores (e.g., EGOS, GCS) is available in this publication (DOI: 10.1080/02699050701727460).
# Load the (raw) data from the table into a plain text file "08_EpiBioSData_Incomplete.csv"
TBI_Data <- read.csv("https://umich.instructure.com/files/720782/download?download_frd=1", na.strings=c("", ".", "NA")) ## 1. read in data
summary(TBI_Data)## id age sex mechanism
## Min. : 1.00 Min. :16.00 Length:46 Length:46
## 1st Qu.:12.25 1st Qu.:23.00 Class :character Class :character
## Median :23.50 Median :33.00 Mode :character Mode :character
## Mean :23.50 Mean :36.89
## 3rd Qu.:34.75 3rd Qu.:47.25
## Max. :46.00 Max. :83.00
##
## field.gcs er.gcs icu.gcs worst.gcs X6m.gose
## Min. : 3 Min. : 3.000 Min. : 0.000 Min. : 0.0 Min. :2.000
## 1st Qu.: 3 1st Qu.: 4.000 1st Qu.: 3.000 1st Qu.: 3.0 1st Qu.:3.000
## Median : 7 Median : 7.500 Median : 6.000 Median : 3.0 Median :5.000
## Mean : 8 Mean : 8.182 Mean : 6.378 Mean : 5.4 Mean :4.805
## 3rd Qu.:12 3rd Qu.:12.250 3rd Qu.: 8.000 3rd Qu.: 7.0 3rd Qu.:6.000
## Max. :15 Max. :15.000 Max. :14.000 Max. :14.0 Max. :8.000
## NA's :2 NA's :2 NA's :1 NA's :1 NA's :5
## X2013.gose skull.fx temp.injury surgery
## Min. :2.000 Min. :0.0000 Min. :0.000 Min. :0.0000
## 1st Qu.:5.000 1st Qu.:0.0000 1st Qu.:0.000 1st Qu.:0.0000
## Median :7.000 Median :1.0000 Median :1.000 Median :1.0000
## Mean :5.804 Mean :0.6087 Mean :0.587 Mean :0.6304
## 3rd Qu.:7.000 3rd Qu.:1.0000 3rd Qu.:1.000 3rd Qu.:1.0000
## Max. :8.000 Max. :1.0000 Max. :1.000 Max. :1.0000
##
## spikes.hr min.hr max.hr acute.sz
## Min. : 1.280 Min. : 0.000 Min. : 12.00 Min. :0.0000
## 1st Qu.: 5.357 1st Qu.: 0.000 1st Qu.: 35.25 1st Qu.:0.0000
## Median : 18.170 Median : 0.000 Median : 97.50 Median :0.0000
## Mean : 52.872 Mean : 3.571 Mean : 241.89 Mean :0.1739
## 3rd Qu.: 57.227 3rd Qu.: 0.000 3rd Qu.: 312.75 3rd Qu.:0.0000
## Max. :294.000 Max. :42.000 Max. :1199.00 Max. :1.0000
## NA's :18 NA's :18 NA's :18
## late.sz ever.sz
## Min. :0.0000 Min. :0.000
## 1st Qu.:0.0000 1st Qu.:0.000
## Median :1.0000 Median :1.000
## Mean :0.5652 Mean :0.587
## 3rd Qu.:1.0000 3rd Qu.:1.000
## Max. :1.0000 Max. :1.000
## # Get information matrix of the data
# 2. create an object of class "missing_data.frame" from the data.frame TBI_data
# Convert to a missing_data.frame
# library("betareg"); library("mi")
mdf <- missing_data.frame(TBI_Data) # warnings about missingness patterns## NOTE: The following pairs of variables appear to have the same missingness pattern.
## Please verify whether they are in fact logically distinct variables.
## [,1] [,2]
## [1,] "icu.gcs" "worst.gcs"
# 3. get description of the "family", "imputation_method", "size", "transformation", "type", "link", or "model" of each incomplete variable
# show(mdf)
# 4. change things: mi::change() method changes the family, imputation method,
# size, type, and so forth of a missing variable. It's called
# before calling mi to affect how the conditional expectation of each
# missing variable is modeled.
mdf <- change(mdf, y = "spikes.hr", what = "transformation", to = "identity")
# The "to" choices include "identity" = no transformation, "standardize" = standardization, "log" = natural logarithm transformation, "logshift" = log(y + a) transformation, where a is a small constant, or "sqrt" = square-root variable transformation. Changing the transformation will correspondingly change the inverse transformation.## id age sex mechanism
## Min. : 1.00 Min. :16.00 Length:46 Length:46
## 1st Qu.:12.25 1st Qu.:23.00 Class :character Class :character
## Median :23.50 Median :33.00 Mode :character Mode :character
## Mean :23.50 Mean :36.89
## 3rd Qu.:34.75 3rd Qu.:47.25
## Max. :46.00 Max. :83.00
##
## field.gcs er.gcs icu.gcs worst.gcs X6m.gose
## Min. : 3 Min. : 3.000 Min. : 0.000 Min. : 0.0 Min. :2.000
## 1st Qu.: 3 1st Qu.: 4.000 1st Qu.: 3.000 1st Qu.: 3.0 1st Qu.:3.000
## Median : 7 Median : 7.500 Median : 6.000 Median : 3.0 Median :5.000
## Mean : 8 Mean : 8.182 Mean : 6.378 Mean : 5.4 Mean :4.805
## 3rd Qu.:12 3rd Qu.:12.250 3rd Qu.: 8.000 3rd Qu.: 7.0 3rd Qu.:6.000
## Max. :15 Max. :15.000 Max. :14.000 Max. :14.0 Max. :8.000
## NA's :2 NA's :2 NA's :1 NA's :1 NA's :5
## X2013.gose skull.fx temp.injury surgery
## Min. :2.000 Min. :0.0000 Min. :0.000 Min. :0.0000
## 1st Qu.:5.000 1st Qu.:0.0000 1st Qu.:0.000 1st Qu.:0.0000
## Median :7.000 Median :1.0000 Median :1.000 Median :1.0000
## Mean :5.804 Mean :0.6087 Mean :0.587 Mean :0.6304
## 3rd Qu.:7.000 3rd Qu.:1.0000 3rd Qu.:1.000 3rd Qu.:1.0000
## Max. :8.000 Max. :1.0000 Max. :1.000 Max. :1.0000
##
## spikes.hr min.hr max.hr acute.sz
## Min. : 1.280 Min. : 0.000 Min. : 12.00 Min. :0.0000
## 1st Qu.: 5.357 1st Qu.: 0.000 1st Qu.: 35.25 1st Qu.:0.0000
## Median : 18.170 Median : 0.000 Median : 97.50 Median :0.0000
## Mean : 52.872 Mean : 3.571 Mean : 241.89 Mean :0.1739
## 3rd Qu.: 57.227 3rd Qu.: 0.000 3rd Qu.: 312.75 3rd Qu.:0.0000
## Max. :294.000 Max. :42.000 Max. :1199.00 Max. :1.0000
## NA's :18 NA's :18 NA's :18
## late.sz ever.sz
## Min. :0.0000 Min. :0.000
## 1st Qu.:0.0000 1st Qu.:0.000
## Median :1.0000 Median :1.000
## Mean :0.5652 Mean :0.587
## 3rd Qu.:1.0000 3rd Qu.:1.000
## Max. :1.0000 Max. :1.000
## image(mdf)
# 6. Perform initial imputation
imputations1 <- mi(mdf, n.iter=10, n.chains=5, verbose=TRUE)
hist(imputations1)
# 7. Extracts several multiply imputed data.frames from "imputations" object
data.frames1 <- complete(imputations1, 5)
# 8. Report a list of "summaries" for each element (imputation instance)
mySummary1 <- lapply(data.frames1, summary)
datatable(data.frame(t(as.matrix(unclass(mySummary1$`chain:1`))), check.names = FALSE, stringsAsFactors = FALSE), caption = htmltools::tags$caption(
style = 'caption-side: bottom; text-align: center;',
'Table: Imputed data: summary(chain:1)'),
extensions = 'Buttons', options = list(dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print')))datatable(data.frame(t(as.matrix(unclass(mySummary1$`chain:5`))), check.names = FALSE, stringsAsFactors = FALSE), caption = htmltools::tags$caption(
style = 'caption-side: bottom; text-align: center;',
'Table: Imputed data: summary(chain:5)'),
extensions = 'Buttons', options = list(dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print')))# 8.a. To cast the imputed numbers as integers (not necessary, but may be useful)
indx <- sapply(data.frames1[[5]], is.numeric) # get the indices of numeric columns
data.frames1[[5]][indx] <- lapply(data.frames1[[5]][indx], function(x) as.numeric(as.integer(x))) # cast each value as integer
# data.frames[[5]]$spikes.hr
# 9. Save results out
write.csv(data.frames1[[5]], "C:\\Users\\Dinov\\Desktop\\TBI_MIData.csv")
# 10. Complete Data analytics functions:
# library("mi")
#lm.mi(); glm.mi(); polr.mi(); bayesglm.mi(); bayespolr.mi(); lmer.mi(); glmer.mi()
# 10.1 Define Linear Regression for multiply imputed dataset - Also see Step (12)
##linear regression for each imputed data set - 5 regression models are fit
fit_lm1 <- glm(ever.sz ~ surgery + worst.gcs + factor(sex) + age, data.frames1$`chain:1`, family = "binomial"); summary(fit_lm1); display(fit_lm1)##
## Call:
## glm(formula = ever.sz ~ surgery + worst.gcs + factor(sex) + age,
## family = "binomial", data = data.frames1$`chain:1`)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.6719009 1.4358336 0.468 0.640
## surgery1 1.1520430 0.7218686 1.596 0.111
## worst.gcs -0.1268804 0.1053932 -1.204 0.229
## factor(sex)Male -0.4243037 0.8568291 -0.495 0.620
## age -0.0001878 0.0199172 -0.009 0.992
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 62.371 on 45 degrees of freedom
## Residual deviance: 59.050 on 41 degrees of freedom
## AIC: 69.05
##
## Number of Fisher Scoring iterations: 4## glm(formula = ever.sz ~ surgery + worst.gcs + factor(sex) + age,
## family = "binomial", data = data.frames1$`chain:1`)
## coef.est coef.se
## (Intercept) 0.67 1.44
## surgery1 1.15 0.72
## worst.gcs -0.13 0.11
## factor(sex)Male -0.42 0.86
## age 0.00 0.02
## ---
## n = 46, k = 5
## residual deviance = 59.1, null deviance = 62.4 (difference = 3.3)# Fit the appropriate model and pool the results (estimates over MI chains)
model_results <- pool(ever.sz ~ surgery + worst.gcs + factor(sex) + age, family = "binomial", data=imputations1, m=5)
display (model_results); summary (model_results) ## bayesglm(formula = ever.sz ~ surgery + worst.gcs + factor(sex) +
## age, data = imputations1, m = 5, family = "binomial")
## coef.est coef.se
## (Intercept) 0.42 1.28
## surgery1 0.92 0.65
## worst.gcs -0.09 0.10
## factor(sex)Male -0.32 0.76
## age 0.00 0.02
## n = 41, k = 5
## residual deviance = 59.5, null deviance = 62.4 (difference = 2.8)##
## Call:
## pool(formula = ever.sz ~ surgery + worst.gcs + factor(sex) +
## age, data = imputations1, m = 5, family = "binomial")
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.422459 1.275751 0.331 0.741
## surgery1 0.916494 0.654397 1.401 0.161
## worst.gcs -0.085375 0.096888 -0.881 0.378
## factor(sex)Male -0.318947 0.761881 -0.419 0.675
## age 0.002199 0.018184 0.121 0.904
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 62.371 on 45 degrees of freedom
## Residual deviance: 59.529 on 41 degrees of freedom
## AIC: 69.529
##
## Number of Fisher Scoring iterations: 6.4# Report the summaries of the imputations
data.frames <- complete(imputations1, 3) # extract the first 3 chains
mySummary2 <-lapply(data.frames1, summary)
datatable(data.frame(t(as.matrix(unclass(mySummary2$`chain:1`))), check.names = FALSE, stringsAsFactors = FALSE), caption = htmltools::tags$caption(
style = 'caption-side: bottom; text-align: center;',
'Table: Imputed data: summary(chain:1)'))# 11. Validation: we now verify whether enough iterations were conducted.
# Validation criteria demands that the mean of each completed variable should
# be similar for each of the k chains (in this case k=5).
# mipply is wrapper for sapply invoked on mi-class objects to compute the col means
round(mipply(imputations1, mean, to.matrix = TRUE), 3)## chain:1 chain:2 chain:3 chain:4 chain:5
## id 23.500 23.500 23.500 23.500 23.500
## age 0.000 0.000 0.000 0.000 0.000
## sex 1.804 1.804 1.804 1.804 1.804
## mechanism 4.261 4.261 4.261 4.261 4.261
## field.gcs -0.053 -0.060 -0.023 0.004 -0.058
## er.gcs -0.063 -0.003 0.002 -0.010 0.044
## icu.gcs -0.013 -0.049 0.012 0.035 0.016
## worst.gcs -0.016 -0.002 0.009 0.006 -0.017
## X6m.gose 0.034 0.021 0.059 -0.035 -0.014
## 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 50.128 34.455 48.623 26.833 48.059
## min.hr -0.166 -0.012 -0.131 -0.145 -0.036
## max.hr 0.050 -0.139 0.060 -0.051 -0.020
## acute.sz 1.174 1.174 1.174 1.174 1.174
## late.sz 1.565 1.565 1.565 1.565 1.565
## ever.sz 1.587 1.587 1.587 1.587 1.587
## missing_field.gcs 0.043 0.043 0.043 0.043 0.043
## missing_er.gcs 0.043 0.043 0.043 0.043 0.043
## missing_icu.gcs 0.022 0.022 0.022 0.022 0.022
## missing_worst.gcs 0.022 0.022 0.022 0.022 0.022
## missing_X6m.gose 0.109 0.109 0.109 0.109 0.109
## missing_spikes.hr 0.391 0.391 0.391 0.391 0.391
## missing_min.hr 0.391 0.391 0.391 0.391 0.391
## missing_max.hr 0.391 0.391 0.391 0.391 0.391# Rhat convergence statistics compares the variance between chains to the variance
# within chains (similar to the ANOVA F-test).
# Rhat Values ~ 1.0 indicate likely convergence,
# Rhat Values > 1.1 indicate that the chains should be run longer
# (use large number of iterations)
Rhats(imputations1, statistic = "moments") # assess the convergence of MI algorithm## mean_field.gcs mean_er.gcs mean_icu.gcs mean_worst.gcs mean_X6m.gose
## 1.142179 1.750576 3.190048 1.567101 2.268377
## mean_spikes.hr mean_min.hr mean_max.hr sd_field.gcs sd_er.gcs
## 1.876283 1.307239 1.421879 1.043464 1.218783
## sd_icu.gcs sd_worst.gcs sd_X6m.gose sd_spikes.hr sd_min.hr
## 1.721829 1.075915 1.442252 1.427767 1.211636
## sd_max.hr
## 1.818859# When convergence is unstable, we can continue the iterations for all chains, e.g.
imputations1 <- mi(imputations1, n.iter=20) # add additional 20 iterations
# To plot the produced mi results, for all missing_variables we can generate
# a histogram of the observed, imputed, and completed data.
# We can compare of the completed data to the fitted values implied by the model
# for the completed data, by plotting binned residuals.
# hist function works similarly as plot.
# image function gives a sense of the missingness patterns in the data
plot(imputations1); hist(imputations1); image(imputations1)
mySummary3 <-lapply(data.frames1, summary)
datatable(data.frame(t(as.matrix(unclass(mySummary3$`chain:1`))), check.names = FALSE, stringsAsFactors = FALSE), caption = htmltools::tags$caption(
style = 'caption-side: bottom; text-align: center;',
'Table: Imputed data: summary(chain:1)'))# 12. Finally, pool over the m = 5 imputed datasets when we fit the "model"
# Pool from across the 4 chains - in order to estimate a linear regression model
# and impact of various predictors
model_results <- pool(ever.sz ~ surgery + worst.gcs + factor(sex) + age, data = imputations1, m = 5 ); display (model_results); summary (model_results)## bayesglm(formula = ever.sz ~ surgery + worst.gcs + factor(sex) +
## age, data = imputations1, m = 5)
## coef.est coef.se
## (Intercept) 0.55 1.30
## surgery1 0.98 0.66
## worst.gcs -0.10 0.09
## factor(sex)Male -0.36 0.77
## age 0.00 0.02
## n = 41, k = 5
## residual deviance = 59.0, null deviance = 62.4 (difference = 3.3)##
## Call:
## pool(formula = ever.sz ~ surgery + worst.gcs + factor(sex) +
## age, data = imputations1, m = 5)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.5540667 1.2991773 0.426 0.670
## surgery1 0.9767509 0.6627950 1.474 0.141
## worst.gcs -0.0999804 0.0948404 -1.054 0.292
## factor(sex)Male -0.3560410 0.7681894 -0.463 0.643
## age 0.0002953 0.0186023 0.016 0.987
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 62.371 on 45 degrees of freedom
## Residual deviance: 59.044 on 41 degrees of freedom
## AIC: 69.044
##
## Number of Fisher Scoring iterations: 6.6## Estimate Std. Error
## (Intercept) 0.5540666851 1.29917731
## surgery1 0.9767508814 0.66279496
## worst.gcs -0.0999803749 0.09484039
## factor(sex)Male -0.3560410168 0.76818943
## age 0.0002952953 0.01860234
1.13.3 Imputation via Expectation-Maximization
Below we present the theory and practice of one specific statistical computing strategy for imputing incomplete datasets.
1.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 devices that can only measure values above \(0.5\), anything below that is missing.
1.13.3.2 General Idea of the EM algorithm
Expectation-Maximization (EM) is an iterative process involving two steps - expectation and maximization, which are applied in tandem. EM can be employed to find parameter estimates using maximum likelihood and is specifically useful when the equations determining the relations of the data-parameters cannot be directly solved. For example, a Gaussian mixture modeling assumes that each data point (\(X\)) has a corresponding latent (unobserved) variable or a missing value (\(Y\)), which may be specified as a mixture of coefficients determining the affinity of the data as a linear combination of Gaussian kernels, determined by a set of parameters (\(\theta\)), e.g., means and variance-covariances. Thus, EM estimation relies on:
- An observed data set \(X\),
- A set of missing (or latent) values \(Y\),
- A parameter \(\theta\), which may be a vector of parameters,
- A likelihood function \(L(\theta | X, Y) =p(X,Y |\theta)\), and
- The maximum likelihood estimate (MLE) of the unknown parameter(s) \(\theta\) that is computed using the marginal likelihood of the observed data:
\[L(\theta | X) = p(X |\theta) =\int { p(X, Y |\theta)dY}.\]
Most of the time, this equation may not be directly solved, e.g., when \(Y\) is missing.
- Expectation step (E step): computes the expected value of the log likelihood function, with respect to the conditional distribution of \(Y\) given \(X\) using the parameter estimates at the previous iteration (or at the position of initialization, for the first iteration), \(\theta_t\): \[Q ( \theta | \theta^{(t)} ) = E_{Y | X, \theta^{(t)} }[ log \left ( L(\theta | X , Y ) \right ];\]
- Maximization step (M step): Determine the parameter, \(\theta\), that maximizes the expectation above, \[\theta^{(t+1)}=\arg\max_{\theta}Q(\theta|\theta^{(t)}).\]
This SOCR EM Activity shows the practical aspects of applying the EM algorithm. Also, in DSPA Chapter 3 we will illustrate the EM method for fitting single distribution models or (linear) mixtures of distributions to data that may represent a blend of heterogeneous observations from multiple different processes.
1.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 there are complications in calculating the MLE. It is applicable for imputing incomplete MAR data, where the missing data mechanism can be ignored and separate parameters may be estimated for each missing feature.
Complete Data: \[Z = \left(\begin{array}{cc} X \\ Y \end{array}\right), ZZ^T = \left(\begin{array}{cc} XX^T & XY^T \\ YX^T & YY^T \end{array}\right),\] where \(X\) is the observed data and \(Y\) is the missing data.
- E-step: (Expectation) Get the expectations of \(Y\) and \(YY^T\) based on observed data, \(X\).
- M-step: (Maximization) Maximize the conditional expectation in E-step to estimate the parameters.
Details: If \(o=obs\) and \(m=mis\) stand for observed and missing, the mean vector, \((\mu_{obs}, \mu_{mis})^T\), and the variance-covariance matrix, \(\Sigma^{(t)} = \left(\begin{array}{cc} \Sigma_{oo} & \Sigma_{om} \\ \Sigma_{mo} & \Sigma_{mm} \end{array}\right)\), are represented by:
\[\mu^{(t)} = \left(\begin{array}{cc} \mu_{obs} \\ \mu_{mis} \end{array}\right),\;\;\;\;\; \Sigma^{(t)} = \left(\begin{array}{cc} \Sigma_{oo} & \Sigma_{om} \\ \Sigma_{mo} & \Sigma_{mm} \end{array}\right)\] E-step:
\[E(Z | X) = \left(\begin{array}{cc} X \\ E(Y|X) \end{array}\right),\;\;\;\;\; E(ZZ^T|X) = \left(\begin{array}{cc} XX^T & XE(Y|X)^T \\ E(Y|X)X^T & E(YY^T|X) \end{array}\right).\]
\[E(Y | X) = \mu_{mis} + \Sigma_{mo}\Sigma_{oo}^{-1}(X - \mu_{obs}).\] \[E(YY^T|X) = (\Sigma_{mm}-\Sigma_{mo}\Sigma_{oo}^{-1}\Sigma_{om})+E(Y|X)E(Y|X)^T.\]
M-step: \[\mu^{(t+1)} = \frac{1}{n}\sum_{i=1}^nE(Z|X).\] \[\Sigma^{(t+1)} = \frac{1}{n}\sum_{i=1}^nE(ZZ^T|X) - \mu^{(t+1)}{\mu^{(t+1)}}^T.\]
1.13.3.4 A simple manual implementation of EM-based imputation
# install.packages(c("gridExtra", "MASS"))
library(ggplot2)
library(gridExtra)
library(MASS)
library(knitr)# 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
# install.packages("e1071")
# 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] 466# 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")
df_mdf <- sim_data.df %>% mutate_if(is.numeric, round, digits = 2)
datatable(df_mdf, caption = htmltools::tags$caption(
style = 'caption-side: bottom; text-align: center;',
'Table: Simulated Data (sim_data.df)'),
extensions = 'Buttons', options = list(dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print')))# datatable(sim_data.df, caption = htmltools::tags$caption(
# style = 'caption-side: bottom; text-align: center;',
# 'Table: Simulated Data (sim_data.df)'),
# extensions = 'Buttons', options = list(dom = 'Bfrtip',
# buttons = c('copy', 'csv', 'excel', 'pdf', 'print')))
# Define the EM imputation method
EM_algorithm <- function(x, tol = 0.001) {
# identify the missing data entries (Boolean indices)
missvals <- is.na(x)
# instantiate 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 inverse-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)
df_mdf <- sim_data.imputed %>% mutate_if(is.numeric, round, digits = 2)
datatable(df_mdf, caption = htmltools::tags$caption(
style = 'caption-side: bottom; text-align: center;',
'Table: EM-Imputed Simulated Data'),
extensions = 'Buttons', options = list(dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print')))1.13.3.5 Plotting the complete and imputed data
Smaller points colored in black represent observed data, and the circle-shapes colored in magenta denote the imputed data.
plot.me <- function(index1, index2){
plot.imputed <- sim_data.imputed[row.names(
subset(sim_data.df, is.na(sim_data.df[, index1]) | is.na(sim_data.df[, index2]))), ]
p = ggplot(sim_data.imputed, aes_string( paste0("X",index1) , paste0("X",index2 ))) +
geom_point(alpha = 0.5, size = 0.7)+theme_bw() +
stat_ellipse(type = "norm", color = "#000099", alpha=0.5) +
geom_point(data = plot.imputed, aes_string( paste0("X",index1) , paste0("X",(index2))),size = 1.5, color = "Magenta", alpha = 0.8)
}
gridExtra::grid.arrange( plot.me(1,2), plot.me(5,6), plot.me(13,20), plot.me(18,19), nrow = 2)
index1=1; index2=5
plot.imputed <- sim_data.imputed[row.names(
subset(sim_data.df, is.na(sim_data.df[, index1]) | is.na(sim_data.df[, index2]))), ]
p = ggplot(sim_data.imputed, aes_string( paste0("X",index1) , paste0("X",index2 ))) +
geom_point(alpha = 0.5, size = 0.7)+theme_bw() +
stat_ellipse(type = "norm", color = "#000099", alpha=0.5) +
geom_point(data = plot.imputed, aes_string( paste0("X",index1) , paste0("X",(index2))),size = 1.5, color = "Magenta", alpha = 0.8)
plot_ly(sim_data.imputed, x = ~X1, y = ~X5, type = "scatter",
mode = "markers") %>%
layout(title='Scatterplot: Improved Water Quality vs. Sanitation Facilities',
xaxis = list (title = 'Water Quality'), yaxis = list (title = 'Sanitation'))1.13.3.6 Validation of
EM-imputation using the R Package Amelia
1.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.
## [1] 200 20## -- Imputation 1 --
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 21
##
## -- Imputation 2 --
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
##
## -- Imputation 3 --
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
##
## -- Imputation 4 --
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
##
## -- Imputation 5 --
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 21 22 23##
## Amelia output with 5 imputed datasets.
## Return code: 1
## Message: Normal EM convergence.
##
## Chain Lengths:
## --------------
## Imputation 1: 21
## Imputation 2: 17
## Imputation 3: 17
## Imputation 4: 15
## Imputation 5: 23- 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)
1.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 the \(4,000\) that we synthetically generated.
# plot.ii3 <- function(index){
# imputed <- sim_data.imputed[is.na(sim_data.df[, index]) , index]
# imputed.amelia <- amelia.imputed.5[is.na(sim_data.df[, index]) , index]
# observed <- sim_data.df[!is.na(sim_data.df[, index]) , index]
# imputed.df <- data.frame(x = c(observed,imputed,imputed.amelia), category = c( rep("obs",length(observed)),rep("simpleImplement",length(imputed)) ,rep("amelia",length(imputed.amelia)) ) )
# p = ggplot(imputed.df, aes(x=x, y =..density..)) +
# geom_density(aes(fill = category),alpha=0.3)+
# theme_bw()
# return(p)
# }
# grid.arrange( plot.ii3(1),plot.ii3(2),plot.ii3(3),plot.ii3(4),plot.ii3(5),
# plot.ii3(6),plot.ii3(7),plot.ii3(8),plot.ii3(9),plot.ii3(10),
# nrow = 5)
library(tidyr)
myPlotly <- function(index){
imputed <- sim_data.imputed[is.na(sim_data.df[, index]) , index]
imputed.amelia <- amelia.imputed.5[is.na(sim_data.df[, index]) , index]
observed <- sim_data.df[!is.na(sim_data.df[, index]) , index]
imputed.df <- data.frame(x = c(observed,imputed,imputed.amelia),
category = c( rep("obs",length(observed)),rep("simpleImplement",length(imputed)),
rep("amelia",length(imputed.amelia)) ) )
df_long <- as.data.frame(cbind(index=c(1:length(imputed.df$x)),
category=imputed.df$category, x=imputed.df$x))
df_wide <- spread(df_long, category, x)
p = plot_ly() %>%
add_lines(x = ~density(as.numeric(df_wide$simpleImplement), na.rm = T)$x,
y= ~density(as.numeric(df_wide$simpleImplement), na.rm = T)$y, name = "EM", mode = 'lines') %>%
add_lines(x = density(as.numeric(df_wide$amelia), na.rm = T)$x,
y= density(as.numeric(df_wide$amelia), na.rm = T)$y, name = "Amelia", mode = 'lines') %>%
add_lines(x = ~density(as.numeric(df_wide$obs), na.rm = T)$x,
y= ~density(as.numeric(df_wide$obs), na.rm = T)$y, name = "Observed", mode = 'lines') %>%
layout(title=sprintf("Distributions: Feature X.%d", index),
xaxis = list(title = 'Measurements'),
yaxis = list(title ="Densities"),
legend = list(title="Distributions", orientation = 'h'))
return(p)
}
# Plot a few features
myPlotly(5)1.14 Parsing web pages and visualizing tabular HTML data
In this section, we will utilize the Earthquakes dataset on SOCR website. It records information about earthquakes that happened between 1969 and 2007 with magnitudes larger than 5 on the Richter scale. Here is how we parse the data on the source webpage and ingest the information into R:
# install.packages("xml2")
library("XML"); library("xml2")
library("rvest")
wiki_url <- read_html("https://wiki.socr.umich.edu/index.php/SOCR_Data_Dinov_021708_Earthquakes")
html_nodes(wiki_url, "#content")## {xml_nodeset (1)}
## [1] <div id="content" class="mw-body" role="main">\n\t\t\t<a id="top"></a>\n\ ...In this dataset, Magt(magnitude type) may be used as a
grouping variable. We will draw a “Longitude vs Latitude” line plot from
this dataset. The function we are using is called ggplot()
under ggplot2. The input type for this function is mostly
data frame. aes() specifies axes.
# library(ggplot2)
# plot4<-ggplot(earthquake, aes(Longitude, Latitude, group=Magt, color=Magt))+
# geom_point(data=earthquake, size=4, mapping=aes(x=Longitude, y=Latitude, shape=Magt))
# plot4 # or plint(plot4)
# https://plotly-r.com/working-with-symbols.html
glyphication <- function (name) {
glyph= vector()
for (i in 1:length(name)){
glyph[i]="triangle-up"
if (name[i]=="Md") { glyph[i]="diamond-open" }
else if (name[i]=="ML") { glyph[i]="circle-open" }
else if (name[i]=="Mw") { glyph[i]="square-open" }
else if (name[i]=="Mx") { glyph[i]="x-open" }
}
return(glyph)
}
earthquake$glyph <- glyphication(earthquake$Magt)
plot_ly(earthquake) %>%
add_markers(x = ~Longitude, y = ~Latitude, type = "scatter", color = ~Magt,
mode = "markers", marker = list(size = ~Depth, color = ~Magt, symbol = ~glyph,
line = list(color = "black",width = 2))) %>%
layout(title="California Earthquakes (1969 - 2007)")The most important line of code has 2 parts. The first part
ggplot(earthquake, aes(Longitude, Latitude, group=Magt, color=Magt))
specifies the setting of the plot: dataset, group and color. The second
part specifies we are going to draw lines between data points. In later
chapters we will frequently use the package ggplot2 and the
structure under this great package is always
function1+function2.
We can visualize the distribution for different variables using
density plots. The following script plots the distribution for Latitude
among different Magnitude types, also using the ggplot()
function combined with geom_density().
We can also compute and display 2D Kernel Density and 3D Surface Plots. Plotting 2D Kernel Density and 3D Surface plots is very important and useful in multivariate exploratory data analytic.
We will use the plot_ly() function under the
plotly package, which takes value from a data frame.
To create a surface plot, we use two vectors: x and y with length m and n respectively. We also need a matrix: z of size \(m\times n\). This z matrix is created from matrix multiplication between x and y.
The kde2d() function is needed for 2D kernel density
estimation.
Here z is an estimate of the kernel density function.
Then we apply plot_ly to the list
kernal_density via the with() function.
Note that we used the option "surface", however you can
experiment with the type option.
Alternatively, one can plot 1D, 2D or 3D plots:
df3D <- data.frame(x=earthquake$Longitude, y=earthquake$Latitude, z=earthquake$Mag)
# Convert he Long (X, Y, Z) Earthquake format data into a Matrix Format
# install.packages("Matrix")
library("Matrix")
matrix_EarthQuakes <- with(df3D, sparseMatrix(i = as.numeric(180-x), j=as.numeric(y), x=z, use.last.ij=T, dimnames=list(levels(x), levels(y))))
dim(matrix_EarthQuakes)## [1] 307 44# colnames(matrix_EarthQuakes) <- seq(from=earthquake$Longitude[1],
# to=earthquake$Longitude[length(earthquake$Longitude)],
# length.out=dim(matrix_EarthQuakes)[2])
# rownames(matrix_EarthQuakes) <- seq(from=earthquake$Latitude[1],
# to=earthquake$Latitude[length(earthquake$Latitude)],
# length.out=dim(matrix_EarthQuakes)[1])
# View(as.matrix(matrix_EarthQuakes))
# view matrix is 2D heatmap:
library("ggplot2"); library("gplots")
# heatmap.2( as.matrix(matrix_EarthQuakes[280:307, 30:44]), Rowv=FALSE, Colv=FALSE, dendrogram='none', cellnote=as.matrix(matrix_EarthQuakes[280:307, 30:44]), notecol="black", trace='none', key=FALSE, lwid = c(.01, .99), lhei = c(.01, .99), margins = c(5, 15 ))
plot_ly(z = ~as.matrix(matrix_EarthQuakes[280:307, 30:44]), type = "heatmap") %>% hide_colorbar()# plot_ly(x=~colnames(matrix_EarthQuakes[280:307, 30:44]),
# y=~rownames(matrix_EarthQuakes[280:307, 30:44]),
# z = ~as.matrix(matrix_EarthQuakes[280:307, 30:44]), type = "heatmap") %>%
# layout(title="California Earthquakes Heatmap",
# xaxis=list(title="Longitude"), yaxis=list(title="Latitude")) %>%
# hide_colorbar()
# Long -180<x<-170, Lat: 30<y<45, Z: 5<Mag<8
matrix_EarthQuakes <- with(df3D, sparseMatrix(i = as.numeric(180+x), j=as.numeric(y), x=z, use.last.ij=TRUE, dimnames=list(levels(x), levels(y))))
mat1 <- as.matrix(matrix_EarthQuakes)
plot_ly(z = ~mat1, type = "surface")1.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 group is more than an order of
magnitude larger than the size of another group has the
potential for bias.
1.15.1 Example 1: Parkinson’s Diseases Study
This Parkinson’s diseases case-study involves neuroimaging, genetics, clinical, and phenotypic data for over 600 volunteers produced multivariate data for 3 cohorts – HC=Healthy Controls(166) , PD=Parkinson’s (434), SWEDD= subjects without evidence for dopaminergic deficit (61).
# update packages
# update.packages()
# load the data: 06_PPMI_ClassificationValidationData_Short.csv
ppmi_data <-read.csv("https://umich.instructure.com/files/330400/download?download_frd=1", header=TRUE)
table(ppmi_data$ResearchGroup)
# binarize the Dx classes
ppmi_data$ResearchGroup <- ifelse(ppmi_data$ResearchGroup == "Control", "Control", "Patient")
attach(ppmi_data)
head(ppmi_data)
# Model-free analysis, classification
# install.packages("crossval")
# install.packages("ada")
# library("crossval")
library(crossval)
library(ada)
#set up adaboosting prediction function
# Define a new classification result-reporting function
my.ada <- function (train.x, train.y, test.x, test.y, negative, formula){
ada.fit <- ada(train.x, train.y)
predict.y <- predict(ada.fit, test.x)
#count TP, FP, TN, FN, Accuracy, etc.
out <- confusionMatrix(test.y, predict.y, negative = negative)
# negative is the label of a negative "null" sample (default: "control").
return (out)
}
# balance cases
# SMOTE: Synthetic Minority Oversampling Technique to handle class imbalance in binary classification.
set.seed(1000)
# https://cran.r-project.org/src/contrib/Archive/unbalanced/
# install.packages('mlr', 'FNN', 'RANN', 'unbalanced') to deal with unbalanced group data
library(unbalanced)
ppmi_data$PD <- ifelse(ppmi_data$ResearchGroup=="Control", 1, 0)
uniqueID <- unique(ppmi_data$FID_IID)
ppmi_data <- ppmi_data[ppmi_data$VisitID==1, ]
ppmi_data$PD <- factor(ppmi_data$PD)
colnames(ppmi_data)
# ppmi_data.1<-ppmi_data[, c(3:281, 284, 287, 336:340, 341)]
n <- ncol(ppmi_data)
output.1 <- ppmi_data$PD
# remove Default Real Clinical subject classifications!
ppmi_data$PD <- ifelse(ppmi_data$ResearchGroup=="Control", 1, 0)
input <- ppmi_data[ , -which(names(ppmi_data) %in% c("ResearchGroup", "PD", "X", "FID_IID"))]
# output <- as.matrix(ppmi_data[ , which(names(ppmi_data) %in% {"PD"})])
output <- as.factor(ppmi_data$PD)
c(dim(input), length(output))
#balance the dataset
data.1<-ubBalance(X= input, Y=output, type="ubSMOTE", percOver=300, percUnder=150, verbose=TRUE)
# percOver = A number that drives the decision of how many extra cases from the minority class are generated (known as over-sampling).
# k = A number indicating the number of nearest neighbors that are used to generate the new examples of the minority class.
# percUnder = A number that drives the decision of how many extra cases from the majority classes are selected for each case generated from the minority class (known as under-sampling)
balancedData<-cbind(data.1$X, data.1$Y)
table(data.1$Y)
nrow(data.1$X); ncol(data.1$X)
nrow(balancedData); ncol(balancedData)
nrow(input); ncol(input)
colnames(balancedData) <- c(colnames(input), "PD")
# check visually for differences between the distributions of the raw (input) and rebalanced data (for only one variable, in this case)
QQ <- qqplot(input[, 5], balancedData [, 5], plot.it=F)
plot_ly(x=~QQ$x, y = ~QQ$y, type="scatter", mode="markers", showlegend=F) %>%
add_lines(x=c(0,0.8), y=c(0,0.8), showlegend=F) %>%
layout(title="QQ-Plot Original vs. Rebalanced Data", xaxis=list(title="original data"),
yaxis=list(title="Rebalanced data"))
###Check balance
## Wilcoxon test
alpha.0.05 <- 0.05
test.results.bin <- NULL # binarized/dichotomized p-values
test.results.raw <- NULL # raw p-values
for (i in 1:(ncol(balancedData)-1)) {
test.results.raw [i] <- wilcox.test(input[, i], balancedData [, i])$p.value
test.results.bin [i] <- ifelse(test.results.raw [i] > alpha.0.05, 1, 0)
print(c("i=", i, "Wilcoxon-test=", test.results.raw [i]))
}
print(c("Wilcoxon test results: ", test.results.bin))
test.results.corr <- stats::p.adjust(test.results.raw, method = "fdr", n = length(test.results.raw))
# where methods are "holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none")
# plot(test.results.raw, test.results.corr)
# zeros (0) are significant independent between-group T-test differences, ones (1) are insignificant
plot_ly(x=~test.results.raw, y = ~test.results.corr, type="scatter", mode="markers", showlegend=F) %>%
add_lines(x=c(0,1), y=c(0,1), showlegend=F) %>%
layout(title="Wilcoxon test results - Original vs. Rebalanced Data", xaxis=list(title="Original"),
yaxis=list(title="Rebalanced"))
# Check the Differences between the rate of significance between the raw and FDR-corrected p-values
test.results.bin <- ifelse(test.results.raw > alpha.0.05, 1, 0)
table(test.results.bin)
test.results.corr.bin <- ifelse(test.results.corr > alpha.0.05, 1, 0)
table(test.results.corr.bin)Notes
- SMOTE oversampling of the minority cohort is via generation of synthetic minority samples within the neighborhoods of observed observations. Thus, new minority instances blend observations in the same class and create clusters around each observed minority element.
- The
percOverparameter (perc.over/100) represents the number of new instances generated for each rare instance in the minority sample, when \(perc.over < 100\), a single instance is generated. For example,percOver=300andpercOver=30would triple (300/100) and leave unchanged (30/100) the size of the minority sample, respectively. - The \(k\) parameter represents the number of neighbors to consider as the aggregate pool that the new examples are generated.
- The
percUnder(perc.under/100) represents the number of “normal” (majority class) instances that are randomly selected for each smoted (synthetically generated) observation. For instance,percUnder=300orpercUnder=30would downsample the majority sample by choosing one-out-of-each-three or all of the majority sample points, respectively.
2 Exploratory Data Analytics (EDA)
In this section, we will see a broad range of simulations and hands-on activities to highlight some of the basic data visualization techniques using R. A brief discussion of alternative visualization methods is followed by demonstrations of histograms, density, pie, jitter, bar, line and scatter plots, as well as strategies for displaying trees and graphs and 3D surface plots. Many of these are also used throughout the textbook in the context of addressing the graphical needs of specific case-studies.
It is practically impossible to cover all options of every different
visualization routine. Readers are encouraged to experiment with each
visualization type, change input data and parameters, explore the
function documentation using R-help (e.g., ?plot), and
search for new R visualization packages and new functionality, which are
continuously being developed.
2.1 General Questions Driving Visualization
- What exploratory visualization techniques are available to visually interrogate my specific data?
- How to examine paired associations and correlations in a multivariate dataset?
2.2 Classification of visualization methods
Scientific data-driven or simulation-driven visualization methods are hard to classify. The following list of criteria can be used for classification:
- Data Type: structured/unstructured, small/large, complete/incomplete, time/space, ASCII/binary, Euclidean/non-Euclidean, etc.
- Task type: Task type is one of the aspects considered in classification of visualization techniques, which provides means of interaction between the researcher, the data and the display software/platform
- Scalability: Visualization techniques are subject
to some limitations, such as the amount of data that a particular
technique can exhibit
- Dimensionality: Visualization techniques can also be classified according to the number of attributes
- Positioning and Attributes: the distribution of attributes on the chart may affect the interpretation of the display representation, e.g., correlation analysis, where the relative distance among the plotted attributes is relevant for observation
- Investigative Need: the specific scientific question or exploratory interest may also determine the type of visualization:
- Examining the composition of the data
- Exploring the distribution of the data
- Contrasting or comparing several data elements, relations, association
- Unsupervised exploratory data mining.
Also, we have the following table for common data visualization methods according to task types:
We chose to introduce common data visualization methods according to this classification criterion, albeit this is not a unique or even broadly agreed upon ontological characterization of exploratory data visualization.
2.3 Composition
In this section, we will see composition plots for different types of variables and data structures.
2.3.1 Histograms and density plots
One of the first few graphs we learned in high school would be
Histogram. In R, the functions hist() or
plot_ly() represent two methods that can be applied to a
vector of values for plotting histograms. The famous 19-th century
statistician Karl
Pearson introduced histograms as graphical representations of the
distribution of a sample of numeric data. The histogram plot uses the
data to infer and display the probability distribution of the underlying
population that the data is sampled from. Histograms are constructed by
selecting a certain number of bins covering the range of values of the
observed process. Typically, the number of bins for a data array of size
\(N\) should be equal to \(\sqrt{N}\). These bins form a partition
(disjoint and covering sets) of the range. Finally, we compute the
relative frequency representing the number of observations that fall
within each bin interval. The histogram just plots a piecewise
step-function defined over the union of the bin interfaces whose height
equals the observed relative frequencies.
# Here `freq=T` shows the frequency for each *x* value and `breaks` controls for the number of bars in our histogram.
# mu <- 15; sd <- 3.7
# set.seed(1234)
# x<-rnorm(100, mean = mu, sd=sd)
# hist(x, freq=F, breaks = 10)
# lines(density(x), lwd=2, col="blue")
# t <- seq(mu-3*sd, mu+3*sd, by=0.01)
# lines(t, dnorm(t,mu,sd), col="magenta") # add the theoretical density line
library(plotly)
N <- 10000
mu <- 15; sd <- 3.7
set.seed(1234)
x <- rnorm(N, mean = mu, sd=sd)
fit <- density(x)
z<-seq(mu-4*sd, mu+4*sd, 0.1) # points from -4 to 4 in 0.1 steps
q<-seq(0.001, 0.999, 0.001) # probability quantile values from 0.1% to 99.9% in 0.1% steps
normDensity <- dnorm(z, mean=15, sd= 3.7)
plot_ly(x = x, type = "histogram", name = "Data Histogram", histnorm = "probability") %>%
add_trace(x = fit$x, y = fit$y, type = "scatter", mode = "lines", opacity=0.1,
fill = "tozeroy", yaxis = "y2", name = "Density (rnorm(100, 15, 3.7))") %>%
add_trace(x = z, y = normDensity, type = "scatter", mode = "lines", opacity=0.1,
fill = "tozeroy", yaxis = "y2", name = "Normal(15, 3.7)") %>%
layout(title='Data Histogram, Density Estimate & Theoretical Model Distribution',
yaxis2 = list(overlaying = "y", side = "right"),
legend = list(orientation = 'h'))The shape of the last histogram we draw is very close to a Normal
distribution (because we sampled from this distribution by
rnorm). Note the superposition of the corresponding Normal
density curve.
2.3.2 Pie Chart
We are all very familiar with pie charts that show us the components of a big “cake”. Although pie charts provide effective simple visualization in certain situations, it may also be difficult to compare segments within a pie chart or across different pie charts. Other plots like bar chart, box or dot plots may be attractive alternatives.
We will use the Letter Frequency Data on SOCR website to illustrate the use of pie charts.
library(rvest)
wiki_url <- read_html("https://wiki.socr.umich.edu/index.php/SOCR_LetterFrequencyData")
html_nodes(wiki_url, "#content")## {xml_nodeset (1)}
## [1] <div id="content" class="mw-body" role="main">\n\t\t\t<a id="top"></a>\n\ ...## Letter English French German
## Length:27 Min. :0.00000 Min. :0.00000 Min. :0.00000
## Class :character 1st Qu.:0.01000 1st Qu.:0.01000 1st Qu.:0.01000
## Mode :character Median :0.02000 Median :0.03000 Median :0.03000
## Mean :0.03667 Mean :0.03704 Mean :0.03741
## 3rd Qu.:0.06000 3rd Qu.:0.06500 3rd Qu.:0.05500
## Max. :0.13000 Max. :0.15000 Max. :0.17000
## Spanish Portuguese Esperanto Italian
## Min. :0.00000 Min. :0.00000 Min. :0.00000 Min. :0.00000
## 1st Qu.:0.01000 1st Qu.:0.00500 1st Qu.:0.01000 1st Qu.:0.00500
## Median :0.03000 Median :0.03000 Median :0.03000 Median :0.03000
## Mean :0.03815 Mean :0.03778 Mean :0.03704 Mean :0.03815
## 3rd Qu.:0.06000 3rd Qu.:0.05000 3rd Qu.:0.06000 3rd Qu.:0.06000
## Max. :0.14000 Max. :0.15000 Max. :0.12000 Max. :0.12000
## Turkish Swedish Polish Toki_Pona
## Min. :0.00000 Min. :0.00000 Min. :0.00000 Min. :0.00000
## 1st Qu.:0.01000 1st Qu.:0.01000 1st Qu.:0.01500 1st Qu.:0.00000
## Median :0.03000 Median :0.03000 Median :0.03000 Median :0.03000
## Mean :0.03667 Mean :0.03704 Mean :0.03704 Mean :0.03704
## 3rd Qu.:0.05500 3rd Qu.:0.05500 3rd Qu.:0.04500 3rd Qu.:0.05000
## Max. :0.12000 Max. :0.10000 Max. :0.20000 Max. :0.17000
## Dutch Avgerage
## Min. :0.00000 Min. :0.00000
## 1st Qu.:0.01000 1st Qu.:0.01000
## Median :0.02000 Median :0.03000
## Mean :0.03704 Mean :0.03741
## 3rd Qu.:0.06000 3rd Qu.:0.06000
## Max. :0.19000 Max. :0.12000We can try to plot the frequency proportion of the 26 English letters using pie and donut charts.
# The left hand side plot is the one without reference table and the right one has the table made by function `legend`.
# par(mfrow=c(1, 2))
# pie(letter$English[1:10], labels=letter$Letter[1:10], col=rainbow(10, start=0.1, end=0.8), clockwise=TRUE, main="First 10 Letters Pie Chart")
# pie(letter$English[1:10], labels=letter$Letter[1:10], col=rainbow(10, start=0.1, end=0.8), clockwise=TRUE, main="First 10 Letters Pie Chart")
# legend("topleft", legend=letter$Letter[1:10], cex=1.3, bty="n", pch=15, pt.cex=1.8, col=rainbow(10, start=0.1, end=0.8), ncol=1)
plot_ly(letter, labels = ~Letter, values = ~English, type = 'pie', name="English",
textposition = 'inside', textinfo = 'label+percent', showlegend = FALSE,
domain = list(row = 0, column = 0)) %>%
add_pie(labels = ~Letter, values = ~Spanish, name = "Spanish",
textposition = 'inside', textinfo = 'label+percent', showlegend = FALSE,
domain = list(row = 0, column = 1)) %>%
add_pie(labels = ~Letter, values = ~Swedish, name = "Swedish",
textposition = 'inside', textinfo = 'label+percent', showlegend = FALSE,
domain = list(row = 1, column = 0)) %>%
add_pie(labels = ~Letter, values = ~Polish, name = "Polish",
textposition = 'inside', textinfo = 'label+percent', showlegend = FALSE,
domain = list(row = 1, column = 1)) %>%
add_annotations(x=0.01, y=0.99,text = "English",showarrow = F, ax = 20, ay = -40) %>%
add_annotations(x=0.58, y=0.99,text = "Spanish",showarrow = F, ax = 20, ay = -40) %>%
add_annotations(x=0.01, y=0.01,text = "Swedish",showarrow = F, ax = 20, ay = -40) %>%
add_annotations(x=0.58, y=0.01,text = "Polish",showarrow = F, ax = 20, ay = -40) %>%
layout(title = 'Pie Charts of English, Spanish, Swedish & Polish Letters',
grid=list(rows=2, columns=2),
xaxis = list(showgrid = FALSE, zeroline = FALSE, showticklabels = FALSE),
yaxis = list(showgrid = FALSE, zeroline = FALSE, showticklabels = FALSE))plot_ly(letter, labels = ~Letter, values = ~German, type = 'pie', name="German",
textposition = 'inside', textinfo = 'label+percent', showlegend = FALSE,
domain = list(row = 0, column = 0), hole = 0.5) %>%
add_pie(labels = ~Letter, values = ~Italian, name = "Italian",
textposition = 'inside', textinfo = 'label+percent', showlegend = FALSE,
domain = list(row = 0, column = 1)) %>%
add_pie(labels = ~Letter, values = ~Dutch, name = "Dutch",
textposition = 'inside', textinfo = 'label+percent', showlegend = FALSE,
domain = list(row = 1, column = 0)) %>%
add_pie(labels = ~Letter, values = ~Esperanto, name = "Esperanto",
textposition = 'inside', textinfo = 'label+percent', showlegend = FALSE,
domain = list(row = 1, column = 1)) %>%
add_annotations(x=0.2, y=0.78,text = "German",showarrow = F, ax = 20, ay = -40) %>%
add_annotations(x=0.8, y=0.78,text = "Italian",showarrow = F, ax = 20, ay = -40) %>%
add_annotations(x=0.2, y=0.21,text = "Dutch",showarrow = F, ax = 20, ay = -40) %>%
add_annotations(x=0.82, y=0.21,text = "Esperanto",showarrow = F, ax = 20, ay = -40) %>%
layout(title = 'Pie Charts of German, Italian, Dutch & Esperanto Letters',
grid=list(rows=2, columns=2),
xaxis = list(showgrid = FALSE, zeroline = FALSE, showticklabels = FALSE),
yaxis = list(showgrid = FALSE, zeroline = FALSE, showticklabels = FALSE))The input type for pie() is a vector of non-negative
numerical quantities. In the pie function we list the data
that we are going to use (positive and numeric), the labels for each of
them, and the colors we want to use for each sector. In the
legend function, we put the location in the first slot and
legend are the labels for colors. cex,
bty, pch, and pt.cex are all
graphic parameters that we have talked about in Chapter
1.
More elaborate pie charts, using the Latin letter data, will be
demonstrated using ggplot later, (Section
7.2.
2.3.3 Heat map
Another common data visualization method is the
heat map. Heat maps can help us visualize the individual
values in a matrix intuitively. It is widely used in genetics research
and financial applications.
We will illustrate the use of heat maps, based on a neuroimaging genetics case-study data about the association (p-values) of different brain regions of interest (ROIs) and genetic traits (SNPs) for Alzheimer’s disease (AD) patients, subjects with mild cognitive impairment (MCI), and normal controls (NC). First, let’s import the data into R. The data are 2D arrays where the rows represent different genetic SNPs, columns represent brain ROIs, and the cell values represent the strength of the SNP-ROI association as probability values (smaller p-values indicate stronger neuroimaging-genetic associations).
AD_Data <- read.table("https://umich.instructure.com/files/330387/download?download_frd=1", header=TRUE, row.names=1, sep=",", dec=".")
MCI_Data <- read.table("https://umich.instructure.com/files/330390/download?download_frd=1", header=TRUE, row.names=1, sep=",", dec=".")
NC_Data <- read.table("https://umich.instructure.com/files/330391/download?download_frd=1", header=TRUE, row.names=1, sep=",", dec=".") Then we load the R packages we need for heat maps (use
install.packages("package name") first if you did not
install them into your computer).
Then we convert the datasets into matrices.
AD_mat <- as.matrix(AD_Data); class(AD_mat) <- "numeric"
MCI_mat <- as.matrix(MCI_Data); class(MCI_mat) <- "numeric"
NC_mat <- as.matrix(NC_Data); class(NC_mat) <- "numeric"We may also want to set up the row (rc) and column (cc) colors for each cohort.
rcAD <- rainbow(nrow(AD_mat), start = 0, end = 1.0); ccAD<-rainbow(ncol(AD_mat), start = 0, end = 1.0)
rcMCI <- rainbow(nrow(MCI_mat), start = 0, end=1.0); ccMCI<-rainbow(ncol(MCI_mat), start=0, end=1.0)
rcNC <- rainbow(nrow(NC_mat), start = 0, end = 1.0); ccNC<-rainbow(ncol(NC_mat), start = 0, end = 1.0)Finally, we got to the point where we can plot heat maps. As we can
see, the input type of heatmap() is a numeric matrix.
# hvAD <- heatmap(AD_mat, col = cm.colors(256), scale = "column", RowSideColors = rcAD, ColSideColors = ccAD, margins = c(2, 2), main="AD Cohort")
# hvMCI <- heatmap(MCI_mat, col = cm.colors(256), scale = "column", RowSideColors = rcMCI, ColSideColors = ccMCI, margins = c(2, 2), main="MCI Cohort")
# hvNC <- heatmap(NC_mat, col = cm.colors(256), scale = "column", RowSideColors = rcNC, ColSideColors = ccNC, margins = c(2, 2), main="NC Cohort")
# if (!require("devtools")) install.packages("devtools")
# devtools::install_github("talgalili/d3heatmap")
# library(d3heatmap)
# d3heatmap(AD_mat, dendrogram = 'both', key = TRUE, col = 'Blues', scale = 'column', key.title = "Legend",
# print.values = T, notecol = 'white') %>%
# hmAxis("x", title = "Imaging Phenotype", location = 'bottom') %>%
# hmAxis("y", title = "Genotype", location = 'left') %>%
# hmCells(font.size = 9, color = 'blue') %>%
# hmLegend(show = T, title = "AD Cohort", location = "tl")
plot_ly(x =~colnames(AD_mat), y = ~rownames(AD_mat), z = ~AD_mat, type = "heatmap") %>%
layout(title="AD Neuroimaging-Genomic Associations (p-values)",
xaxis=list(title="ROI Imaging Biomarkers"), yaxis=list(title="SNPs"))# d3heatmap(MCI_mat, dendrogram = 'both', key = TRUE, col = 'Blues', scale = 'column', key.title = "Legend",
# print.values = T, notecol = 'white') %>%
# hmAxis("x", title = "Imaging Phenotype", location = 'bottom') %>%
# hmAxis("y", title = "Genotype", location = 'left') %>%
# hmCells(font.size = 9, color = 'blue') %>%
# hmLegend(show = T, title = "MCI Cohort", location = "tl")
plot_ly(x =~colnames(MCI_mat), y = ~rownames(MCI_mat), z = ~MCI_mat, type = "heatmap") %>%
layout(title="MCI Neuroimaging-Genomic Associations (p-values)",
xaxis=list(title="ROI Imaging Biomarkers"), yaxis=list(title="SNPs"))# d3heatmap(NC_mat, dendrogram = 'both', key = TRUE, col = 'Blues', scale = 'column', key.title = "Legend",
# print.values = T, notecol = 'white') %>%
# hmAxis("x", title = "Imaging Phenotype", location = 'bottom') %>%
# hmAxis("y", title = "Genotype", location = 'left') %>%
# hmCells(font.size = 9, color = 'blue') %>%
# hmLegend(show = T, title = "Normal Cohort", location = "tl")
plot_ly(x =~colnames(NC_mat), y = ~rownames(NC_mat), z = ~NC_mat, type = "heatmap") %>%
layout(title="(Normal) HC Neuroimaging-Genomic Associations (p-values)",
xaxis=list(title="ROI Imaging Biomarkers"), yaxis=list(title="SNPs"))In the heatmap() function the first argument is for
matrices we want to use. col is the color scheme;
scale is a character indicating if the values should be
centered and scaled in either the row direction or the column direction,
or none (“row”, “column”, and “none”); RowSideColors and
ColSideColors creates the color names for horizontal side
bars.
The differences between the AD, MCI and NC heat maps are suggestive of variations of genetic traits or alternative brain regions that may be affected in the three clinically different cohorts.
2.4 Comparison
Plots used for comparing different individuals, groups of subjects, or multiple units represent another set of popular exploratory visualization tools.
2.4.1 Paired Scatter Plots
Scatter plots use the 2D Cartesian plane to display a graph indexed
by a pair of variables. 2D points in the graph represent values
associated with the two variables corresponding to the two coordinate
axes. The position of each 2D point is determined by the values of the
first and second variables, tracked on the horizontal and vertical axes.
If no clear dependent variable exists, either variable can be plotted on
the X axis and the corresponding scatter plot will illustrate the degree
of correlation (not necessarily causation) between two variables.
Although we will mostly demonstrate the use of plot_ly(),
which provides dynamic and interactive charts, many basic graphs,
including scatter plots, can be rendered using the R function
plot(x, y).
N <- 50
ind <- c(1:N)
x<-runif(N)
y<-runif(N)
z<-runif(N)
hoverText <- paste0("Point ", ind, ": (", round(x, 3), ",", round(y, 3), ")")
# plot(x, y, main="Scatter Plot")
plot_ly(x=~x[1:20], y=~y[1:20], type="scatter", size=2, name=ind[1:20],
color=~z[1:20], mode="markers", text = hoverText[1:20]) %>%
layout(title="Random Scatterplot", xaxis=list(title="X"), yaxis=list(title="Y")) %>%
hide_colorbar()# `qplot()` is another way to plot fancy scatter plots. We can manage the colors and sizes of dots. The input type for `qplot()` is a data frame. In the following example, larger *x* will have larger dot sizes. We also grouped the data as 10 points per group.
#
# library(ggplot2)
# cat <- rep(c("A", "B", "C", "D", "E"), 10)
# plot.1 <- qplot(x, y, geom="point", size=5*x, color=cat, main="GGplot with Relative Dot Size and Color")
# print(plot.1)Now let’s draw a paired scatter plot with 5 variables.
# The input type for `pairs()` function is a matrix or data frame.
# pairs(data.frame(x, y, z))
N=1000
w<-rnorm(N)
u<-rpois(N, lambda = 1.7)
# generate some random categorical labels for all N observations
class <- sample( LETTERS[1:3], N, replace=TRUE, prob=c(0.2, 0.5, 0.3))
df <- as.data.frame(cbind(x=x,y=y,z=z,w=w,u=u, class=class))
pl_colorscale=list(c(0.0, '#19d3f3'), c(0.333, '#19d3f3'), c(0.333, '#e763fa'), c(0.666, '#e763fa'),
c(0.666, '#636efa'), c(1, '#636efa'))
axis = list(showline=FALSE, zeroline=FALSE, gridcolor='#ffff', ticklen=4)
plot_ly(df) %>%
add_trace(type = 'splom', dimensions = list( list(label='X', values=~x), list(label='Y', values=~y),
list(label='Z', values=~z), list(label='w', values=~w), list(label='U', values=~u)),
text=~class,
marker = list(color = as.integer(df$class), colorscale = pl_colorscale,
size = 7, line = list(width = 1, color = 'rgb(230,230,230)')
)
) %>%
layout(
title= 'Random Data Pairs Plot', hovermode='closest', dragmode= 'select',
plot_bgcolor='rgba(240,240,240, 0.95)',
xaxis=list(domain=NULL, showline=F, zeroline=F, gridcolor='#ffff', ticklen=4),
yaxis=list(domain=NULL, showline=F, zeroline=F, gridcolor='#ffff', ticklen=4),
xaxis2=axis, xaxis3=axis, xaxis4=axis,yaxis2=axis, yaxis3=axis, yaxis4=axis)This is an interactive scatter plot where you can select/subset some observations in any of the plots and see their associations with other variables across all pairs plots.
Let’s see a real word data example. First, we can import the Mental Health Services Survey Data into R, which is on the class website. This survey data covers \(10,374\) mental health facilities across the US, the District of Columbia, and US Territories with 237 variables about various facility characteristics. A subset of 10 variables is included in this dataset with all 10,374 cases. Two of the facilitate characteristics involve (1) supp, representing the number of specialty and support services available at the mental health facility; and (2) qual, which is the number of quality indicators present at the mental health facility.
data1 <- read.table('https://umich.instructure.com/files/399128/download?download_frd=1', header=T)
head(data1)## STFIPS majorfundtype FacilityType Ownership Focus PostTraum GLBT num
## 1 southeast 1 5 2 1 0 0 5
## 2 southeast 3 5 3 1 0 0 4
## 3 southeast 1 6 2 1 1 1 9
## 4 greatlakes NA 2 2 1 0 0 7
## 5 rockymountain 1 5 2 3 0 0 9
## 6 mideast NA 2 2 1 0 0 8
## qual supp
## 1 NA NA
## 2 15 4
## 3 15 NA
## 4 14 6
## 5 18 NA
## 6 14 NAWe can see from head() that there are a lot of
NA’s in the dataset and the pairs plot (splom)
automatically ignores these (and posts a warning message).
# plot(data1[, 9], data1[, 10], pch=20, col="red", main="qual vs supp")
# pairs(data1[, 5:10])
plot_ly(data1, x=~qual, y=~supp, type="scatter", size=2, name=STFIPS,
color=~num, mode="markers", text = STFIPS) %>%
layout(title="2010 National Mental Health Services Survey: Support Services vs. Quality Indicators Scatterplot",
xaxis=list(title="Support Services"), yaxis=list(title="Quality Indicators")) %>%
hide_colorbar()plot_ly(data1) %>%
add_trace(type = 'splom', dimensions = list( list(label='FacilityType', values=~FacilityType ),
list(label='Ownership', values=~Ownership), list(label='Focus', values=~Focus),
list(label='PostTraum', values=~PostTraum), list(label='num', values=~num)),
text=~STFIPS,
marker = list(color = as.integer(qual), colorscale = pl_colorscale,
size = 7, line = list(width = 1, color = qual)
)
) %>%
layout(
title= '2010 National Mental Health Services Survey Pairs Plot (color=qual)', hovermode='closest', dragmode= 'select',
plot_bgcolor='rgba(240,240,240, 0.95)',
xaxis=list(domain=NULL, showline=F, zeroline=F, gridcolor='#ffff', ticklen=4),
yaxis=list(domain=NULL, showline=F, zeroline=F, gridcolor='#ffff', ticklen=4),
xaxis2=axis, xaxis3=axis, xaxis4=axis,yaxis2=axis, yaxis3=axis, yaxis4=axis)The first plot shows the relation between supp (support services) and qual (quality indicators). The more elaborate pairs plot illustrates multiple bivariate relations that can be interactively explored by selecting points in any of the plots, where points are color-coded by the quality indicator variable.
To see this trend model (loess(supp ~ qual) exposing the
trajectory of the support-services to quality relationship. This
locally estimated scatterplot smoothing (LOESS) model
represents a nonlinear smoothing regression.
# plot.2 <- qplot(qual, supp, data = data1, geom = c("point", "smooth"))
# print(plot.2)
# extract only the complete cases
library(dplyr)
df1 <- data1 %>% filter_at(vars(qual,supp), all_vars(!is.na(.)))
ll.smooth = loess(df1$supp ~ df1$qual, span=0.7)
ll.pred = predict(ll.smooth, se = TRUE)
ll.df = data.frame(x=ll.smooth$x, fit=ll.pred$fit, lb=ll.pred$fit-(1.96*ll.pred$se),
ub=ll.pred$fit+(1.96*ll.pred$se))
ll.df = ll.df[order(ll.df$df1.qual),]
plot_ly(x=df1$qual, y=df1$supp, type="scatter", mode="markers", name="Data") %>%
add_lines(x=df1$qual, y=ll.pred$fit, name="Mean", line=list(color="gray", width=4)) %>%
add_ribbons(x=ll.df$df1.qual, ymin=ll.df$lb, ymax=ll.df$ub, name="95% CI",
line=list(opacity=0.4, width=1, color="lightgray")) %>%
layout(title = "LOESS Model (Supp ~ Qual) with Confidence Band",
xaxis=list(title="Quality Indicator"), yaxis=list(title="Supporting Services"))You can also use the human height and weight dataset or the knee pain dataset to illustrate some interesting scatter plots.
2.4.2 Jitter plot
Jitter plot can help us deal with the overplot issue when we have
many points in the data. The function we will be using is still in the
package ggplot2 called position_jitter().
Still we use the earthquake data for example. We will compare the
differences with and without the position_jitter()
function.
# library("xml2"); library("rvest")
wiki_url <- read_html("https://wiki.socr.umich.edu/index.php/SOCR_Data_Dinov_021708_Earthquakes")
html_nodes(wiki_url, "#content")## {xml_nodeset (1)}
## [1] <div id="content" class="mw-body" role="main">\n\t\t\t<a id="top"></a>\n\ ...earthquake <- html_table(html_nodes(wiki_url, "table")[[2]])
# plot6.1<-ggplot(earthquake, aes(Depth, Latitude, group=Magt, color=Magt))+geom_point()
# plot6.2<-ggplot(earthquake, aes(Depth, Latitude, group=Magt, color=Magt))+geom_point(position = position_jitter(w = 0.3, h = 0.3), alpha=0.5)
# print(plot6.1)
# print(plot6.2)
# Note that with option `alpha=0.5` the "crowded" places are darker than the places with only one data point.
# Sometimes, we need to add text to these points, i.e., add label in `aes` or add `geom_text`. It looks messy.
# ggplot(earthquake, aes(Depth, Latitude, group=Magt, color=Magt,label=rownames(earthquake)))+
# geom_point(position = position_jitter(w = 0.3, h = 0.3), alpha=0.5)+geom_text()
# Let's try to fix the overlap of points and labels. We need to add `check_overlap` in `geom_text` and adjust the positions of the text labels with respect to the points.
#
# ```{r warning=FALSE, message=FALSE, error=FALSE}
# ggplot(earthquake, aes(Depth, Latitude, group=Magt, color=Magt,label=rownames(earthquake)))+
# geom_point(position = position_jitter(w = 0.3, h = 0.3), alpha=0.5)+
# geom_text(check_overlap = T,vjust = 0, nudge_y = 0.5, size = 2,angle = 45)
#
# # Or you can simply use the text to denote the positions of points.
# ggplot(earthquake, aes(Depth, Latitude, group=Magt, color=Magt,label=rownames(earthquake)))+
# geom_text(check_overlap = T,vjust = 0, nudge_y = 0, size = 3,angle = 45)
# # Warning: check_overlap will not show those overlapped points. Thus, if you need an analysis at the level of every instance, do not use it.
glyphication <- function (name) {
glyph= vector()
for (i in 1:length(name)){
glyph[i]="triangle-up"
if (name[i]=="Md") { glyph[i]="diamond-open" }
else if (name[i]=="ML") { glyph[i]="circle-open" }
else if (name[i]=="Mw") { glyph[i]="square-open" }
else if (name[i]=="Mx") { glyph[i]="x-open" }
}
return(glyph)
}
earthquake$glyph <- glyphication(earthquake$Magt)
plot_ly(earthquake) %>%
add_markers(x = ~Longitude, y = ~Latitude, type = "scatter", color = ~Magt,
mode = "markers", marker = list(size = ~Depth, color = ~Magt, symbol = ~glyph,
line = list(color = ~Magt, width = 3))) %>%
layout(title="California Earthquakes (1969 - 2007)")2.4.3 Bar Plots
Bar plots, or bar charts, represent group data with rectangular bars. There are many variants of bar charts for comparison among categories. Typically, either horizontal or vertical bars are used where one of the axes shows the compared categories and the other axis represents a discrete value. It’s possible, and sometimes desirable, to plot bar graphs including bars clustered by groups.
In R we can use plotly or barplot() for
barplots with inputs either vectors or matrices. The
ggplot2::diamonds dataset is comprised of \(53,940\) diamond records (rows) with 10
observed characteristics: price ($326–$18,823); carat (diamond weight);
cut (quality of the cut); color (D (best) to J (worst)); clarity (I1
(worst), …, IF (best)); x, and z length in mm; depth (total depth
percentage = z/mean(x, y) = 2*z/(x + y)); and table (diamond width of
top).
We can add error-bars to each bar to indicate a statistical variability. T
# bar <- barplot(m <- rowMeans(x) * 10, ylim=c(0, 10))
# stdev <- sd(t(x[1:4, ]))
# arrows(bar, m, bar, m + stdev, length=0.15, angle = 90)
plot_ly(ggplot2::diamonds, y = ~log(price), color=~cut, type = "box") %>%
layout(title = "Boxplot of Diamond (log) Price by Cut",
xaxis=list(title="Diamond Cut"))plot_ly(ggplot2::diamonds, x= ~clarity, y = ~log(price), color=~color, type = "box") %>%
layout(boxmode = "group", title = "Grouped Boxplot of Diamond (log) Price by Clarity and Color",
legend=list(title=list(text='<b> Diamond Color </b>')),
xaxis=list(title="Diamond Clarity"))Let’s look at a more complex example. We utilize the dataset Case_04_ChildTrauma for illustration. This case study examines associations between post-traumatic psychopathology and service utilization by trauma-exposed children.
data2 <- read.table('https://umich.instructure.com/files/399129/download?download_frd=1', header=T)
attach(data2)
head(data2)## id sex age ses race traumatype ptsd dissoc service
## 1 1 1 6 0 black sexabuse 1 1 17
## 2 2 1 14 0 black sexabuse 0 0 12
## 3 3 0 6 0 black sexabuse 0 1 9
## 4 4 0 11 0 black sexabuse 0 1 11
## 5 5 1 7 0 black sexabuse 1 1 15
## 6 6 0 9 0 black sexabuse 1 0 6We have two character variables. Our goal is to draw a bar plot
comparing the means of age and service among
different races in this study and we want to add standard deviation for
each bar. The first thing to do is delete the two character columns.
Remember the input for barplot() is numerical vector or
matrix. However, we will need race information for classification. Thus,
we store it in a different dataset.
Then, we are ready to separate groups and get group means.
data2.df <- as.data.frame(data2)
Blacks <- data2[which(data2$race=="black"), ]
Other <- data2[which(data2$race=="other"), ]
Hispanic <- data2[which(data2$race=="hispanic"), ]
White <- data2[which(data2$race=="white"), ]
B <- c(mean(Blacks$age), mean(Blacks$service))
O <- c(mean(Other$age), mean(Other$service))
H <- c(mean(Hispanic$age), mean(Hispanic$service))
W <- c(mean(White$age), mean(White$service))
x <- cbind(B, O, H, W)
x## B O H W
## [1,] 9.165 9.12 8.67 8.950000
## [2,] 9.930 10.32 9.61 9.911667Until now, we had a numerical matrix for the means available for plotting. Now, we can compute a second order statistics - standard deviation, and plot it along with the means, to illustrate the amount of dispersion for each variable.
# bar <- barplot(x, ylim=c(0, max(x)+2.0), beside=TRUE,
# legend.text = c("age", "service") , args.legend = list(x = "right"))
# text(labels=round(as.vector(as.matrix(x)), 2), x=seq(1.4, 21, by=1.5), #y=as.vector(as.matrix(x[1:2, ]))+0.3)
# y=11.5)
#
# m <- x; stdev <- sd(t(x))
# arrows(bar, m, bar, m + stdev, length=0.15, angle = 90)
# Here, we want the y margin to be little higher than the greatest value (`ylim=c(0, max(x)+2.0)`) because we need to leave space for value labels. Now we can easily notice that Hispanic trauma-exposed children are the youngest in terms of average age and they are less likely to utilize services like primary care, emergency room, outpatient therapy, outpatient psychiatrist, etc.
# Diamonds Dataset example
# data_mean <- ddply(diamonds, c("clarity", "cut"), summarize, price = mean(price))
# data_sd <- ddply(diamonds, c("clarity", "cut"), summarize, price = sd(price))
# data2 <- data.frame(data_mean, sd=data_sd$price)
#
# plot_ly(data = data2[which(data2$cut == 'Ideal'), ], x = ~clarity, y = ~price, type = 'bar',
# name = 'Cut=Ideal', error_y = ~list(array = sd, color = '#000000')) %>%
# add_trace(data = data2[which(data2$cut == 'Premium'), ], name = 'Cut=Premium') %>%
# add_trace(data = data2[which(data2$cut == 'Very Good'), ], name = 'Cut=Very Good') %>%
# add_trace(data = data2[which(data2$cut == 'Good'), ], name = 'Cut=Good') %>%
# add_trace(data = data2[which(data2$cut == 'Fair'), ], name = 'Cut=Fair') %>%
# layout(title="Statistical Barplots (Diamonds Dataset)",
# legend=list(title=list(text='<b> Diamond Cuts </b>')))
library(plyr)
data_mean <- ddply(data2, c("traumatype", "race"), summarise, service = mean(service))
data_sd <- ddply(diamonds, c("traumatype", "race"), summarise, service = sd(service))
data2 <- data.frame(data_mean, sd=data_sd$service)
plot_ly(data = data2[which(data2$race == 'black'), ], x = ~traumatype, y = ~service, type = 'bar',
name = 'Black', error_y = ~list(array = sd, color = '#000000')) %>%
add_trace(data = data2[which(data2$race == 'hispanic'), ], name = 'Hispanic') %>%
add_trace(data = data2[which(data2$race == 'other'), ], name = 'Other') %>%
add_trace(data = data2[which(data2$race == 'white'), ], name = 'White') %>%
layout(title="Statistical Barplots (Child Trauma Dataset)",
legend=list(title=list(text='<b> Race </b>')))Another way to plot bar plots is to use ggplot() in the
ggplot package. This kind of bar plots are quite different from the one
we introduced previously. It plots the counts of character variables
rather than the means of numerical variables. It takes the values from a
data.frame. Unlike barplot(), drawing bar
plots using ggplot2 requires remaining character variables
in the original data frame.
library(ggplot2)
#data2 <- read.table('https://umich.instructure.com/files/399129/download?download_frd=1', header=T)
ggplot(data2, aes(race, fill=race)) + geom_bar()+facet_grid(. ~ traumatype) 
This plot helps us to compare the occurrence of different types of child-trauma among different races.
2.4.4 Trees and Graphs
In general, a graph is an ordered pair \(G = (V, E)\) of vertices (\(V\)). i.e., nodes or points, and a set edges (\(E\)), arcs or lines connecting pairs of nodes in \(V\). A tree is a special type of acyclic graph that does not include looping paths. Visualization of graphs is critical in many biosocial and health studies and we will see examples throughout this textbook.
In Chapter 3 and Chapter 8 we will learn more about how to build tree models and other clustering methods, and in Chapter 22, we will discuss deep learning and neural networks, which intrinsically represent AI decision graphs.
This section will be focused on displaying tree graphs. We will use 02_Nof1_Data.csv for this demonstration.
data3<- read.table("https://umich.instructure.com/files/330385/download?download_frd=1", sep=",", header = TRUE)
head(data3)## ID Day Tx SelfEff SelfEff25 WPSS SocSuppt PMss PMss3 PhyAct
## 1 1 1 1 33 8 0.97 5.00 4.03 1.03 53
## 2 1 2 1 33 8 -0.17 3.87 4.03 1.03 73
## 3 1 3 0 33 8 0.81 4.84 4.03 1.03 23
## 4 1 4 0 33 8 -0.41 3.62 4.03 1.03 36
## 5 1 5 1 33 8 0.59 4.62 4.03 1.03 21
## 6 1 6 1 33 8 -1.16 2.87 4.03 1.03 0We use hclust to build the hierarchical cluster model.
hclust takes only inputs that have dissimilarity structure
as produced by dist(). Also, we use the ave()
method for agglomeration and plot our first tree graph.
When we have no limit for maximum cluster groups, we will get the
above graph, which is miserable to look at. Luckily, cutree
will help us to set limitations to the number of clusters.
cutree() takes a hclust object and returns a
vector of group indicators for all observations.
require(graphics)
mem <- cutree(hc, k = 10)
# mem; # to print the hierarchical tree labels for each case
# which(mem==5) # to identify which cases belong to class/cluster 5
# To see the number of Subjects in which cluster:
# table(cutree(hc, k=5))Then, we can get the mean of each variable within groups by the following for loop.
Now we can plot the new tree graph with 10 groups. With
members=table(mem) option, the matrix is taken to be a
dissimilarity matrix between clusters instead of dissimilarities between
singletons and members giving the number of observations per
cluster.
hc1 <- hclust(dist(cent), method = "ave", members = table(mem))
plot(hc1, hang = -1, main = "Re-start from 10 clusters")
2.4.5 Correlation Plots
The corrplot package enables the graphical display of a
correlation matrix, and confidence intervals, along with some tools for
matrix reordering. There are seven visualization methods (parameter
method) in the corrplot package, named “circle”, “square”,
“ellipse”, “number”, “shade”, “color”, “pie”.
Let’s use 03_NC_SNP_ROI_Assoc_P_values.csv again to investigate the associations among SNPs using correlation plots.
The corrplot() function we will be using takes
correlation matrix only. So we need to get the correlation matrix of our
data first via the cor() function.
# install.packages("corrplot")
library(corrplot)
NC_Associations_Data <- read.table("https://umich.instructure.com/files/330391/download?download_frd=1", header=TRUE, row.names=1, sep=",", dec=".")
M <- cor(NC_Associations_Data)
M[1:10, 1:10]## P2 P5 P9 P12 P13 P14
## P2 1.00000000 -0.05976123 0.99999944 -0.05976123 0.21245299 -0.05976123
## P5 -0.05976123 1.00000000 -0.05976131 -0.02857143 0.56024640 1.00000000
## P9 0.99999944 -0.05976131 1.00000000 -0.05976131 0.21248635 -0.05976131
## P12 -0.05976123 -0.02857143 -0.05976131 1.00000000 -0.05096471 -0.02857143
## P13 0.21245299 0.56024640 0.21248635 -0.05096471 1.00000000 0.56024640
## P14 -0.05976123 1.00000000 -0.05976131 -0.02857143 0.56024640 1.00000000
## P15 -0.08574886 0.69821536 -0.08574898 -0.04099594 0.36613665 0.69821536
## P16 -0.08574886 0.69821536 -0.08574898 -0.04099594 0.36613665 0.69821536
## P17 -0.05976123 -0.02857143 -0.05976131 -0.02857143 -0.05096471 -0.02857143
## P18 -0.05976123 -0.02857143 -0.05976131 -0.02857143 -0.05096471 -0.02857143
## P15 P16 P17 P18
## P2 -0.08574886 -0.08574886 -0.05976123 -0.05976123
## P5 0.69821536 0.69821536 -0.02857143 -0.02857143
## P9 -0.08574898 -0.08574898 -0.05976131 -0.05976131
## P12 -0.04099594 -0.04099594 -0.02857143 -0.02857143
## P13 0.36613665 0.36613665 -0.05096471 -0.05096471
## P14 0.69821536 0.69821536 -0.02857143 -0.02857143
## P15 1.00000000 1.00000000 -0.04099594 -0.04099594
## P16 1.00000000 1.00000000 -0.04099594 -0.04099594
## P17 -0.04099594 -0.04099594 1.00000000 -0.02857143
## P18 -0.04099594 -0.04099594 -0.02857143 1.00000000We will discover the difference among different methods under
corrplot.
# par specs c(bottom, left, top, right) which gives the margin size specified in inches
corrplot(M, method = "square", title = "square", tl.cex = 0.5, tl.col = 'black', mar=c(1, 1, 1, 1))
corrplot(M, method = "ellipse", title = "ellipse", tl.cex = 0.5, tl.col = 'black', mar=c(1, 1, 1, 1))
corrplot(M, type = "upper", tl.pos = "td",
method = "circle", tl.cex = 0.5, tl.col = 'black',
order = "hclust", diag = FALSE, mar=c(1, 1, 0, 1))
The shades are different and darker dots represent high correlation of the two variables corresponding to the x and y axes.
2.5 Relationships
2.5.1 Line plots using
ggplot
Line charts display a series of data points, e.g., observed intensities (\(Y\)) over time (\(X\)), by connecting them with straight-line segments. These can be used to either track temporal changes of a process or compare the trajectories of multiple cases, time series or subjects over time, space, or state.
In this section, we will utilize the Earthquakes dataset on SOCR website. It records information about earthquakes that occurred between 1969 and 2007 with magnitudes larger than 5 on the Richter scale.
# library("xml2"); library("rvest")
wiki_url <- read_html("https://wiki.socr.umich.edu/index.php/SOCR_Data_Dinov_021708_Earthquakes")
html_nodes(wiki_url, "#content")## {xml_nodeset (1)}
## [1] <div id="content" class="mw-body" role="main">\n\t\t\t<a id="top"></a>\n\ ...In this dataset, we set Magt(magnitude type) as groups.
We will draw a “Depth 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 and aes() specifies aesthetic mappings of how
variables in the data are mapped to visual properties (aesthetics) of
the geom objects, e.g., lines.
library(ggplot2)
plot4 <- ggplot(earthquake, aes(Longitude, Latitude, group=Magt, color=Magt))+
# Either draw lines
# geom_line()
# or, alternatively, we can draw glyphs/points
geom_point(data=earthquake, size=4, mapping=aes(x=Longitude, y=Latitude, shape=Magt))
plot4 # or print(plot4)
The first part
ggplot(earthquake, aes(Depth, Latitude, group=Magt, color=Magt))
in the code specifies the setting of the plot: dataset, group and color.
The second part specifies we are going to draw (points or) lines between
data points. In later chapters, we will frequently use the package
ggplot2 and the structure under this great package is
always function1+function2.
2.5.2 Density Plots
We can visualize the distribution for different variables using density plots.
The following segment of R code plots the distribution for latitude
among different earthquake
magnitude types. Also, it is using the ggplot()
function but combined with geom_density().
# library("ggplot2")
ggplot(earthquake, aes(Latitude, group=Magt, newsize=2))+geom_density(aes(color=Magt), size = 2) +
theme(legend.position = 'right',
legend.text = element_text(color= 'black', size = 12, face = 'bold'),
legend.key = element_rect(size = 0.5, linetype='solid'),
legend.key.size = unit(1.5, 'lines'))
Note how the green magt type (Local (ML) earthquakes)
has a peak at latitude \(37.5\), which
represents 37-38 degrees
North.
2.6 Distributions
Recall that there is a duality between theoretical and empirical mass, density, and distribution functions. Earlier, we saw the relations between these using the (continuous) Normal distribution, let’s now look at the (discrete) Poisson distribution. The graph below plots (1) the histogram of a sample of 1,000 Poisson(1) random observations (light blue color), (2) the theoretical density/mass function (magenta color), and (3) a smooth continuous (Gaussian) kernel density estimation based on the random sample (blue color). More interactive plots of univariate distributions and multivariate distributions are available online.
set.seed(1234)
poisson_sample <- rpois(1000, 1)
# slightly offset the histogram bins to align with mass function
hist_breakes <- c(-0.5, 0.5, 1.5, 2.5, 3.5, 6.5)
# hist(poisson_sample, freq=F, breaks = hist_breakes, col="light blue", lwd=2, ylim = c(0, 0.45))
# lines(density(poisson_sample, kernel = "gaussian"), lwd=2, col="blue")
# t <- seq(0, 6, by=0.01)
# lines(t, dpois(t,1), type="h", col="magenta", lwd=6) # add the theoretical density line
# legend(3,0.3, legend=c("Sample histogram (n=1,000)", "Theoretical mass function",
# "Gaussian kernel density estimate"),
# bty = "n", box.lty=0, col=c("light blue", "magenta", "blue"), lty=1, lwd=3)
h <-hist(poisson_sample, breaks = hist_breakes, plot = F)
t <- seq(0, 6, by=0.01)
Pois <- density(poisson_sample, kernel = "gaussian")
plot_ly(x = h$mids, y = h$density, type = "bar", name="Sample Histogram") %>%
add_lines(x=t, y=dpois(t,1), type="scatter", mode="lines",
name="(Theoretical) Poisson Mass Function") %>%
add_lines(x=Pois$x, y=Pois$y,
type="scatter", mode="lines",
name="Gaussian kernel density estimate (sample)") %>%
layout(bargap=0.1, title="Histogram (Simulated Poisson Data)",
legend = list(orientation = 'h'))2.6.1 Data Modeler
A common task in data-driven inference involves the fitting of appropriate distribution models to specific observed data elements (features). In general, as there are uncountably many possible distributions that can be used as models for various types of processes, this is a difficult task. The Probability Distributome Project (see Distributome Navigator) provides a deeper understanding of the notion of a probability distribution and the relations between various distributions.
We will demonstrate the concept of a data modeler by
using crystallographic
data from the Ivanova Lab
at the University of Michigan, which includes the crystal spectra of
9
length samples and 9
width samples. For both, the length and width spectra, the 9
features include “AC1338”, “AC1432”, “AC1593”, “AC1679”, “AC1860”,
“AC1874”, “AC1881”, “AC1903”, and “Rec” (these represent different
samples). Notice that the nine spectra are not congruent, different
features have different sampling rates. We will employ the fitdistrplus
R-package to estimate the parameters of 3 complementary
distributions, however, there are many alternative packages that can
also be used.
2.6.1.1 Loading the spectral crystallography data
The data include two separate signals capturing the spectral length and the width of the crystallographic sample.
- Dec 2019 crystallography spectral data
- crystallography Length data are here
- crystallography Width data are here
# You may choose which of the 2 CSV files (width or length) to work with
crystallography_Length_data <- read.csv(file = "https://umich.instructure.com/files/11653615/download?download_frd=1",
header=TRUE)
crystallography_Width_data <- read.csv(file = "https://umich.instructure.com/files/11653614/download?download_frd=1",
header=TRUE)
crystallography_data <- crystallography_Length_data
# crystallography_data <- crystallography_Width_data
# Get the feature names (IDs)
colNames <- colnames(crystallography_data); colNames## [1] "AC1338" "AC1432" "AC1593" "AC1679" "AC1860" "AC1874" "AC1881" "AC1903"
## [9] "Rec"2.6.1.2 Feature distributions
Let’s plot the histograms of each of the nine features.
# plot all histograms
library(tidyr)
# library(ggplot2)
# # or `library(tidyverse)`
#
# crystallography_data %>% gather() %>% head()
# # key value
# #1 AC1338 70.547
# #2 AC1338 40.448
# #3 AC1338 47.212
# #4 AC1338 91.468
# #5 AC1338 79.088
# #6 AC1338 132.319
# #...
# crystallography_data %>% gather() %>% tail()
# # key value
# #5872 Rec 68.479
# #5873 Rec 41.047
# #5874 Rec 47.546
# #5875 Rec 98.558
# #5876 Rec 52.956
# #5877 Rec 82.470
#
# ggplot(gather(crystallography_data), aes(value)) +
# geom_histogram(bins = 20) +
# facet_wrap(~key, scales = 'free_x')
crystalCompleteData <- crystallography_data[complete.cases(crystallography_data), ]
df_crystal <- apply(crystalCompleteData, 2, density, kernel="gaussian", bw=15)
df <- data.frame(x = unlist(lapply(df_crystal, "[[", "x")),
y = unlist(lapply(df_crystal, "[[", "y")),
sample = rep(names(df_crystal), each = length(df_crystal[[1]]$x)))
plot_ly(df, x = ~x, y = ~y, color = ~sample, type = "scatter", mode = "lines") %>%
layout(title='Crystallography Sample Densities',
legend=list(title=list(text='<b> Samples </b>')),
xaxis=list(title='X'), yaxis=list(title='Density'))2.6.1.3 Fitting single-feature univariate distribution models
We will fit Weibull, Gamma, and Log-Normal distribution models to each feature in the data.
# install.packages("fitdistrplus")
library(fitdistrplus)
col_num <- dim(crystallography_data)[2]; col_num## [1] 9# Store the Weibull, Gamma, and Log-Normal Distribution models for the 9 features
fit_W <- vector(mode = "list", length = col_num)
fit_G <- vector(mode = "list", length = col_num)
fit_LN <- vector(mode = "list", length = col_num)
for(i in 1:col_num) {
data_no_NA <- crystallography_data[complete.cases(crystallography_data[, i]), i]
length(data_no_NA)
fit_W[[i]] <- fitdist(data_no_NA, "weibull"); summary(fit_W[i])
fit_G[[i]] <- fitdist(data_no_NA, "gamma"); summary(fit_G[i])
fit_LN[[i]] <- fitdist(data_no_NA, "lnorm"); summary(fit_LN[i])
}
# extract the model parameters
W_mod_p1_name = array(dim=c(col_num,2)); dim(W_mod_p1_name) # param name## [1] 9 2## [1] 9 2## [1] 9 2## [1] 9 2## [1] 9 2## [1] 9 2# Compute the mean (m) and standard deviation (sd) for each model distribution
W_mod_mean = array(dim=c(col_num,1)); length(W_mod_mean) # Weibull mean or mode## [1] 9## [1] 9## [1] 9## [1] 9## [1] 9## [1] 9for(i in 1:col_num) {
W_mod_p1_name[i, 1] <- names(fit_W[[i]]$estimate[1]) # Weibull "shape"
W_mod_p1_val[i, 1] <- fit_W[[i]]$estimate[[1]]
W_mod_p1_name[i, 2] <- names(fit_W[[i]]$estimate[2]) # Weibull "scale"
W_mod_p1_val[i, 2] <- fit_W[[i]]$estimate[[2]]
W_mod_mean[i] = W_mod_p1_val[i, 2] * gamma(1+1/W_mod_p1_val[i, 1]) # Weibull mean
W_mod_mean[i] = W_mod_p1_val[i, 2] *
((W_mod_p1_val[i, 1]-1)/W_mod_p1_val[i, 1])^(1/W_mod_p1_val[i, 1]) # Weibull mode
W_mod_sd[i] = W_mod_p1_val[i, 2]*sqrt(gamma(1+2/W_mod_p1_val[i, 1])-
(gamma(1+1/W_mod_p1_val[i, 1]))^2) # Weibull SD
G_mod_p1_name[i, 1] <- names(fit_G[[i]]$estimate[1]) # Gamma "shape"
G_mod_p1_val[i, 1] <- fit_G[[i]]$estimate[[1]]
G_mod_p1_name[i, 2] <- names(fit_G[[i]]$estimate[2]) # Gamma "scale"
G_mod_p1_val[i, 2] <- fit_G[[i]]$estimate[[2]]
G_mod_mean[i] = G_mod_p1_val[i, 1] / G_mod_p1_val[i, 2] # Gamma mean
G_mod_mean[i] = (G_mod_p1_val[i, 1]-1) / G_mod_p1_val[i, 2] # Gamma mode
G_mod_sd[i] = sqrt(G_mod_p1_val[i, 1]) / G_mod_p1_val[i, 2] # Gamma SD
LN_mod_p1_name[i, 1] <- names(fit_LN[[i]]$estimate[1]) # Log-normal "shape"
LN_mod_p1_val[i, 1] <- fit_LN[[i]]$estimate[[1]]
LN_mod_p1_name[i, 2] <- names(fit_LN[[i]]$estimate[2]) # Log-normal "scale"
LN_mod_p1_val[i, 2] <- fit_LN[[i]]$estimate[[2]]
LN_mod_mean[i] = exp(LN_mod_p1_val[i, 1]+ (LN_mod_p1_val[i, 2])^2/2) # Log-normal mean
LN_mod_mean[i] = exp(LN_mod_p1_val[i, 1] - LN_mod_p1_val[i, 2]^2) # Log-normal mean
LN_mod_sd[i] = sqrt((exp(LN_mod_p1_val[i, 2]^2)-1)*
exp(2*LN_mod_p1_val[i, 1]+LN_mod_p1_val[i, 2]^2)) # Log-normal SD
}
# Check results, just for one model
str(fit_W[[1]])## List of 17
## $ estimate : Named num [1:2] 2.12 96.21
## ..- attr(*, "names")= chr [1:2] "shape" "scale"
## $ method : chr "mle"
## $ sd : Named num [1:2] 0.074 2.251
## ..- attr(*, "names")= chr [1:2] "shape" "scale"
## $ cor : num [1:2, 1:2] 1 0.328 0.328 1
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:2] "shape" "scale"
## .. ..$ : chr [1:2] "shape" "scale"
## $ vcov : num [1:2, 1:2] 0.00548 0.05464 0.05464 5.06895
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:2] "shape" "scale"
## .. ..$ : chr [1:2] "shape" "scale"
## $ loglik : num -2308
## $ aic : num 4621
## $ bic : num 4629
## $ n : int 453
## $ data : num [1:453] 70.5 40.4 47.2 91.5 79.1 ...
## $ distname : chr "weibull"
## $ fix.arg : NULL
## $ fix.arg.fun: NULL
## $ dots : NULL
## $ convergence: int 0
## $ discrete : logi FALSE
## $ weights : NULL
## - attr(*, "class")= chr "fitdist"2.6.1.4 Visual inspection
Let’s examine graphically the quality of the fitted distribution models. We’ll plot the histograms of the features, the fitted probability densities, and the corresponding cumulative distribution functions (CDF) and compare them to their sample counterparts.
windows(width=20, height=8)
par(mfrow=c(3,3))
for(i in 1:col_num) {
# W_mod_p1_name[i] <- names(fit_W[[i]]$estimate[1])
# W_mod_p1_val[i] <- fit_W[[1]]$estimate[[1]]
plot.legend <- c(sprintf("Weibull(%s=%s,%s=%s) (m=%s,sd=%s)",
W_mod_p1_name[i, 1], format(W_mod_p1_val[i, 1], digits=2),
W_mod_p1_name[i, 2], format(W_mod_p1_val[i, 2], digits=2),
format(W_mod_mean[i], digits=2),
format(W_mod_sd[i], digits=2)),
sprintf("Gamma(%s=%s,%s=%s) (m=%s,sd=%s)",
G_mod_p1_name[i, 1], format(G_mod_p1_val[i, 1], digits=2),
G_mod_p1_name[i, 2], format(G_mod_p1_val[i, 2], digits=2),
format(G_mod_mean[i], digits=2),
format(G_mod_sd[i], digits=2)),
sprintf("Log-normal(%s=%s,%s=%s) (m=%s,sd=%s)",
LN_mod_p1_name[i, 1], format(LN_mod_p1_val[i, 1], digits=2),
LN_mod_p1_name[i, 2], format(LN_mod_p1_val[i, 2], digits=2),
format(LN_mod_mean[i], digits=2),
format(LN_mod_sd[i], digits=2)))
denscomp(list(fit_W[[i]], fit_G[[i]], fit_LN[[i]]), legendtext = plot.legend,
xlegend = "topright", ylegend ="right",
main=sprintf("Width: Feature: %s: Histogram & Model Densities", colnames(crystallography_data)[i]))
abline(v = format(W_mod_mean[i], digits=2), col = "red", lty=1)
abline(v = format(G_mod_mean[i], digits=2), col = "green", lty=2)
abline(v = format(LN_mod_mean[i], digits=2), col = "blue", lty=3)
# cdfcomp (list(fit_w, fit_g, fit_ln), legendtext = plot.legend)
# qqcomp (list(fit_w, fit_g, fit_ln), legendtext = plot.legend)
# ppcomp (list(fit_w, fit_g, fit_ln), legendtext = plot.legend)
}
# Plot histograms and CDF (cumulative distribution function) models
windows(width=20, height=12)
par(mfrow=c(3,3))
for(i in 1:col_num) {
plot.legend <- c(sprintf("Weibull(%s=%s,%s=%s) (m=%s,sd=%s)",
W_mod_p1_name[i, 1], format(W_mod_p1_val[i, 1], digits=2),
W_mod_p1_name[i, 2], format(W_mod_p1_val[i, 2], digits=2),
format(W_mod_mean[i], digits=2),
format(W_mod_sd[i], digits=2)),
sprintf("Gamma(%s=%s,%s=%s) (m=%s,sd=%s)",
G_mod_p1_name[i, 1], format(G_mod_p1_val[i, 1], digits=2),
G_mod_p1_name[i, 2], format(G_mod_p1_val[i, 2], digits=2),
format(G_mod_mean[i], digits=2),
format(G_mod_sd[i], digits=2)),
sprintf("Log-normal(%s=%s,%s=%s) (m=%s,sd=%s)",
LN_mod_p1_name[i, 1], format(LN_mod_p1_val[i, 1], digits=2),
LN_mod_p1_name[i, 2], format(LN_mod_p1_val[i, 2], digits=2),
format(LN_mod_mean[i], digits=2),
format(LN_mod_sd[i], digits=2)))
cdfcomp(list(fit_W[[i]], fit_G[[i]], fit_LN[[i]]), legendtext = plot.legend,
xlegend = "bottomright", ylegend ="right",
main=sprintf("Width: Feature: %s: Aggregate Hist & Model CDFs", colnames(crystallography_data)[i]))
}
Below is the plot_ly() version of the model fit for one
case.
pl_list <- list()
for(i in 1:col_num) {
# W_mod_p1_name[i] <- names(fit_W[[i]]$estimate[1])
# W_mod_p1_val[i] <- fit_W[[1]]$estimate[[1]]
plot.legend <- c(sprintf("Weibull(%s=%s,%s=%s) (m=%s,sd=%s)",
W_mod_p1_name[i, 1], format(W_mod_p1_val[i, 1], digits=2),
W_mod_p1_name[i, 2], format(W_mod_p1_val[i, 2], digits=2),
format(W_mod_mean[i], digits=2),
format(W_mod_sd[i], digits=2)),
sprintf("Gamma(%s=%s,%s=%s) (m=%s,sd=%s)",
G_mod_p1_name[i, 1], format(G_mod_p1_val[i, 1], digits=2),
G_mod_p1_name[i, 2], format(G_mod_p1_val[i, 2], digits=2),
format(G_mod_mean[i], digits=2),
format(G_mod_sd[i], digits=2)),
sprintf("Log-normal(%s=%s,%s=%s) (m=%s,sd=%s)",
LN_mod_p1_name[i, 1], format(LN_mod_p1_val[i, 1], digits=2),
LN_mod_p1_name[i, 2], format(LN_mod_p1_val[i, 2], digits=2),
format(LN_mod_mean[i], digits=2),
format(LN_mod_sd[i], digits=2)))
# x <- dweibull(10000, shape=fit_W[[i]]$estimate[1], scale =fit_W[[i]]$estimate[2])
# fit <- density(x)
z <- seq(from=min(fit_W[[i]]$data), max(fit_W[[i]]$data), 0.1) # points from -4 to 4 in 0.1 steps
weibullDens <- dweibull(z, shape=fit_W[[i]]$estimate[1], scale =fit_W[[i]]$estimate[2])
gammaDens <- dgamma(z, shape=fit_G[[i]]$estimate[1], rate =fit_G[[i]]$estimate[2])
logNormalDens <- dlnorm(z, meanlog=fit_LN[[i]]$estimate[1], sdlog =fit_LN[[i]]$estimate[2])
# z<-seq(from=min(fit_W[[i]]$data), to=max(fit_W[[i]]$data), 0.1) # Range points in 0.1 steps
pl_list[[i]] <-
plot_ly(x=~fit_W[[i]]$data, name=~colnames(crystallography_data)[i], showlegend = FALSE,
marker = list(color = "transparent", line = list(color = "darkgray", width = 2)),
type="histogram", mode="markers", opacity=0.9, nbinsx=20, histnorm="probability") %>%
# add models
add_trace(x=z, y=15*weibullDens, type="scatter", mode="lines", opacity=0.5, name=plot.legend[1],
line = list(color = "red", width = 2)) %>%
add_trace(x=z, y=15*gammaDens, type="scatter", mode="lines", opacity=0.5, name=plot.legend[2],
line = list(color = "green", width = 2)) %>%
add_trace(x=z, y=15*logNormalDens, type="scatter", mode="lines", opacity=0.5, name=plot.legend[3],
line = list(color = "blue", width = 2)) %>%
# add vertical mean lines
add_segments(x=W_mod_mean[i], y=0, xend=W_mod_mean[i], yend=0.2, name="Weibull mean", color="red") %>%
add_segments(x=G_mod_mean[i], y=0, xend=G_mod_mean[i], yend=0.2, name="Gamma mean", color="green") %>%
add_segments(x=LN_mod_mean[i], y=0, xend=LN_mod_mean[i], yend=0.2, name="LogNormal mean", color="blue") %>%
layout(title = sprintf("Width: Feature: %s: Histogram & Model Densities", colnames(crystallography_data)[i]),
xaxis = list(title = colnames(crystallography_data)[i]), yaxis = list(title = "Density"),
bargap=0.1) %>% hide_colorbar()
}
pl_list %>% plotly::subplot(nrows = 3) %>% layout(title="Mixture Modeling of Crystallography Data (Interactive Plot)") 2.6.1.5 Quantitative summaries
Often, it’s useful to export the numerical results of the models. This may include various distribution characteristics like measure of centrality (e.g., mean, median, mode), measures of dispersion, and metrics of the model performance (e.g., Kolmogorov-Smirnov test).
# Save the summary outputs (mode & SD) across 9 samples, 3 models and 2 measures into a dataframe
df_matrix = array(dim=c(col_num,3*2*2)); dim(df_matrix) ## [1] 9 12for(i in 1:col_num) {
data1 <- crystallography_data[complete.cases(crystallography_data[, i]), i]
df_matrix[i, 1] = format(W_mod_mean[i], digits=2) # Weibull mode
df_matrix[i, 2] = format(W_mod_sd[i], digits=2) # Weibull SD
ks_W <- ks.test(data1, "pweibull", scale=W_mod_p1_val[i, 2], shape=W_mod_p1_val[i, 1])
df_matrix[i, 3] = format(ks_W$statistic[[1]], digits=4) # KS-test-stat Weibull
df_matrix[i, 4] = format(ks_W$p.value, digits=5) # KS-test-p-value Weibull
df_matrix[i, 5] = format(G_mod_mean[i], digits=2) # Gamma mode
df_matrix[i, 6] = format(G_mod_sd[i], digits=2) # Gamma SD
ks_G <- ks.test(data1, "pgamma", rate=G_mod_p1_val[i, 2], shape=G_mod_p1_val[i, 1])
df_matrix[i, 7] = format(ks_G$statistic[[1]], digits=4) # KS-test-stat Gamma
df_matrix[i, 8] = format(ks_G$p.value, digits=5) # KS-test-p-value Gamma
df_matrix[i, 9] = format(LN_mod_mean[i], digits=2) # Log-normal mode
df_matrix[i, 10] = format(LN_mod_sd[i], digits=2) # Log-normal SD
ks_LN <- ks.test(data1, "plnorm", sdlog=LN_mod_p1_val[i, 2], meanlog=LN_mod_p1_val[i, 1])
df_matrix[i, 11] = format(ks_LN$statistic[[1]], digits=4) # KS-test-stat Log-normal
df_matrix[i, 12] = format(ks_G$p.value, digits=5) # KS-test-p-value Log-normal
}
df_summary <- as.data.frame(df_matrix, row.names=colNames)
colnames(df_summary) <- c("Weibull_mode", "Weibull_sd","Weibull_KS.test.stat", "Weibull_KS.p.val",
"Gamma_mode", "Gamma_sd","Gamma_KS.test.stat", "Gamma_KS.p.val",
"Lognormal_mode", "Lognormal_sd","Lognormal_KS.test.stat", "Lognormal_KS.p.val")
df_summary## Weibull_mode Weibull_sd Weibull_KS.test.stat Weibull_KS.p.val Gamma_mode
## AC1338 71 42 0.0411 0.4284 64
## AC1432 75 40 0.07218 0.047982 69
## AC1593 81 54 0.05572 0.10341 75
## AC1679 81 49 0.0462 0.36208 73
## AC1860 78 45 0.06798 0.088752 73
## AC1874 75 42 0.06495 0.032324 68
## AC1881 72 58 0.0821 0.00069318 70
## AC1903 80 48 0.07426 0.059275 73
## Rec 76 41 0.05729 0.027524 68
## Gamma_sd Gamma_KS.test.stat Gamma_KS.p.val Lognormal_mode Lognormal_sd
## AC1338 42 0.02878 0.84738 57 48
## AC1432 38 0.03942 0.63424 63 40
## AC1593 52 0.03823 0.4885 67 58
## AC1679 49 0.03222 0.80172 64 56
## AC1860 42 0.03691 0.74826 67 45
## AC1874 41 0.03431 0.61239 61 45
## AC1881 55 0.05289 0.073267 63 60
## AC1903 47 0.06417 0.14456 66 51
## Rec 40 0.03865 0.28357 62 44
## Lognormal_KS.test.stat Lognormal_KS.p.val
## AC1338 0.05412 0.84738
## AC1432 0.0315 0.63424
## AC1593 0.03584 0.4885
## AC1679 0.03622 0.80172
## AC1860 0.03832 0.74826
## AC1874 0.03334 0.61239
## AC1881 0.0294 0.073267
## AC1903 0.04493 0.14456
## Rec 0.03565 0.283572.6.1.6 Mixture distribution data modeling
Earlier, we discussed the expectations maximization (EM) algorithm for parameter estimation. Now, we will illustrate the use of EM to estimate the mixture weights and the distribution parameters needed to obtain mixture-distribution data models.
For each sample, we fit a mixture distribution of \(k=3\) (different number of distribution models, which is predefined). The specific types of mixtures for each of the 9 samples are indicated below.
2.6.1.7 Mixture-distribution model fitting and parameter estimation
We will use the R package mixtools to obtain the EM estimates of the mixture distribution weights and the corresponding distribution parameters.
# crystallography_data <- read.csv(file = "https://umich.instructure.com/files/13375767/download?download_frd=1",
# header=TRUE)
# crystallography_data <- read.csv(file = "https://umich.instructure.com/files/11653615/download?download_frd=1",
# header=TRUE)
# install.packages("mixtools")
library(mixtools)
col_num <- dim(crystallography_data)[2]; col_num## [1] 9# Fit mixture models
capture.output(
for(i in 1:col_num) { # remove all non-numeric elements (if any)
# data_no_NA <- unlist(Filter(is.numeric, crystallography_data[complete.cases(crystallography_data[, i]), i]))
data_no_NA <- crystallography_data[complete.cases(crystallography_data[, i]), i]
length(data_no_NA)
fit_W[[i]] <- weibullRMM_SEM(data_no_NA, k=df_sampleMixtureParam[1,i], verb=F)
# summary(fit_W[i])
fit_G[[i]] <- gammamixEM(data_no_NA, k=df_sampleMixtureParam[1,i], verb=F)
# summary(fit_G[i])
fit_LN[[i]] <- normalmixEM(data_no_NA, k=df_sampleMixtureParam[1,i], verb=F)
# summary(fit_LN[i])
},
file='NUL'
)
# plot(fit_LN[[1]], which=2)
# lines(density(crystallography_data[complete.cases(crystallography_data[, 1]), 1]), lty=2, lwd=2)2.6.1.8 Plotting the mixture distribution models
We will define custom plots for the mixtures of Gamma,
Weibull, and Normal distributions. Alternatively, we
can also use some of the mixtools::plot() function to
display mixture distribution models.
# Custom design of Gamma-Mixture Model plot
gammaMM.plot <- function(mix.object, k = 2, main = "") { # mix.object <- fit_G[[i]]
data_no_NA <- crystallography_data[complete.cases(crystallography_data[, i]), i]
d3 <- function(x) { # construct the mixture using the estimated parameters
mix.object$lambda[1]*dgamma(x, shape=mix.object$gamma.pars[1,1], 1/mix.object$gamma.pars[2,1]) +
mix.object$lambda[2]*dgamma(x, shape=mix.object$gamma.pars[1,2], 1/mix.object$gamma.pars[2,2]) +
mix.object$lambda[3]*dgamma(x, shape=mix.object$gamma.pars[1,3], 1/mix.object$gamma.pars[2,3])
}
x <- seq(min(data_no_NA), max(data_no_NA), 0.001)
hist(data_no_NA, col="pink", freq=F, breaks=10, main = main, xlab="Intensities")
lines(x, d3(x), lwd=3, col="black", xlim=c(4,23), ylim=c(0, 0.25))
mixColors <- colorRampPalette(c("blue", "red"))(k)
for (i in 1:k) {
d = function(x) { # construct each of the Gamma components using the estimated parameters
mix.object$lambda[i]*dgamma(x, shape=mix.object$gamma.pars[1, i], 1/mix.object$gamma.pars[2,i])
}
lines(x, d(x), lwd=3, col=mixColors[i])
}
}
# Custom design of Weibull-Mixture Model plot
weibullMM.plot <- function(mix.object, k = 2, main = "") { # mix.object <- fit_W[[i]]
data_no_NA <- crystallography_data[complete.cases(crystallography_data[, i]), i]
d3 <- function(x) { # construct the mixture using the estimated parameters
mix.object$lambda[1]*dweibull(x, shape=mix.object$shape[1], scale=mix.object$scale[1]) +
mix.object$lambda[2]*dweibull(x, shape=mix.object$shape[2], scale=mix.object$scale[2]) +
mix.object$lambda[3]*dweibull(x, shape=mix.object$shape[3], scale=mix.object$scale[3])
}
x <- seq(min(data_no_NA), max(data_no_NA), 0.001)
hist(data_no_NA, col="pink", freq=F, breaks=15, main = main, xlab="Intensities")
lines(x, d3(x), lwd=3, col="black", xlim=c(4,23), ylim=c(0, 0.25))
mixColors <- colorRampPalette(c("blue", "red"))(k)
for (i in 1:k) {
d = function(x) { # construct each of the Weibull components using the estimated parameters
mix.object$lambda[i]*dweibull(x, shape=mix.object$shape[i], scale=mix.object$scale[i])
}
lines(x, d(x), lwd=3, col=mixColors[i])
}
}
# Custom design of Normal-Mixture Model plot
normalMM.plot <- function(mix.object, k = 2, main = "") { # mix.object <- fit_LN[[i]]
data_no_NA <- crystallography_data[complete.cases(crystallography_data[, i]), i]
d3 <- function(x) { # construct the mixture using the estimated parameters
mix.object$lambda[1]*dnorm(x, mean=mix.object$mu[1], sd=mix.object$sigma[1]) +
mix.object$lambda[2]*dnorm(x, mean=mix.object$mu[2], sd=mix.object$sigma[2]) +
mix.object$lambda[3]*dnorm(x, mean=mix.object$mu[3], sd=mix.object$sigma[3])
}
x <- seq(min(data_no_NA), max(data_no_NA), 0.001)
hist(data_no_NA, col="pink", freq=F, breaks=20, main = main, xlab="Intensities", xlim = c(4,180), ylim = c(0.0, 0.02))
lines(x, d3(x), lwd=3, col="black")
mixColors <- colorRampPalette(c("blue", "red"))(k)
for (i in 1:k) {
d = function(x) { # construct each of the Normal components using the estimated parameters
mix.object$lambda[i]*dnorm(x, mean=mix.object$mu[i], sd=mix.object$sigma[i])
}
lines(x, d(x), lwd=3, col=mixColors[i])
}
}Next, we will display the three alternative mixture distribution models overlaid on the sample histograms of each of the nine samples.
# Plot Mixture Models and Report model parameter estimates
# for(i in 1:col_num) { # uncomment this to plot all 9 samples
for(i in 1:2) { # this only plots the first 2 samples to save space
weibullMM.plot(fit_W[[i]], df_sampleMixtureParam[1,i],
paste0("Mixture of ", df_sampleMixtureParam[1, sampleColNames[i]],
" Weibull Models of ", sampleColNames[i]))
#plot(fit_W[[i]], density=TRUE, whichplots = 2,
# main2=paste0("Mixture of ", df_sampleMixtureParam[1, sampleColNames[i]],
# " Weibull Models of ", sampleColNames[i]), xlab2="Intensities")
gammaMM.plot(fit_G[[i]], df_sampleMixtureParam[1,i],
paste0("Mixture of ", df_sampleMixtureParam[1, sampleColNames[i]],
" Gamma Models of ", sampleColNames[i]))
normalMM.plot(fit_LN[[i]], df_sampleMixtureParam[1,i],
paste0("Mixture of ", df_sampleMixtureParam[1, sampleColNames[i]],
" Normal Models of ", sampleColNames[i]))
}
2.6.1.9 Reporting model parameter estimates
For each of the 9 samples in this dataset) and each of the 3 types of mixture distribution models (Weibull, Gamma, and Normal) we will summarize:
- lambda: The weights (impacts) of each of the 3 mixture components to the overall mixture model,
- parameters: of each mixture distribution component, mean and sd,
- loglik: the overall mixture distribution log-likelihood value.
# Generate the summary DF
getSummaryTable <- function (crystalSampleIndex) {
mat <- matrix(0, nrow = 3, ncol = 10)
# Weibull estimates for all 3 model components
# For Weibull Dist mean and SD see: https://en.wikipedia.org/wiki/Weibull_distribution
mat[1,1] <- round(fit_W[[crystalSampleIndex]]$lambda[1],3) # lambda
mat[1,2] <- round(fit_W[[crystalSampleIndex]]$scale[1] *
gamma(1+1/fit_W[[crystalSampleIndex]]$shape[1]),3) # mean
mat[1,3] <- round(fit_W[[crystalSampleIndex]]$scale[1] *
sqrt(gamma(1+2/fit_W[[crystalSampleIndex]]$shape[1])-
(gamma(1+1/fit_W[[crystalSampleIndex]]$shape[1]))^2),3) # sd
mat[1,4] <- round(fit_W[[crystalSampleIndex]]$lambda[2],3) # lambda
mat[1,5] <- round(fit_W[[crystalSampleIndex]]$scale[2] *
gamma(1+1/fit_W[[crystalSampleIndex]]$shape[2]),3) # mean
mat[1,6] <- round(fit_W[[crystalSampleIndex]]$scale[2] *
sqrt(gamma(1+2/fit_W[[crystalSampleIndex]]$shape[2])-
(gamma(1+1/fit_W[[crystalSampleIndex]]$shape[2]))^2),3) # sd
mat[1,7] <- round(fit_W[[crystalSampleIndex]]$lambda[3],3) # lambda
mat[1,8] <- round(fit_W[[crystalSampleIndex]]$scale[3] *
gamma(1+1/fit_W[[crystalSampleIndex]]$shape[3]),3) # mean
mat[1,9] <- round(fit_W[[crystalSampleIndex]]$scale[3] *
sqrt(gamma(1+2/fit_W[[crystalSampleIndex]]$shape[3])-
(gamma(1+1/fit_W[[crystalSampleIndex]]$shape[3]))^2),3) # sd
mat[1,10] <- round(fit_W[[crystalSampleIndex]]$loglik,3) # Log-lik
# Gamma estimates for all 3 model components
# For Gamma dist mean & SD see: https://en.wikipedia.org/wiki/Gamma_distribution
mat[2,1] <- round(fit_G[[crystalSampleIndex]]$lambda[1],3) # lambda
mat[2,2] <- round(fit_G[[crystalSampleIndex]]$gamma.pars[1,1]*
fit_G[[crystalSampleIndex]]$gamma.pars[2,1],3) # mean
mat[2,3] <- round(sqrt(fit_G[[crystalSampleIndex]]$gamma.pars[1,1])*
fit_G[[crystalSampleIndex]]$gamma.pars[2,1],3) # SD
mat[2,4] <- round(fit_G[[crystalSampleIndex]]$lambda[2],3) # lambda
mat[2,5] <- round(fit_G[[crystalSampleIndex]]$gamma.pars[1,2]*
fit_G[[crystalSampleIndex]]$gamma.pars[2,2],3) # mean
mat[2,6] <- round(sqrt(fit_G[[crystalSampleIndex]]$gamma.pars[1,2])*
fit_G[[crystalSampleIndex]]$gamma.pars[2,2],3) # sd
mat[2,7] <- round(fit_G[[crystalSampleIndex]]$lambda[3],3) # lambda
mat[2,8] <- round(fit_G[[crystalSampleIndex]]$gamma.pars[1,3]*
fit_G[[crystalSampleIndex]]$gamma.pars[2,3],3) # mean
mat[2,9] <- round(sqrt(fit_G[[crystalSampleIndex]]$gamma.pars[1,3])*
fit_G[[crystalSampleIndex]]$gamma.pars[2,3],3) # sd
mat[2,10] <- round(fit_G[[crystalSampleIndex]]$loglik,3) # Log-lik
# Normal estimates for all 3 model components
mat[3,1] <- round(fit_LN[[crystalSampleIndex]]$lambda[1],3) # lambda
mat[3,2] <- round(fit_LN[[crystalSampleIndex]]$mu[1],3) # shape
mat[3,3] <- round(fit_LN[[crystalSampleIndex]]$sigma[1],3) # scale
mat[3,4] <- round(fit_LN[[crystalSampleIndex]]$lambda[2],3) # lambda
mat[3,5] <- round(fit_LN[[crystalSampleIndex]]$mu[2],3) # shape
mat[3,6] <- round(fit_LN[[crystalSampleIndex]]$sigma[2],3) # scale
mat[3,7] <- round(fit_LN[[crystalSampleIndex]]$lambda[3],3) # lambda
mat[3,8] <- round(fit_LN[[crystalSampleIndex]]$mu[3],3) # shape
mat[3,9] <- round(fit_LN[[crystalSampleIndex]]$sigma[3],3) # scale
mat[3,10] <- round(fit_LN[[crystalSampleIndex]]$loglik,3) # Log-lik
return(as.data.frame(mat))
}
# render the summary DT tables
library("DT")Below we summarize the mixture-distribution models just for the first two crystallographic features.
2.6.1.9.1 AC1338 Report (Case 1)
2.6.1.9.2 AC1432 Report (Case 2)
2.6.2 2D Kernel Density and 3D Surface Plots
Density estimation is the process of using observed data to compute an estimate of the underlying process’ probability density function. There are several approaches to obtain density estimation, but the most basic technique is to use a rescaled histogram.
Plotting 2D Kernel Density and 3D Surface plots is very important and useful in multivariate exploratory data analytics.
We will use the plot_ly() function in the
plotly package, which works with data frame objects.
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.
To plot the 2D Kernel Density estimation plot we will use the
eruptions data from the “Old Faithful” geyser in Yellowstone National
Park, Wyoming stored under geyser. Also,
kde2d() function is needed for 2D kernel density
estimation.
## [1] 0.8333333 0.9275510 1.0217687 1.1159864 1.2102041## [1] 43.00000 44.32653 45.65306 46.97959 48.30612## [,1] [,2] [,3] [,4] [,5]
## [1,] 9.068691e-13 4.238943e-12 1.839285e-11 7.415672e-11 2.781459e-10
## [2,] 1.814923e-12 8.473636e-12 3.671290e-11 1.477410e-10 5.528260e-10
## [3,] 3.428664e-12 1.599235e-11 6.920273e-11 2.780463e-10 1.038314e-09
## [4,] 6.114498e-12 2.849475e-11 1.231748e-10 4.942437e-10 1.842547e-09
## [5,] 1.029643e-11 4.793481e-11 2.070127e-10 8.297218e-10 3.088867e-09Here z=t(x)%*%y. Then we apply plot_ly to
the list kd using the with() function.
Note we used the option "surface".
For 3D surfaces, we have a built-in dataset in R called
volcano. It records the volcano height at location x, y
(longitude, latitude). Because z is always made from x
and y, we can simply specify z to get the complete
surface plot.
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 100 100 101 101 101 101 101 100 100 100
## [2,] 101 101 102 102 102 102 102 101 101 101
## [3,] 102 102 103 103 103 103 103 102 102 102
## [4,] 103 103 104 104 104 104 104 103 103 103
## [5,] 104 104 105 105 105 105 105 104 104 103
## [6,] 105 105 105 106 106 106 106 105 105 104
## [7,] 105 106 106 107 107 107 107 106 106 105
## [8,] 106 107 107 108 108 108 108 107 107 106
## [9,] 107 108 108 109 109 109 109 108 108 107
## [10,] 108 109 109 110 110 110 110 109 109 108
2.6.3 Multiple 2D image surface plots
#install.packages("jpeg") ## if necessary
library(jpeg)
# Get an image file downloaded (default: MRI_ImageHematoma.jpg)
img_url <- "https://umich.instructure.com/files/1627149/download?download_frd=1"
img_file <- tempfile(); download.file(img_url, img_file, mode="wb")
img <- readJPEG(img_file)
file.info(img_file)## size
## C:\\Users\\dinov\\AppData\\Local\\Temp\\1\\RtmpyMeyyo\\file845c1a016589 8019
## isdir
## C:\\Users\\dinov\\AppData\\Local\\Temp\\1\\RtmpyMeyyo\\file845c1a016589 FALSE
## mode
## C:\\Users\\dinov\\AppData\\Local\\Temp\\1\\RtmpyMeyyo\\file845c1a016589 666
## mtime
## C:\\Users\\dinov\\AppData\\Local\\Temp\\1\\RtmpyMeyyo\\file845c1a016589 2023-06-30 10:14:13
## ctime
## C:\\Users\\dinov\\AppData\\Local\\Temp\\1\\RtmpyMeyyo\\file845c1a016589 2023-06-30 10:14:13
## atime
## C:\\Users\\dinov\\AppData\\Local\\Temp\\1\\RtmpyMeyyo\\file845c1a016589 2023-06-30 10:14:13
## exe
## C:\\Users\\dinov\\AppData\\Local\\Temp\\1\\RtmpyMeyyo\\file845c1a016589 no## [1] TRUEimg <- img[, , 1] # extract the first channel (from RGB intensity spectrum) as a univariate 2D array
# install.packages("spatstat")
# package spatstat has a function blur() that applies a Gaussian blur
library(spatstat)
img_s <- as.matrix(blur(as.im(img), sigma=10)) # the smoothed version of the image
z2 <- img_s + 1 # abs(rnorm(1, 1, 1)) # Upper confidence surface
z3 <- img_s - 1 # abs(rnorm(1, 1, 1)) # Lower confidence limit
# Plot the image surfaces
p <- plot_ly(z=img, type="surface", showscale=FALSE) %>%
add_trace(z=z2, type="surface", showscale=FALSE, opacity=0.98) %>%
add_trace(z=z3, type="surface", showscale=FALSE, opacity=0.98)
p # Plot the mean-surface along with lower and upper confidence services.
2.6.4 3D and 4D Visualizations
Many datasets have intrinsic multi-dimensional characteristics. For instance, the human body is a 3D solid of matter (3 spatial dimensions can be used to describe the position of every component, e.g., sMRI volume) that changes over time (the fourth dimension, e.g., fMRI hypervolumes).
The SOCR BrainViewer shows how to use a web-browser to visualize 2D cross-sections of 3D volumes, display volume-rendering, and show 1D (e.g., 1-manifold curves embedded in 3D) and 2D (e.g., surfaces, shapes) models jointly into the same 3D scene.
We will now illustrate an example of 3D/4D visualization in
R using the packages brainR
and rgl. This
code is included as it runs well in interactive R sessions.
However, it is suppressed during HTML knitting
(eval=FALSE), as rgl causes some browser-OS
combinations to fail while loading the resulting HTML file.
# install.packages("brainR") ## if necessary
library(brainR)
# Test data: https://socr.umich.edu/HTML5/BrainViewer/data/TestBrain.nii.gz
brainURL <- "https://socr.umich.edu/HTML5/BrainViewer/data/TestBrain.nii.gz"
brainFile <- file.path(tempdir(), "TestBrain.nii.gz")
download.file(brainURL, dest=brainFile, quiet=TRUE)
brainVolume <- readNIfTI(brainFile, reorient=FALSE)
brainVolDims <- dim(brainVolume); brainVolDims
# try different levels at which to construct contour surfaces (10 fast)
# lower values yield smoother surfaces # see ?contour3d
contour3d(brainVolume, level = 20, alpha = 0.1, draw = TRUE)
# multiple levels may be used to show multiple shells
# "activations" or surfaces like hyper-intense white matter
# This will take 1-2 minutes to rend!
contour3d(brainVolume, level = c(10, 120), alpha = c(0.3, 0.5),
add = TRUE, color=c("yellow", "red"))
# create text for orientation of right/left
text3d(x=brainVolDims[1]/2, y=brainVolDims[2]/2, z = brainVolDims[3]*0.98, text="Top")
text3d(x=brainVolDims[1]*0.98, y=brainVolDims[2]/2, z = brainVolDims[3]/2, text="Right")
### render this on a webpage and view it!
#browseURL(paste("file://",
# writeWebGL_split(dir= file.path(tempdir(),"webGL"),
# template = system.file("my_template.html", package="brainR"),
# width=500), sep=""))Below we provide some additional 3D/4D PET, sMRI, and fMRI volumes in *.nii.gz format:
- sMRI (3D real-valued structural MRI volume)
- fMRI (4D real-valued functional MRI hyper-volume)
- PET (3D perfusion Positron Emission Tomography volume).
For 4D fMRI time-series, we can load the hypervolumes similarly and then display some lower dimensional projections.
# See examples here: https://cran.r-project.org/web/packages/oro.nifti/vignettes/nifti.pdf
# and here: http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0089470
library(oro.nifti)
fMRIURL <- "https://socr.umich.edu/HTML5/BrainViewer/data/fMRI_FilteredData_4D.nii.gz"
fMRIFile <- file.path(tempdir(), "fMRI_FilteredData_4D.nii.gz")
download.file(fMRIURL, dest=fMRIFile, quiet=TRUE)
(fMRIVolume <- readNIfTI(fMRIFile, reorient=FALSE))## NIfTI-1 format
## Type : nifti
## Data Type : 4 (INT16)
## Bits per Pixel : 16
## Slice Code : 0 (Unknown)
## Intent Code : 0 (None)
## Qform Code : 1 (Scanner_Anat)
## Sform Code : 0 (Unknown)
## Dimension : 64 x 64 x 21 x 180
## Pixel Dimension : 4 x 4 x 6 x 3
## Voxel Units : mm
## Time Units : sec# dimensions: 64 x 64 x 21 x 180 ; 4mm x 4mm x 6mm x 3 sec
fMRIVolDims <- dim(fMRIVolume); fMRIVolDims## [1] 64 64 21 180## [1] 180# Plot the 4D array of imaging data in a 5x5 grid of images
# The first three dimensions are spatial locations of the voxel (volume element) and the fourth dimension is time for this functional MRI (fMRI) acquisition.
image(fMRIVolume, zlim=range(fMRIVolume)*0.95)
h <- hist(fMRIVolume, plot = F)
plot_ly(x = h$mids, y = h$density, type = "bar") %>%
layout(bargap=0.1, title="fMRI Histogram")# Plot an orthographic display of the fMRI data using the axial plane containing the left-and-right thalamus to approximately center the crosshair vertically
orthographic(fMRIVolume, xyz=c(34,29,10), zlim=range(fMRIVolume)*0.9)
stat_fmri_test <- ifelse(fMRIVolume > 15000, fMRIVolume, NA)
h <- hist(stat_fmri_test, plot = F)
plot_ly(x = h$mids, y = h$density, type = "bar") %>%
layout(bargap=0.1, title="fMRI Histogram (high intensities)")## [1] 64 64 21 180
# overlay(fMRIVolume, stat_fmri_test[,,,5], zlim.x=range(fMRIVolume)*0.95)
# To examine the time course of a specific 3D voxel (say the one at x=30, y=30, z=10):
# plot(fMRIVolume[30, 30, 10,], type='l', main="Time Series of 3D Voxel \n (x=30, y=30, z=10)", col="blue")
x1 <- c(1:180)
y1 <- loess(fMRIVolume[30, 30, 10,]~ x1, family = "gaussian")
# lines(x1, smooth(fMRIVolume[30, 30, 10,]), col = "red", lwd = 2)
# lines(ksmooth(x1, fMRIVolume[30, 30, 10,], kernel = "normal", bandwidth = 5), col = "green", lwd = 3)
# legend("bottomright", legend=c("(raw) fMRI", "smooth(fMRI)", "ksmooth(fMRI"),
# col=c("blue", "red", "green"), lty=1, cex=0.8,
# y.intersp=0.8)
plot_ly(x = x1, y = fMRIVolume[30, 30, 10,],
name="Raw fMRI", type = 'scatter', mode = 'lines') %>%
add_trace(y = smooth(fMRIVolume[30, 30, 10,]), name = 'loess fMRI') %>%
add_trace(y = ksmooth(x1, fMRIVolume[30, 30, 10,], kernel="normal", bandwidth = 5)$y, name='kSmooth fMRI') %>%
layout(title="Time Series of 3D Voxel (x=30, y=30, z=10)", legend = list(orientation = 'h'))Chapter 12 provides more details about longitudinal and time-series data analysis.
Finally, DSPA Appendix 3 includes details about classification, representation, modeling, and visualization of parametric and implicit, open and closed manifolds.
3 Appendix
3.1 Importing Data from SQL Databases
We can also import SQL databases into R. First, we need
to install and load the RODBC (R Open Database Connectivity)
package.
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.
3.2 Additional
R scripts
The code below was used to generate some of the 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)3.3 Case-Study 11 - Traumatic Brain Injury (TBI)
The data is available in the Canvas case-studies folder.
# load data CaseStudy11_TBI.xlsx
tmp = tempfile(fileext = ".xlsx")
download.file(url = "https://umich.instructure.com/files/416270/download?download_frd=1", destfile = tmp, mode="wb")
df_TBI <- openxlsx::read.xlsx(xlsxFile = tmp, sheet = "Sheet1", skipEmptyRows = TRUE)
dim(df_TBI)## [1] 46 19Preprocess the data and plot the clustering dendrogram.
# install.packages("dendextend")
library(dendextend)
# Clean the data first (missing values, characters, etc.)
na_strings <- c("NA", ".")
df_TBI_clean <- df_TBI %>% naniar::replace_with_na_all(condition = ~.x %in% na_strings)
df_TBI_clean <- as.data.frame(df_TBI_clean[, -c(3:4)])
df_TBI_clean <- df_TBI_clean %>% tidyr::drop_na ()
dim(df_TBI_clean) # [1] 23 17## [1] 23 17rownames(df_TBI_clean) <- as.character(df_TBI_clean[ ,1])
df_TBI_clean <- df_TBI_clean[, -1]
df_TBI_clean <- as.data.frame(sapply(df_TBI_clean, as.numeric))
df_TBI_clean <- df_TBI_clean[, c("age", "2013.gose", "skull.fx", "temp.injury", "surgery", "acute.sz")]
df_TBI_clean <- as.data.frame(scale(df_TBI_clean))
hc <- hclust(dist(df_TBI_clean), "ave")
dend <- as.dendrogram(hc)
plot_dendro(dend, height = 600) %>%
layout(xaxis = list(range = c(-1, 5))) %>%
hide_legend() %>%
highlight(persistent = TRUE, dynamic = TRUE)# cutree(hc, k = 2)
# alternatively specify the height, which is, the value of the criterion associated with the
# clustering method for the particular agglomeration -- cutree(hc, h= 10)
table(cutree(hc, h= 3)) # cluster distribution##
## 1 2 3 4 5 6
## 6 10 1 3 1 2To identify the number of cases for varying number of clusters
# To identify the number of cases for varying number of clusters we can combine calls to cutree and table
# in a call to sapply -- to see the sizes of the clusters for $2\ge k \ge 10$ cluster-solutions:
# numbClusters=4;
myClusters = sapply(2:5, function(numbClusters)table(cutree(hc, numbClusters)))
names(myClusters) <- paste("Number of Clusters=", 2:5, sep = "")
myClusters## $`Number of Clusters=2`
##
## 1 2
## 19 4
##
## $`Number of Clusters=3`
##
## 1 2 3
## 6 13 4
##
## $`Number of Clusters=4`
##
## 1 2 3 4
## 6 11 4 2
##
## $`Number of Clusters=5`
##
## 1 2 3 4 5
## 6 11 3 1 2Inspect which SubjectIDs are in which clusters:
##
## 1 2
## 19 4groups.k.2 <- cutree(hc, k = 2)
sapply(unique(groups.k.2), function(g) rownames(df_TBI_clean)[groups.k.2 == g])## [[1]]
## [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "11" "12" "14" "15" "16" "17"
## [16] "18" "19" "20" "21"
##
## [[2]]
## [1] "10" "13" "22" "23"Let’s see which Age and which Surgery cohorts fall within each of the derived cluster labels. Remember that all variables are scaled, so they represent standardized variable values!
groups.k.3 <- cutree(hc, k = 3)
sapply(unique(groups.k.3), function(g) df_TBI_clean$age[groups.k.3 == g])## [[1]]
## [1] -0.8625007 0.3227597 -0.4258258 -1.2367934 -1.1744113 0.6346703
##
## [[2]]
## [1] 1.19610942 1.00896305 -1.36155766 -0.80011855 -0.48820793 0.01084907
## [7] 0.13561331 -0.98726492 -0.23867943 2.44375190 1.50802004 0.19799544
## [13] 1.38325579
##
## [[3]]
## [1] -0.1762973 -1.0496470 0.1979954 -0.2386794## [[1]]
## [1] -1.219804 0.784160 -1.219804 -1.219804 -1.219804 0.784160
##
## [[2]]
## [1] 0.784160 0.784160 0.784160 -1.219804 0.784160 0.784160 0.784160
## [8] 0.784160 0.784160 -1.219804 0.784160 0.784160 -1.219804
##
## [[3]]
## [1] -1.219804 0.784160 -1.219804 0.784160# Note that there may be dependencies between some variables
fit <- lm(`2013.gose` ~ age, data = df_TBI_clean)
plot_ly(df_TBI_clean, x = ~age, y = ~`2013.gose`, type = 'scatter', mode = "markers", name="Data") %>%
add_lines(x = ~age, y = fit$fitted.values, mode = "lines", name="Linear Model") %>%
layout(title=paste0("Correlation(2013.gose,age) = ", round(cor(df_TBI_clean$`2013.gose`, df_TBI_clean$age),3)))##
## groups.k.3 -1.21980437173918 0.7841599532609
## 1 4 2
## 2 3 10
## 3 2 2To characterize the clusters, we can look at cluster summary statistics, like the median, of the variables that were used to perform the cluster analysis. These can be broken down by the groups identified by the cluster analysis. The aggregate function will compute stats (e.g., median) on many variables simultaneously. To look at the median values for the variables we’ve used in the cluster analysis, broken up by the cluster groups:
## Group.1 age 2013.gose skull.fx temp.injury surgery acute.sz
## 1 1 -0.6441632 0.7779885 -0.2178222 -1.646252 -1.2198044 -0.448746
## 2 2 0.1356133 -0.1637871 0.7841600 0.581030 0.7841600 -0.448746
## 3 3 -0.2074884 -0.1637871 0.7841600 0.581030 -0.2178222 2.1315443.4 Some additional
ggplot examples
3.4.1 Housing Price Data
This example uses the SOCR Home Price Index data of 19 major city in US from 1991-2009.
library(rvest)
# draw data
wiki_url <- read_html("https://wiki.socr.umich.edu/index.php/SOCR_Data_Dinov_091609_SnP_HomePriceIndex")
hm_price_index<- html_table(html_nodes(wiki_url, "table")[[1]])
head(hm_price_index)## # A tibble: 6 × 23
## Index Year Month `AZ-Phoenix` `CA-LosAngeles` `CA-SanDiego` `CA-SanFrancisco`
## <int> <int> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 1 1991 Janu… 65.3 95.3 83.1 71.2
## 2 2 1991 Febr… 65.3 94.1 81.9 70.3
## 3 3 1991 March 64.6 92.8 80.9 69.6
## 4 4 1991 April 64.4 92.8 80.7 69.5
## 5 5 1991 May 64.4 93.4 81.4 70.1
## 6 6 1991 June 64.9 94.2 82.2 70.8
## # ℹ 16 more variables: `CO-Denver` <dbl>, `DC-Washington` <dbl>,
## # `FL-Miami` <dbl>, `FL-Tampa` <dbl>, `GA-Atlanta` <dbl>, `IL-Chicago` <dbl>,
## # `MA-Boston` <dbl>, `MI-Detroit` <dbl>, `MN-Minneapolis` <dbl>,
## # `NC-Charlotte` <dbl>, `NV-LasVegas` <dbl>, `NY-NewYork` <dbl>,
## # `OH-Cleveland` <dbl>, `OR-Portland` <dbl>, `WA-Seattle` <dbl>,
## # `Composite-10` <dbl>period <- lubridate::parse_date_time(paste(hm_price_index$Year, hm_price_index$Month), "ym")
hm_price_index <- hm_price_index[, c(-1,-2, -3)]
hm_price_index$Date <- period
library(reshape2)
hm_index_melted = melt(hm_price_index, id.vars='Date') #a common trick for plot, wide -> long format
# ggplot(data=hm_index_melted, aes(x=Date, y=value, color=variable)) +
# geom_line(size=1.5) + ggtitle("HomePriceIndex:1991-2009")
plot_ly(hm_index_melted, x=~Date, y=~value, color=~variable,
type="scatter", mode="lines+markers") %>%
layout(title="US Housing Price Index (1991-2009)", yaxis=list(title="HPI"), legend=list(orientation = 'h'))3.4.2 Modeling the home price index data
#Linear regression and predict
hm_price_index$pred = predict(lm(`CA-SanFrancisco` ~ `CA-LosAngeles`, data=hm_price_index))
# ggplot(data=hm_price_index, aes(x = `CA-LosAngeles`)) +
# geom_point(aes(y = `CA-SanFrancisco`)) +
# geom_line(aes(y = pred), color='Magenta', size=2) + ggtitle("PredictHomeIndex SF - LA")
plot_ly(hm_price_index, x=~`CA-LosAngeles`, y=~`CA-SanFrancisco`, color=~`Composite-10`,
type="scatter", mode="lines+markers", name="HPI Data") %>%
add_lines(x = ~`CA-LosAngeles`, y = hm_price_index$pred, mode = "lines", name="Linear Model") %>%
layout(title="LA (SoCal) vs. FS (NoCal)", yaxis=list(title="Los Angeles"),
yaxis=list(title="San Francisco"), legend=list(orientation = 'h'))## Warning: line.color doesn't (yet) support data arrays
## Warning: line.color doesn't (yet) support data arrays
## Warning: line.color doesn't (yet) support data arrays
## Warning: line.color doesn't (yet) support data arraysLet’s examine some popular ggplot graphs.
## # A tibble: 6 × 6
## `IL-Chicago` `MA-Boston` `MI-Detroit` `MN-Minneapolis` `NC-Charlotte`
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 70.0 65.0 58.2 64.2 73.3
## 2 70.5 64.2 57.8 64.2 73.3
## 3 70.6 63.6 57.6 64.2 72.8
## 4 71.1 63.4 57.8 64.3 72.9
## 5 71.4 63.8 58.4 64.8 73.3
## 6 71.7 64.2 58.9 65.0 73.5
## # ℹ 1 more variable: `NV-LasVegas` <dbl>colnames(pairs) <- c("Atlanta", "Chicago", "Boston", "Detroit", "Minneapolis", "Charlotte")
ggpairs(pairs) # you can define the plot design by specifying "upper", "lower", "diag", etc. 
3.4.3 Map of the neighborhoods of Los Angeles (LA)
This example interrogates data of 110 LA neighborhoods, which includes measures of education, income and population demographics.
Here, we select the Longitude and Latitude as the axes, mark these 110 Neighborhoods according to their population, fill out those points according to the income of each area, and label each neighborhood.
library(rvest)
library(ggplot2)
#draw data
wiki_url <- read_html("https://wiki.socr.umich.edu/index.php/SOCR_Data_LA_Neighborhoods_Data")
html_nodes(wiki_url, "#content")## {xml_nodeset (1)}
## [1] <div id="content" class="mw-body" role="main">\n\t\t\t<a id="top"></a>\n\ ...LA_Nbhd_data <- html_table(html_nodes(wiki_url, "table")[[2]])
#display several lines of data
head(LA_Nbhd_data); ## # A tibble: 6 × 15
## LA_Nbhd Income Schools Diversity Age Homes Vets Asian Black Latino White
## <chr> <int> <int> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Adams_Nor… 29606 691 0.6 26 0.26 0.05 0.05 0.25 0.62 0.06
## 2 Arleta 65649 719 0.4 29 0.29 0.07 0.11 0.02 0.72 0.13
## 3 Arlington… 31423 687 0.8 31 0.31 0.05 0.13 0.25 0.57 0.05
## 4 Atwater_V… 53872 762 0.9 34 0.34 0.06 0.2 0.01 0.51 0.22
## 5 Baldwin_H… 37948 656 0.4 36 0.36 0.1 0.05 0.71 0.17 0.03
## 6 Bel-Air 208861 924 0.2 46 0.46 0.13 0.08 0.01 0.05 0.83
## # ℹ 4 more variables: Population <int>, Area <dbl>, Longitude <dbl>,
## # Latitude <dbl>theme_set(theme_grey())
#treat ggplot as a variable
#When claim "data", we can access its column directly e.g., "x = Longitude"
plot1 = ggplot(data=LA_Nbhd_data, aes(x=LA_Nbhd_data$Longitude, y=LA_Nbhd_data$Latitude))
#you can easily add attribute, points, label(e.g., :text)
plot1 + geom_point(aes(size=Population, fill=LA_Nbhd_data$Income), pch=21, stroke=0.2, alpha=0.7, color=2)+
geom_text(aes(label=LA_Nbhd_data$LA_Nbhd), size=1.5, hjust=0.5, vjust=2, check_overlap = T)+
scale_size_area() + scale_fill_distiller(limits=c(range(LA_Nbhd_data$Income)), palette='RdBu', na.value='white', name='Income') +
scale_y_continuous(limits=c(min(LA_Nbhd_data$Latitude), max(LA_Nbhd_data$Latitude))) +
coord_fixed(ratio=1) + ggtitle('LA Neighborhoods Scatter Plot (Location, Population, Income)') 
Observe that some areas (e.g., Beverly Hills) have disproportionately higher incomes and notice that the resulting plot resembles this plot
.
3.4.4 Latin letter frequency in different languages
This example uses ggplot to interrogate the SOCR
Latin letter frequency data.
library(rvest)
wiki_url <- read_html("https://wiki.socr.umich.edu/index.php/SOCR_LetterFrequencyData")
letter<- html_table(html_nodes(wiki_url, "table")[[1]])
summary(letter)## Letter English French German
## Length:27 Min. :0.00000 Min. :0.00000 Min. :0.00000
## Class :character 1st Qu.:0.01000 1st Qu.:0.01000 1st Qu.:0.01000
## Mode :character Median :0.02000 Median :0.03000 Median :0.03000
## Mean :0.03667 Mean :0.03704 Mean :0.03741
## 3rd Qu.:0.06000 3rd Qu.:0.06500 3rd Qu.:0.05500
## Max. :0.13000 Max. :0.15000 Max. :0.17000
## Spanish Portuguese Esperanto Italian
## Min. :0.00000 Min. :0.00000 Min. :0.00000 Min. :0.00000
## 1st Qu.:0.01000 1st Qu.:0.00500 1st Qu.:0.01000 1st Qu.:0.00500
## Median :0.03000 Median :0.03000 Median :0.03000 Median :0.03000
## Mean :0.03815 Mean :0.03778 Mean :0.03704 Mean :0.03815
## 3rd Qu.:0.06000 3rd Qu.:0.05000 3rd Qu.:0.06000 3rd Qu.:0.06000
## Max. :0.14000 Max. :0.15000 Max. :0.12000 Max. :0.12000
## Turkish Swedish Polish Toki_Pona
## Min. :0.00000 Min. :0.00000 Min. :0.00000 Min. :0.00000
## 1st Qu.:0.01000 1st Qu.:0.01000 1st Qu.:0.01500 1st Qu.:0.00000
## Median :0.03000 Median :0.03000 Median :0.03000 Median :0.03000
## Mean :0.03667 Mean :0.03704 Mean :0.03704 Mean :0.03704
## 3rd Qu.:0.05500 3rd Qu.:0.05500 3rd Qu.:0.04500 3rd Qu.:0.05000
## Max. :0.12000 Max. :0.10000 Max. :0.20000 Max. :0.17000
## Dutch Avgerage
## Min. :0.00000 Min. :0.00000
## 1st Qu.:0.01000 1st Qu.:0.01000
## Median :0.02000 Median :0.03000
## Mean :0.03704 Mean :0.03741
## 3rd Qu.:0.06000 3rd Qu.:0.06000
## Max. :0.19000 Max. :0.12000## # A tibble: 6 × 14
## Letter English French German Spanish Portuguese Esperanto Italian Turkish
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 a 0.08 0.08 0.07 0.13 0.15 0.12 0.12 0.12
## 2 b 0.01 0.01 0.02 0.01 0.01 0.01 0.01 0.03
## 3 c 0.03 0.03 0.03 0.05 0.04 0.01 0.05 0.01
## 4 d 0.04 0.04 0.05 0.06 0.05 0.03 0.04 0.05
## 5 e 0.13 0.15 0.17 0.14 0.13 0.09 0.12 0.09
## 6 f 0.02 0.01 0.02 0.01 0.01 0.01 0.01 0
## # ℹ 5 more variables: Swedish <dbl>, Polish <dbl>, Toki_Pona <dbl>,
## # Dutch <dbl>, Avgerage <dbl>## [1] 13.08# require(reshape)
# library(scales)
# dtm = melt(letter[, -14], id.vars = c('Letter'))
# p = ggplot(dtm, aes(x = Letter, y = value, fill = variable)) +
# geom_bar(position = "fill", stat = "identity") +
# scale_y_continuous(labels = percent_format())+ggtitle('Pie Chart')
# #or exchange
# #p = ggplot(dtm, aes(x = variable, y = value, fill = Letter)) + geom_bar(position = "fill", stat = "identity") + scale_y_continuous(labels = percent_format())
# p
# #gg pie plot actually is stack plot + polar coordinate
# p + coord_polar()
reshape2::melt(letter, id.vars='Letter') %>%
plot_ly(x = ~Letter, y = ~value, type = 'bar',
name = ~variable, color = ~variable) %>%
layout(yaxis = list(title = 'Count'), barmode = 'stack')You can see some additional Latin Letters plots here.