Introduction: Random Variables - Cats example

Loading the data from the MASS package

# Getting the data
cats <- MASS::cats

Examining the data:

Let’s start by looking at the data using slice_sample(data). Why slice_sample() instead of another function like head()? It’s not uncommon for the data sets available in R or online to have the data sorted in some way, so the first 6 rows might all be similar. In fact, the cats are sorted by sex, then Bwt, then Hwt.

# first six sorted rows
head(cats)

##   Sex Bwt Hwt
## 1   F 2.0 7.0
## 2   F 2.0 7.4
## 3   F 2.0 9.5
## 4   F 2.1 7.2
## 5   F 2.1 7.3
## 6   F 2.1 7.6

# six random rows
slice_sample(cats, n = 6)

##   Sex Bwt  Hwt
## 1   F 2.1  7.6
## 2   F 2.9 10.1
## 3   M 2.5 11.0
## 4   M 3.6 15.0
## 5   M 2.9 11.8
## 6   M 3.2 12.3

We can use summary(data) to summarize each column individually, including see the number of missing values by column (if any).

# Summary of each column
summary(cats)

##  Sex         Bwt             Hwt       
##  F:47   Min.   :2.000   Min.   : 6.30  
##  M:97   1st Qu.:2.300   1st Qu.: 8.95  
##         Median :2.700   Median :10.10  
##         Mean   :2.724   Mean   :10.63  
##         3rd Qu.:3.025   3rd Qu.:12.12  
##         Max.   :3.900   Max.   :20.50

There are 47 female cats (‘F’) and 97 male cats (‘M’). summary() also gives the mean and 5 number summary for each numeric (quantitative) column.

A better way is to visualize the data with either a histogram or density plot. I prefer density plots since there isn’t a need to specify the bin width or bin number, but either one works:

# Stacking both Bwt and Hwt into one column called cats_long
cats_long <- 
  cats |>
  dplyr::select(
    sex = Sex,
    `body weight (kg)` = Bwt,
    `heart weight (g)` = Hwt
  ) |>
  # adding an id column
  mutate(
    .before = 1,
    cat_id = row_number()
  ) |>
  pivot_longer(
    cols = c(-sex, -cat_id),
    names_to = 'feature',
    values_to = 'weight'
  )

# Creating the density plots
  ggplot(
    data = cats_long,
    mapping = aes(x = weight)
  ) +
  geom_density(
    fill = 'orange'
  ) + 
  facet_wrap(
    facets = vars(feature),
    nrow = 2,
    scales = 'free'
  ) + 
  labs(
    title = 'Body and heart weight of 144 cats',
    x = NULL
  ) + 
  theme_bw() +
  # Making the density plot 'sit' on the x-axis
  scale_y_continuous(
    expand = c(0, 0, 0.05, 0)
  )

Both heart and body weight of the cats are unimodal and right skewed. Body weight has a min of 2 kg, max of about 4 kg, and a mode of around 2.3 kgs

Heart weight has a min of about 6 g, max of about 10 g, and a mode of 10 g.

Summary stats

When using probability distributions, we often need to estimate the parameters. For numeric data, common probability distributions used have a mean parameter, \(\mu\), and a variance parameter, \(\sigma^2\), or standard deviation, \(\sigma\).

We estimate the parameters from the data using sample statistics:

\[\hat{\mu} = \bar{y} = \frac{1}{n}\sum_{i = 1}^n y_i\]

\[\hat{\sigma}^2 = s^2 = \frac{1}{n-1}\sum_{i = 1}^n (y_i-\bar{y})^2\]

\[\hat{\sigma} = s = \sqrt{\frac{1}{n-1}\sum_{i = 1}^n (y_i-\bar{y})^2}\]

To speed up the calculating the mean and standard deviation, we’ll use the cats_long data set created in the previous code chunk:

slice_sample(cats_long, n = 10)

## # A tibble: 10 × 4
##    cat_id sex   feature          weight
##     <int> <fct> <chr>             <dbl>
##  1    135 M     heart weight (g)   17.2
##  2     34 F     heart weight (g)   10.2
##  3      9 F     heart weight (g)    8.3
##  4     56 M     body weight (kg)    2.2
##  5     24 F     body weight (kg)    2.3
##  6     45 F     heart weight (g)   10.1
##  7     30 F     body weight (kg)    2.3
##  8    132 M     body weight (kg)    3.5
##  9     84 M     body weight (kg)    2.7
## 10    114 M     heart weight (g)   14.3

Now we’ll calculate the parameter estimates using the summarize() function in the dplyr package:

cats_long |>
  summarize(
    .by = feature,            # separate parameters for each RV
    avg = sum(weight) / n(),  # n() counts the number of rows used in the subset
    avg2 = mean(weight),      # or we can use the mean() function
    
    # Calculating the variability parameters
    variance = (sum((weight - mean(weight))^2)/(n() - 1)),
    variance2 = var(weight),
    stdev = sqrt(variance),
    stdev2 = sd(weight)
  )

## # A tibble: 2 × 7
##   feature            avg  avg2 variance variance2 stdev stdev2
##   <chr>            <dbl> <dbl>    <dbl>     <dbl> <dbl>  <dbl>
## 1 body weight (kg)  2.72  2.72    0.236     0.236 0.485  0.485
## 2 heart weight (g) 10.6  10.6     5.93      5.93  2.43   2.43

We can estimate the covariance and correlation:

\[\widehat{cov}(X,Y) = \hat{\sigma}_{xy} = \frac{1}{n-1} \sum_{i = 1}^n(x_i - \bar{x})(y_i -\bar{y}) \]

\[r = \frac{\hat{\sigma}_{xy}}{s_x s_y} = \frac{\sum_{i = 1}^n(x_i - \bar{x})(y_i -\bar{y})}{\sqrt{\sum_{i = 1}^n(x_i - \bar{x})^2}\sqrt{\sum_{i = 1}^n(y_i - \bar{y})^2}}\]

cats |>
  summarize(
    # Calculating the covariance by "hand" and using the R function
    covariance  = sum((Bwt - mean(Bwt)) * (Hwt - mean(Hwt)))/(n() - 1),
    covariance2 = cov(Bwt, Hwt),
    
    # Calculating the correlation 'by hand' and using the R function
    correlation  = sum((Bwt - mean(Bwt)) * (Hwt - mean(Hwt))) / 
                   sqrt(sum((Bwt - mean(Bwt))^2) * sum((Hwt - mean(Hwt))^2)),
    correlation2 = covariance / (sd(Bwt) * sd(Hwt)),
    correlation3 = cor(Bwt, Hwt)
  )

##   covariance covariance2 correlation correlation2 correlation3
## 1  0.9501127   0.9501127   0.8041274    0.8041274    0.8041274