Topics discussed in Lecture # 3

Needed R Packages

# Load pckgs for this R session

# -- General 
library(fs)      # file/directory interactions
library(here)    # tools find your project's files, based on working directory
library(paint) # paint data.frames summaries in colour
library(janitor) # tools for examining and cleaning data
library(dplyr)   # {tidyverse} tools for manipulating and summarizing tidy data 
library(forcats) # {tidyverse} tool for handling factors
library(openxlsx) # Read, Write and Edit xlsx Files
library(flextable) # Functions for Tabular Reporting
# -- Statistics
library(rstatix) # Pipe-Friendly Framework for Basic Statistical Tests
library(lmtest) # Testing Linear Regression Models 
library(broom) # Convert Statistical Objects into Tidy Tibbles
#library(tidymodels) # not installed on this machine
library(performance) # Assessment of Regression Models Performance 
# -- Plotting
library(ggplot2) # Create Elegant Data Visualisations Using the Grammar of Graphics

Importing Dataset 1 (NHANES)

# Check my working directory location
# here::here()

# Use `here` in specifying all the subfolders AFTER the working directory 
nhanes_samp <- read.csv(file = here::here("practice", "data_input", "03_datasets",
                                      "nhanes.samp.csv"), 
                          header = TRUE, # 1st line is the name of the variables
                          sep = ",", # which is the field separator character.
                          na.strings = c("?","NA" ), # specific MISSING values  
                          row.names = NULL) 

• Adapting the function here to match your own folder structure

NHANES Variables and their description

Importing Dataset 2 (PREVEND)

# Check my working directory location
# here::here()

# Use `here` in specifying all the subfolders AFTER the working directory 
prevend_samp <- read.csv(file = here::here("practice", "data_input", "03_datasets",
                                      "prevend.samp.csv"), 
                          header = TRUE, # 1st line is the name of the variables
                          sep = ",", # which is the field separator character.
                          na.strings = c("?","NA" ), # specific MISSING values  
                          row.names = NULL) 

PREVEND Variables and their description

Importing Dataset 3 (FAMuSS)

# Check my working directory location
# here::here()

# Use `here` in specifying all the subfolders AFTER the working directory 
famuss <- read.csv(file = here::here("practice", "data_input", "03_datasets",
                                      "famuss.csv"), 
                          header = TRUE, # 1st line is the name of the variables
                          sep = ",", # which is the field separator character.
                          na.strings = c("?","NA" ), # specific MISSING values  
                          row.names = NULL) 

FAMuSS Variables and their description

Explore relationships between two variables

Approaches for summarizing relationships between two variables vary depending on variable types…

• Two numerical variables
• Two categorical variables
• One numerical variable and one categorical variable

Two variables \(x\) and \(y\) are

• positively associated if \(y\) increases as \(x\) increases.
• negatively associated if \(y\) decreases as \(x\) increases.

Two numerical variables (plot)

Height and weight (taken from the nhanes_samp dataset) are positively associated.

• notice we can also use the generic base R function plot for a quick scatter plot

# rename for convenience
nhanes <- nhanes_samp %>% 
  janitor::clean_names()

# basis plot 
plot(nhanes$height, nhanes$weight,
     xlab = "Height (cm)", ylab = "Weight (kg)", cex = 0.8)  

Two numerical variables: correlation (with stats::cor)

is.numeric(nhanes$height) 

[1] TRUE

is.numeric(nhanes$weight)

[1] TRUE

# using `stats` package
stats::cor(x = nhanes$height, y =  nhanes$weight, 
    # argument for dealing with missing values
    use = "pairwise.complete.obs",
    method = "pearson")

[1] 0.4102269

Two numerical variables: correlation (with stats::cor.test)

• Here we compute the Pearson rho (parametric), with the function cor.test (the same we used for testing paired samples)
  • implicitely takes care on NAs

# using `stats` package 
cor_test_result <- cor.test(x = nhanes$height, y =  nhanes$weight, 
                            method = "pearson")

# looking at the cor estimate
cor_test_result[["estimate"]][["cor"]]

[1] 0.4102269

• The function ggpubr::ggscatter gives us all in one (scatter plot + \(r\) (“R”))! 🤯

library("ggpubr") # 'ggplot2' Based Publication Ready Plots
ggpubr::ggscatter(nhanes, x = "height", y = "weight", 
                  cor.coef = TRUE, cor.method = "pearson", #cor.coef.coord = 2,
                  xlab = "Height (in)", ylab = "Weight (lb)")

Spearman rank-order correlation

The Spearman’s rank-order correlation is the nonparametric version of the Pearson correlation.

Spearman’s correlation coefficient, (\(ρ\), also signified by \(rs\)) measures the strength and direction of association between two ranked variables.

• used when 2 variables have a non-linear relationship
• excellent for ordinal data (when Pearson’s is not appropriate), i.e. Likert scale items

To compute it, we simply calculate Pearson’s correlation of the rankings of the raw data (instead of the data).

Spearman rank-order correlation (example)

tabyl(nhanes$education)

 nhanes$education   n percent valid_percent
        8th Grade  32   0.064    0.06412826
   9 - 11th Grade  68   0.136    0.13627255
     College Grad 157   0.314    0.31462926
      High School  94   0.188    0.18837675
     Some College 148   0.296    0.29659319
             <NA>   1   0.002            NA

tabyl(nhanes$health_gen)

 nhanes$health_gen   n percent valid_percent
         Excellent  47   0.094    0.10444444
              Fair  53   0.106    0.11777778
              Good 177   0.354    0.39333333
              Poor  11   0.022    0.02444444
             Vgood 162   0.324    0.36000000
              <NA>  50   0.100            NA

nhanes <- nhanes %>% 
  # reorder education
  mutate (edu_ord = factor (education, 
                            levels = c("8th Grade", "9 - 11th Grade",
                                       "High School", "Some College",
                                       "College Grad" , NA))) %>%  
  # create edu_rank 
  mutate (edu_rank = as.numeric(edu_ord)) %>% 
  # reorder health education
  mutate (health_ord = factor (health_gen, 
                            levels = c( NA, "Poor", "Fair",
                                       "Good", "Vgood",
                                       "Excellent"))) %>%
  # create health_rank 
  mutate (health_rank = as.numeric(health_ord))

Spearman rank-order correlation (example), cont.

table(nhanes$edu_ord, useNA = "ifany" )

     8th Grade 9 - 11th Grade    High School   Some College   College Grad 
            32             68             94            148            157 
          <NA> 
             1 

table(nhanes$edu_rank, useNA = "ifany" )

   1    2    3    4    5 <NA> 
  32   68   94  148  157    1 

table(nhanes$health_ord, useNA = "ifany" )

     Poor      Fair      Good     Vgood Excellent      <NA> 
       11        53       177       162        47        50 

table(nhanes$health_rank,  useNA = "ifany" )

   1    2    3    4    5 <NA> 
  11   53  177  162   47   50 

Spearman rank-order correlation (example cont.)

After setting up the variables in the correct (numerical rank) format, now we can actually compute it: + same function call stats::cor.test + but specifying argument method = "spearman"

# -- using `stats` package 
cor_test_result_sp <- cor.test(x = nhanes$edu_rank,
                               y = nhanes$health_rank, 
                               method = "spearman", 
                               exact = FALSE) # removes the Ties message warning 
# looking at the cor estimate
cor_test_result_sp

    Spearman's rank correlation rho

data:  nhanes$edu_rank and nhanes$health_rank
S = 10641203, p-value = 1.915e-10
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.2946493 

# -- only print Spearman rho 
#cor_test_result_sp[["estimate"]][["rho"]]

Two categorical variables (plot)

mycolors_contrast <- c("#9b2339", "#E7B800","#239b85", "#85239b", "#9b8523","#23399b", "#d8e600", "#0084e6","#399B23",  "#e60066" , "#00d8e6",  "#005ca1", "#e68000")

## genotypes as columns
genotype.race = matrix(table(famuss$actn3.r577x, famuss$race), ncol=3, byrow=T)
colnames(genotype.race) = c("CC", "CT", "TT")
rownames(genotype.race) = c("African Am", "Asian", "Caucasian", "Hispanic", "Other")

# using generic base::barplot
graphics::barplot(genotype.race, col = mycolors_contrast[1:5], ylim=c(0,300), width=2)
legend("topright", inset=c(.05, 0), fill=mycolors_contrast[1:5], 
       legend=rownames(genotype.race))

Two categorical variables (contingency table)

# levels of actn3.r577x
table(famuss$actn3.r577x)

 CC  CT  TT 
173 261 161 

# contingency table to summarize race and actn3.r577x
addmargins(table(famuss$race, famuss$actn3.r577x))

            
              CC  CT  TT Sum
  African Am  16   6   5  27
  Asian       21  18  16  55
  Caucasian  125 216 126 467
  Hispanic     4  10   9  23
  Other        7  11   5  23
  Sum        173 261 161 595

Two categorical variables (contingency table prop)

# adding row proportions
addmargins(prop.table(table(famuss$race, famuss$actn3.r577x), margin =  1))

            
                    CC        CT        TT       Sum
  African Am 0.5925926 0.2222222 0.1851852 1.0000000
  Asian      0.3818182 0.3272727 0.2909091 1.0000000
  Caucasian  0.2676660 0.4625268 0.2698073 1.0000000
  Hispanic   0.1739130 0.4347826 0.3913043 1.0000000
  Other      0.3043478 0.4782609 0.2173913 1.0000000
  Sum        1.7203376 1.9250652 1.3545972 5.0000000

# adding column proportions
addmargins(prop.table(table(famuss$race, famuss$actn3.r577x),margin =  2))

            
                     CC         CT         TT        Sum
  African Am 0.09248555 0.02298851 0.03105590 0.14652996
  Asian      0.12138728 0.06896552 0.09937888 0.28973168
  Caucasian  0.72254335 0.82758621 0.78260870 2.33273826
  Hispanic   0.02312139 0.03831418 0.05590062 0.11733618
  Other      0.04046243 0.04214559 0.03105590 0.11366392
  Sum        1.00000000 1.00000000 1.00000000 3.00000000

Chi Squared test of independence

The Chi-squared test is a hypothesis test used to determine whether there is a relationship between two categorical variables.

• categorical vars. can have nominal or ordinal measurement scale
• the observed frequencies are compared with the expected frequencies and their deviations are examined.

# Chi-squared test
# (Test of association to see if 
# H0: the 2 cat var (race  & actn3.r577x ) are independent
# H1: the 2 cat var are correlated in __some way__

tab <- table(famuss$race, famuss$actn3.r577x)
test_chi <- chisq.test(tab)

the obtained result (test_chi) is a list of objects…

Chi Squared test of independence (cont)

Within test_chi results there are:

# Observed frequencies
test_chi$observed

            
              CC  CT  TT
  African Am  16   6   5
  Asian       21  18  16
  Caucasian  125 216 126
  Hispanic     4  10   9
  Other        7  11   5

# Expected frequencies
round(test_chi$expected  , digits = 1 )

            
                CC    CT    TT
  African Am   7.9  11.8   7.3
  Asian       16.0  24.1  14.9
  Caucasian  135.8 204.9 126.4
  Hispanic     6.7  10.1   6.2
  Other        6.7  10.1   6.2

Chi Squared test of independence (results)

• Recall that:
  • \(H_{0}\): the 2 cat. var. are independent
  • \(H_{1}\): the 2 cat. var. are correlated in some way
• The result of Chi-Square test represents a comparison of the above two tables (observed v. expected):
  • p-value = 0.01286 smaller than α = 0.05 so we REJECT the null hypothesis (i.e. there’s likely an association between race and ACTN3 gene)

test_chi

    Pearson's Chi-squared test

data:  tab
X-squared = 19.4, df = 8, p-value = 0.01286

Computing Cramer’s V after test of independence

Computing Cramer’s V after test of independence (2 ways)

# Compute Creamer's V by hand
 
# inputs 
chi_calc <- test_chi$statistic
n <- nrow(famuss) # N of obd 
n_r <- nrow(test_chi$observed) # number of rows in the contingency table
n_c <- ncol(test_chi$observed) # number of columns in the contingency table

# Cramer’s V
sqrt(chi_calc / (n*min(n_r -1, n_c -1)) )

X-squared 
0.1276816 

# Cramer’s V with rstatix
rstatix::cramer_v(test_chi$observed)

[1] 0.1276816

Chi Squared test of goodness of fit

In some cases the Chi-square test examines whether or not an observed frequency distribution matches an expected theoretical distribution.

Here, we are conducting a type of Chi-square Goodness of Fit Test which:

• serves to test whether the observed distribution of a categorical variable differs from your expectations
• interprets the statistic based on the discrepancies between observed and expected counts

Chi Squared test of goodness of fit (example)

Chi Squared test of goodness of fit (example)

# Subset the vectors of frequencies from the 2 rows  
observed <- c(27,  55,  467, 46)
expected <- c(76.2,  5.95, 478.38,  34.51)

# Calculate Chi-Square statistic manually 
chi_sq_statistic <- sum((observed - expected)^2 / expected) 
df <- length(observed) - 1 
p_value <- 1 - pchisq(chi_sq_statistic, df) 

# Print results 
chi_sq_statistic

[1] 440.2166

df

[1] 3

p_value 

[1] 0

Visualize the data: BMI and age

We are mainly looking for a “vaguely” linear shape here

• ggplot2 gives us a visual confirmation with geom_point()
• Essentially, geom_smooth() adds a trend line over an existing plot
  • inside the function, we have different options with the method argument (default is LOESS (locally estimated scatterplot smoothing))
  • with method = lm we get the linear best fit (the least squares regression line) & its 95% CI

ggplot(nhanes, aes (x = age, 
                          y = bmi)) + 
  geom_point() + 
  geom_smooth(method = lm,  
              #se = FALSE
              )

Linear regression model

The lm() function is used to fit linear models has the following generic structure:

lm(y ~ x, data)

where:

• the 1st argument y ~ x specifies the variables used in the model (here the model regresses a response variable \(y\) against an explanatory variable \(x\).
• The 2nd argument data is used only when the dataframe name is not already specified in the first argument.

Linear regression models syntax

The following example shows fitting a linear model that predicts BMI from age (in years) using data from nhanes adult sample (individuals 21 years of age or older from the NHANES data).

# fitting linear model
lm(nhanes$bmi ~ nhanes$age)

# or equivalently...
lm(bmi ~ age, data = nhanes)

Call:
lm(formula = bmi ~ age, data = nhanes)

Coefficients:
(Intercept)          age  
   28.40113      0.01982  

• Running the function creates an object (of class lm) that contains several components (model coefficients, etc), either directly displayed or accessible with summary() notation or specific functions.

Linear regression models syntax

# name the model object
lr_model <- lm(bmi ~ age, data = nhanes)

# extract model output elements
lr_model$coefficients
lr_model$residuals
lr_model$fitted.values

Linear regression models interpretation: coefficients

• The model tests the null hypothesis \(H_{0}\) that a coefficient is 0
• coefficients outputs are: estimate, std. error, t-statistic, and p-value correspondent to the t-statistic for:
  • intercept
  • explanatory variable(s) slope
• In regression, the population parameter of interest is typically the slope parameter
  • in this model, age doesn’t appear significantly ≠ 0

summary(lr_model)$coefficients 

               Estimate Std. Error   t value      Pr(>|t|)
(Intercept) 28.40112932 0.96172389 29.531480 2.851707e-111
age          0.01981675 0.01824641  1.086063  2.779797e-01

Linear regression models interpretation: Coefficients 2

For the the estimated coefficients of the regression model, we get:

• Estimate = the average increase in the response variable associated with a one unit increase in the predictor variable, (assuming all other predictor variables are held constant).
• Std. Error = a measure of the uncertainty in our estimate of the coefficient.
• t value = the t-statistic for the predictor variable, calculated as (Estimate) / (Std. Error).
• Pr(>|t|) = the p-value that corresponds to the t-statistic. If less than some alpha level (e.g. 0.05). the predictor variable is said to be statistically significant.

Linear regression models outputs: fitted values

Here we see \(\hat{y}_i\), i.e. the fitted \(y\) value for the \(i\)-th individual

fit_val <- lr_model$fitted.values

# print the first 6 elements
head(fit_val)

       1        2        3        4        5        6 
29.39197 29.33252 29.31270 28.95600 29.39197 29.17398 

Linear regression models outputs: residuals

Here we see \(e_i = y_i - \hat{y}_i\), i.e. the residual value for the \(i\)-th individual

resid_val <- lr_model$residuals 

# print the first 6 elements
head(resid_val)

          1           2           3           4           5           6 
-1.49196704  0.06748322 -3.96270002 -3.15599844 -2.49196704  3.75601726 

Linear regression model’s fit: Residual standard error

# Residual Standard error (Like Standard Deviation)

# ---  inputs 
# sample size
n =length(lr_model$residuals)
# n of parameters in the model
k = length(lr_model$coefficients)-1 #Subtract one to ignore intercept
# degrees of freedom of the the residuals 
df_resid = n-k-1
# Squared Sum of Errors
SSE =sum(lr_model$residuals^2) # 22991.19

# --- Residual Standard Error
ResStdErr <- sqrt(SSE/df_resid)  # 6.815192
ResStdErr

[1] 6.815192

Linear regression model’s fit: : \(R^2\) and \(Adj. R^2\)

# --- R^2
summary(lr_model)$r.squared

[1] 0.00237723

# --- Adj. R^2
summary(lr_model)$adj.r.squared

[1] 0.0003618303

Linear regression model’s fit: : F statistic

The F-statistic indicates whether the regression model provides a better fit to the data than a model that contains no independent variables. In essence, it tests if the regression model as a whole is useful.

# extract only F statistic 
summary(lr_model)$fstatistic 

     value      numdf      dendf 
  1.179533   1.000000 495.000000 

# define function to extract overall p-value of model
overall_p <- function(my_model) {
    f <- summary(my_model)$fstatistic
    p <- pf(f[1],f[2],f[3],lower.tail=F)
    attributes(p) <- NULL
    return(p)
}

# extract overall p-value of model
overall_p(lr_model)

[1] 0.2779797

Given the p-value is > 0.05, this indicate that the predictor variable is not useful for predicting the value of the response variable.

Linear regression diagnostic plots: residuals 1/4

ASSUMPTION 1: there exists a linear relationship between the independent variable, x, and the dependent variable, y

For an observation \((x_i, y_i)\), where \(\hat{y}_i\) is the predicted value according to the line \(\hat{y} = b_0 + b_1x\), the residual is the value \(e_i = y_i - \hat{y}_i\)

• A linear (e.g. lr_model) is a particularly good fit for the data when the residual plot shows random scatter above and below the horizontal line.
  • (In this R plot, we look for a red line that is fairly straight)

# residual plot
plot(lr_model, which = 1 )

• We use the argument which in the function plot so we see the plots one at a time.

Linear regression diagnostic plots: normality of residuals 2/4

ASSUMPTION 2: The residuals of the model are normally distributed

With the quantile-quantile plot (Q-Q) we can checking normality of the residuals.

# quantile-quantile plot
plot(lr_model, which = 2 )