R Workshop

Module 5: Statistical modeling with R (1)

2018-04-18
Bobae Kang
(Bobae.Kang@illinois.gov)

Source: Vadlo.com

summary(data)

• Base R's summary() function is a quick way to get descriptive statistics on each columan of a tabular data object.
  • For numerical columns, we get minimum, 1st quartile, median, mean, 3rd quartie and maximum values, as well as the count of missing value (NA).
  • For categorical variables (e.g. factor), we can the count of each level as well as the missing value (NA).
• summary() is also used to get the “summary” of a fitted model object (e.g. lm) as well, as we will see shortly.

• Central tendency
• Variability
• Shape
• Outliers

• In statistics, central tendency is concerned with the “center” of a distribution. The following three are the common measures of central tendency:
  • “mean” for the arithmetic mean
  • “median” for the 50th percentile
  • “mode” for the most frequent value
• R has built-in functions for the first two: mean() and median()

• Variability (or dispersion) in statistics is concerned with the extent to which a distribution is spreaded out.
• There exist various measures for variability and R has built-in functions for many of them:
  • range: min(), max(), range()
  • percentiles: fivenum(), IQR(), quantile()
  • variance: var(), sd()

• The skewness and kurtosis (fatness/thinness) of a distribution are called the “shape” of the distribution.
• The moments package offers the following functions to meaasure the shape of a distribution:
  • skewness()
  • kurtosis()
  • See moments package documentation

• Outliers are observations or data points that are far from most other observations and disproportinately affect key summary statistics
• The outliers package offers useful functions to detect and measure outliers in data.
  • See outliers package documentation

• Frequency table is one of the most common ways to summarize the distribution of a catagorical variable
• Frequency tables can be made using table functions
  • table() for generating frequency tables
  • prop.table() for tables of proportions
  • xtabs() for creating frequency tables using formula
  • ftable() for creating “flat” contingency tables

table(...)
prop.table(x, margin = NULL)
ftable(x)
xtabs(formula, data, ...)

• table() takes one or more data vectors of same length
  • Each input data can be named
  • Use as.data.frame() to turn a table into a data frame
• prop.table() and ftable() takes a table object
• xtabs use formula to generate a frequency table
  • If a data frame is provided as the data input, its column names can be used directly in formula

Source: Vadlo.com

t.test(x, y = NULL, alternative = c("two.sided", "less", "greater"),
       mu = 0, conf.level = 0.95, ...)
t.test(formula, data, subset, na.action, ...)

• The null hypothesis assumes normal distribution
• The tested null hypothesis is no difference in mean and the alternative hypothesis can one of the three options, with the default of “two.sided”
• One sample t-test is conducted with only x input is given
  • The mean of x is compared against the mu value (default 0)
• Two sample t-test can be conducted by either supplying x and y inputs or using formula (y ~ x) with data input

wilcox.test(x, y = NULL, alternative = c("two.sided", "less", "greater"),
            mu = 0, conf.level = 0.95, ...)
wilcox.test(formula, data, subset, na.action, ...)

• The null hypothesis does NOT assume normality
• One sample Wilcoxon test (rank-sum test) is conducted with only x input is provided
• Two sample Wilcoxon test (singed-rank test) can be conducted by either supplying x and y inputs or using formula (y ~ x) with data input

aov(formula, data, qr = TRUE, ...)

• formula is of the following format: value ~ subgroup
• summary() has a method for aov class that returns a refined print result
• model.tables() has a method for aov class for reporting table of means or table of effects

# Test of equal proprtions
prop.test(x, n, p, ...)

# Chi-square test of indepndence/goodness-of-fit
chisq.test(x, y = NULL, ...)

# Shapiro-Wilk test of normality
shapiro.test(x)

# One- or two-sample Kolmogorov-Smirnov test
ks.test(x, y, ..., alternative = "two.sided", exact = NULL)

# F-test to compare variances
var.test(x, y, ratio = 1, alternative = "two.sided",
         conf.level = 0.95, ...)
var.test(formula, data, subset, na.action, ...)

# Test of correlation between paired samples
cor.test(x, y, alternative = "two.sided", conf.level = 0.95,
         method = c("pearson", "kendell", "spearman"), ...)

\[ y_i = \beta_0 + \beta_1 x_i + \varepsilon_i \]
\[ \boldsymbol{\text{y}} = \boldsymbol{\text{X}}^{\text{T}}\boldsymbol{\beta} + \boldsymbol{\varepsilon} \]

lm(formula, data, subset, na.action, ...)

• lm() is used to fit linear models, according to the formula input.
• data is an optional argument, which can take a data frame
  • If data input is provided, its column names can be used directly in the formula

formula <- y ~ x1 + x2

• The response variable is placed on the left-hand side of the tilde symbol (~), and the explanatory variables on the right-hand side
  • Adding multiple explanatory variables can be done using + symbol
• Simply use y ~ . to include all other columns additively as explanatory variables

lmfit <- lm(formula, data)

print(lmfit) # not recommended
summary(lmfit)

• lm() output is an object of lm class
• Simply printing an lm object gives only minimal information on the fitted model
• summary() function with a model input (e.g. lm) prints a more comprehensive model summary

plot(lmfit, which = c(1:3, 5), ...)

• Base R plot() function has a method for lm objects
• Total 6 plots are generated
  • “Residuals vs Fitted”
  • “Normal Q-Q”
  • “Scale-Location”
  • “Cook's distance”
  • “Residuals vs Leverage”
  • “Cook's dist vs Leverage \( h_{ii}/(1-h_ii) \)”

logLik(object, ...)
AIC(object, ..., k = 2)
BIC(object, ...) # equivalent to AIC with k = log(n)

• There are measures of model fit other than \( R^2 \).
• For a model fitted by maximum likelihood algorithm:
  • logLik() for log-liklehood value
  • AIC() for Akaike Information Criterion (with default k)
  • BIC() for Bayesian Information Criterion
• Read this Wikipedia article for more on maximum likelihood estimation

• Transformations
• Interaction terms
• Polynomials
• Polynomial splines

scale(y) ~ log(x1) + sqrt(x2) + ...

• We can transform formula terms on either side of ~ by applying a function to each term
• Common transformations include:
  • log(): natural log (cf. log2 and log10)
  • exp(): exponentiate
  • sqrt(): square root
  • scale(): rescale data such that the mean is 0 and the standard deviation is 1

y ~ x1 + x2 + x1:x2
y ~ x1 * x2

• x1:x2 takes the interactions of all terms in x1 and x2
• x1*x2 is equivalent to x1:x2 and the original terms x1 and x2
• It is recommended to use an interaction term with all the original terms

y ~ ploy(x, degree = n)
y ~ x + I(x^2) + I(x^3) + ...

• poly() generates the n-th degree polynimals of the input x
• ploy(x, 3) is equivalent to using x + I(x^2) + I(x^3)
  • I() ensures the ^ with a symbol is used as an arithmetic operator

library(splines) # part of R "base pacakges"
y ~ ns(x, ...)

y ~ ns(x, ...)

y ~ bs(x, degree = n, ...)

y ~ bs(x, degree = n, ...)

• ns() is used for the natural cubic splines
• bs() is used for the n-th degree polynomial splines
  • equivalent to bs(x, degree = 3)

• broom is a tidyverse package for “tidying up” statistical model outputs in R, i.e., converting them into tidy data frames.
• broom functions:
  • glance()
  • tidy()
  • augment()
• Visit broom GitHub repository and “broom as well as dplyr” vignette for more

glance(x, ...)
tidy(x, ...)
augment(x, ...)

• glance() takes a model object and returns a concise one-row summary of the model.
• tidy() takes a model object and returns a data.frame with the coeffcient estimation results
• augment() takes a model object and returns a data.frame for each observation with additional columns such as predictions, residuals and cluster assignments

• modelr is a tidyverse package that provides helper functions for statistical modeling in R
  • partitioning and sampling
  • model quality metrics
  • interacting with models
• Visit modelr GitHub repository and the “Model Basics” chapter in R for Data Science for more

resample(data, idx)
resample_partition(data, p)
bootstrap(data, n, id = ".id")

• resample take a data frame and a vector of indices to create a pointer object, which can be either turned into a data frame or used directly in modeling functions (e.g. lm())
• resample_partition takes a data and a named numeric vector to specify partitioning. Its output is a list of random partitions.
• bootstrap takes a data frame and an integer to specify the number of bootstrap replicates to generate. The output is a data frame of bootstrap samples column and id column

rmse(model, data)
mae(model, data)
qae(model, data, probs = c(0.05, 0.25, 0.5, 0.75, 0.95))
rsquare(model, data)

• rmse() returns the root mean squared error
• mae() returns the mean abosolute error
• qae() returns quantiles of absolute error
• rsquare() returns the \( R^2 \) value, i.e. the variance of the predictions divided by the variance of the reponse

add_predictions(data, model, var = "pred")
add_residuals(data, model, var = "resid")

• add_predictions() adds a new column to a data frame for predicted values based on a model
• add_residuals() adds a new column to a data frame for residuals based on a model

\[ \text{E}[\boldsymbol{\text{Y}}] = \mu = g^{-1}(\boldsymbol{\text{X}\beta}) \]
\[ \text{Var}[\boldsymbol{\text{Y}}] = \text{V}[\mu] = \text{V}[g^{-1}(\boldsymbol{\text{X}}\boldsymbol{\beta})] \]

• Generalized linear models (GLMs) extend the linear model to fit response variables with error distribution other than the Gaussian (normal) distribution
• Three components of GLM:
  • A probability distribution from the exponential family
  • A linear predictor, \( \eta = \boldsymbol{\text{X}\beta} \)
  • A link function \( g \), such that \( \text{E}[\boldsymbol{\text{Y}}] = \mu = g^{-1}(\eta) \)
• See this Wikipedia article for more on GLM and the model components

glm(formula, family = gaussian, data, ...)

• formula is the same formula we have seen all along
• family input specifies the link function \( g \)
  • the default is gaussian(link = "identity") for the linear model
• glm() function returns a glm class object
  • Use summary() to inspect the model estimation results

glm(formula, family = binomial, data, ...)

• A logit model, i.e. logistic regression, can be fitted using glm() with family = binomial
• The response variable must range between 0 and 1 (an error will be thrown otherwise)

glm(formula, family = poisson, data, ...)

• A possion model can be fitted using glm() with family = poisson
• The response variable must not include negative values (an error will be thrown otherwise)

There are thrid-party packages to provide additional fuctions and/or different interfaces to existing statistical analysis functions. Although I have little exposure to them, the following two might be worth checking:

• psych package
  • See this package vignette by the package author
• jmv package
  • See this package documentation page

• Kabacoff, R. “Statistics”. Quick-R.
• Lane D. et al. Online Statistics: A Multimedia Course of Study.
• Prabhakaran, S. r-statistics.co.
• UCLA Institute for Digital Research and Education. “R”.
• University of Cincinnati. “Descriptive analytics”. UC Business Analytic R Programming Guide.
• University of Cincinnati. “Predictive analytics”. UC Business Analytic R Programming Guide.
• Wollschlaeger, D. RExRepos.
• Yau, C. “Elementary Statistics with R”. R Tutorial.

Source: Giphy