BANA 7051 - Study of Wine Quality Part B

Load the Data:

set.seed(7052)
library("tidyverse")

red_wine <- read_delim("https://mailuc-my.sharepoint.com/:x:/g/personal/brancaey_mail_uc_edu/IQBKSlUQ9WdUQZHv4AHIBIXSAc6vGQExmH61pfwWOuLki7I?e=LWB97W&download=1")

Step 1:

1a) Estimate \(\hat\mu\) using the sample mean \(\bar x\)

density_mean <- mean(red_wine$density)
print(density_mean)

[1] 0.9967467

The sample mean of density for the 1,599 wines in Portugal is 0.997, this allows us to infer what the population mean might be.

1b) We can calculate the standard error with \(s\div\sqrt(n)\)

n <- length(red_wine$density)
density_sd <- sd(red_wine$density)
density_error <- density_sd / sqrt(n)
print(density_error)

[1] 4.71981e-05

The standard error of density is very close to zero at approximately 0.0000472.

1c) CLT tells us that we can build a confidence interval using \(\mu \pm 2\sigma\)

upper_bound <- density_mean + (2*density_error)
lower_bound <- density_mean - (2*density_error)
c(lower_bound,upper_bound)

[1] 0.9966523 0.9968411

The 95% confidence interval is (0.99665, 0.99684). This means that we are 95% confident that the true mean of density lies between 0.99665 and 0.99684.

1d) Bootstrap Analysis

set.seed(7052)
density.mu.hat.set <- NULL
for (k in 1:2000) {
  wine.bootstrap <- sample(red_wine$density, size = n, replace = T)
  mu.hat <- mean(wine.bootstrap)
  density.mu.hat.set[k] <- mu.hat
}

sd(density.mu.hat.set)

[1] 4.755423e-05

quantile(density.mu.hat.set, probs = c(0.025, 0.975))

     2.5%     97.5% 
0.9966530 0.9968366 

The bootstrap calculated error is 0.0000476, with a 95% confidence interval of (0.9966530, 0.9968366).

The CLT calculated error is 0.0000472, with a 95% confidence interval of (0.9966523, 0.9968411), if we repeated the resampling process repeatedly, about 95% of the time those intervals would contain the true mean.

This example ended up being so close to the sample mean that the end values are practically identical. This clearly tells us that the bootstrap analysis and the CLT are both strong methods for estimating the standard error of the sample mean.

1e) MLE

Looking at the distribution of density shows that it is generally normally distributed, meaning that a normal distribution would be a good fit for modeling density.

hist(red_wine$density,
     main = "Distribution of Density",
     xlab = "Density",
     ylab = "Frequency",
     col = "red",
     border = "black",
     breaks = 30)

Since this is a normal distribution, we calculated the MLE with a normal distribution.

library(stats4)

density_scaled <- red_wine$density * 10

minuslog.lik <- function(mu, sigma) {
  log.lik <- 0
  for (i in 1:length(red_wine$density)) {
    log.lik <- log.lik+log(dnorm(density_scaled[i],mean=mu,sd=sigma))
  }
  return(-log.lik)
}

est <- mle(minuslog=minuslog.lik,
           start=list(mu=mean(density_scaled),
                      sigma=sd(density_scaled)))

round(coef(summary(est))[, 1:2] / 10, 5)

      Estimate Std. Error
mu     0.99675      5e-05
sigma  0.00189      3e-05

The estimated \(\hat{\mu}\) is 0.997 with a standard error of 0.00005 and a \(\hat{\sigma}\) value of 0.002 with a standard error of 0.00003.

Step 2:

2a) Estimate \(\mu\) using the sample mean for residual sugar.

sugar_mean <- mean(red_wine$`residual sugar`)
print(sugar_mean)

[1] 2.538806

The sample mean \(\bar x\) = 2.5388

2b) Can CLT be used to quantify the variability? can CLT be used to obtain a 95% confidence interval?

Yes, We can use the CLT to quantify the variability of the estimate and create a confidence interval. The distribution of the mean will approximately follow a normal distribution regardless of the skew. We use this approximately normally distributed variability of the sample mean to build the confidence interval.

n <- length(red_wine$`residual sugar`) 
sugar_sd <- sd(red_wine$`residual sugar`) 
sugar_error <- sugar_sd / sqrt(n) 
print(sugar_error) 

[1] 0.03525922

The standard error of the sample mean is ~0.0353.

sugar_mean <- mean(red_wine$`residual sugar`)  
upper_bound <- sugar_mean + (2*sugar_error) 
lower_bound <- sugar_mean - (2*sugar_error) 
c(lower_bound,upper_bound)

[1] 2.468287 2.609324

The 95% confidence interval is ~(2.468, 2.609). We are 95% confident the true population mean for residual sugar falls within this interval.

2c) Bootstrap Method

set.seed(7052)
sugar.mu.hat.set <- NULL
for (k in 1:2000) {
  wine.bootstrap <- sample(red_wine$`residual sugar`, size = n, replace = T)
  mu.hat <- mean(wine.bootstrap)
  sugar.mu.hat.set[k] <- mu.hat
}

sd(sugar.mu.hat.set)

[1] 0.03553346

CI <- quantile(sugar.mu.hat.set, probs = c(0.025, 0.975))
CI

    2.5%    97.5% 
2.471567 2.609863 

The estimated variability is ~0.0355 with a confidence interval of ~(2.472, 2.610)

hist(sugar.mu.hat.set, freq=FALSE)
lines(density(sugar.mu.hat.set), lwd=5, col = 'blue')

The density distribution is symmetric and is normally distributed. If we compare how far each quantile is from the mean, we can also see exactly how symmetric it is.

lower_bound <- sugar_mean - CI[1]
upper_bound <- CI[2] - sugar_mean

c(lower_bound, upper_bound)

      2.5%      97.5% 
0.06723812 0.07105769 

We can see that they are not perfectly symmetric, but the standard error will trend toward a normal distribution consistent with the Central Limit Theorem.

2d) MLE

We see that the distribution of residual sugar is highly right skewed with most values clustered around 0-5.

hist(red_wine$`residual sugar`,
     main = "Distribution of Residual Sugar",
     xlab = "Residual Sugar",
     ylab = "Frequency",
     col = "red",
     border = "black",
     breaks = 30)

This means that a different distribution is needed for a better fit of the data. A log normal distribution should be used instead of a normal distribution because it is able to match the shape of the skew and only includes values 0 and greater which works for residual sugar.

The log normal distribution has 2 parameters, \(\mu\) and \(\sigma\).

library(stats4)

minuslog.lik <- function(mu, sigma) {
  log.lik <- 0
  for (i in 1:length(red_wine$`residual sugar`)) {
    log.lik <- log.lik + log(dlnorm(x = red_wine$`residual sugar`[i], meanlog = mu, sdlog = sigma))
  }
  return(-log.lik)
}

est.lognorm <- mle(minuslog = minuslog.lik, 
                   start = list(mu = log(mean(red_wine$`residual sugar`)), 
                                sigma = log(sd(red_wine$`residual sugar`))))

round(coef(summary(est.lognorm))[, 1:2],3)

      Estimate Std. Error
mu       0.850      0.009
sigma    0.357      0.006

The MLE estimate of these parameters is 0.85 for \(\hat\mu\) with a standard error of 0.009 and 0.357 for \(\hat\sigma\) with a standard error of 0.006.

Step 3:

red_wine$excellent <- ifelse(red_wine$quality >= 7, 1, 0)

3a) Confidence Interval

excellent_mean <- mean(red_wine$excellent)  
n <- length(red_wine$excellent) 
excellent_sd <- sd(red_wine$excellent) 
excellent_error <- excellent_sd / sqrt(n) 
upper_bound <- excellent_mean + (2*excellent_error) 
lower_bound <- excellent_mean - (2*excellent_error) 
c(lower_bound,upper_bound)

[1] 0.1185751 0.1528445

The 95% confidence interval for the new variable “excellent” is ~(0.1186, 0.1529). This means that we are 95% confident that the true population proportion of excellent wines in Northern Portugal lies with the interval of 11.86% to 15.29%.

3b) Bootstrap Method

set.seed(7052)
excellent.mu.hat.set <- NULL
for (k in 1:2000) {
  wine.bootstrap <- sample(red_wine$excellent, size = n, replace = T)
  mu.hat <- mean(wine.bootstrap)
  excellent.mu.hat.set[k] <- mu.hat
}

sd(excellent.mu.hat.set)

[1] 0.008559933

quantile(excellent.mu.hat.set, probs = c(0.025, 0.975))

     2.5%     97.5% 
0.1181989 0.1532208 

The bootstrap estimated 95% confidence interval for \(\hat{p}\) is (0.1182, 0.1532)

3c) Compare the Intervals

The intervals are off by only (0.0004, 0.0003), which could be attributed to random error, as standard error will tend towards a normal distribution under the CLT. This again shows that the bootstrap and CLT methods are comparable when creating confidence intervals for estimating the mean of population p, even though this variable is binary. This is probably because the sample is so large.

3d) MLE and Standard Error

Since this is a binary variable, we use a binomial distribution to calculate the MLE.

library(stats4)

minuslog.lik <- function(p) {
  log.lik <- 0
  for (i in 1:length(red_wine$excellent)) {
    log.lik <- log.lik + log(dbinom(x = red_wine$excellent[i], size = 1, prob = p))
  }
  return(-log.lik)
}

est.p <- mle(minuslog = minuslog.lik, start = list(p = 0.5))
round(coef(summary(est.p))[, 1:2],3)

  Estimate Std. Error 
     0.136      0.009 

The estimated MLE of \(\hat{p}\) is 0.136 with a standard error of 0.009. This means that we estimate about 13.6% of wines from northern Portugal would be rated excellent in quality with a range of possible estimated values within 1.8% overall, this is a relatively precise estimate.