Preface

For all analysis, a significance threshold of $\alpha = 0.05$ is employed.

Problem One

Part A: In Vitro Experiment

Task 1.1

Before anything else, I need to load in and visualise the data so that I can get a handle on everything. Then I can make some minor modifications and/or transformations depending on my needs. For this data set I wanted to be able to visualise the mean and and standard error in my plots so I used dplyr to generate a tibble containing these values.

vitro_data <- read.csv("in_vitro_felodipine.csv")
vitro_data

##    treatment replicate metabolised_felodipine
## 1          N         1                   64.9
## 2          N         2                   75.0
## 3          N         3                   62.6
## 4          N         4                   59.4
## 5          N         5                   69.1
## 6          C         1                   33.0
## 7          C         2                   44.4
## 8          C         3                   33.3
## 9          C         4                   46.6
## 10         C         5                   43.1
## 11         G         1                   20.4
## 12         G         2                    9.5
## 13         G         3                   17.7
## 14         G         4                   21.7
## 15         G         5                   41.0
## 16         B         1                   17.6
## 17         B         2                   16.1
## 18         B         3                   10.3
## 19         B         4                    9.5
## 20         B         5                   15.3

avg_vitro <- vitro_data %>%
  group_by(treatment) %>%
  summarise(mean_se(metabolised_felodipine))
avg_vitro

## # A tibble: 4 × 4
##   treatment     y  ymin  ymax
##   <chr>     <dbl> <dbl> <dbl>
## 1 B          13.8  12.1  15.4
## 2 C          40.1  37.2  43.0
## 3 G          22.1  16.9  27.2
## 4 N          66.2  63.5  68.9

str(vitro_data)

## 'data.frame':    20 obs. of  3 variables:
##  $ treatment             : chr  "N" "N" "N" "N" ...
##  $ replicate             : int  1 2 3 4 5 1 2 3 4 5 ...
##  $ metabolised_felodipine: num  64.9 75 62.6 59.4 69.1 33 44.4 33.3 46.6 43.1 ...

The Design:

treatment
- Factor with four levels
replicate
- Numeric, integer
metabolised_felodipine
- Numeric, double

My first instinct with these data is to employ geom_point() so that the raw data isn’t hidden.

ggplot(vitro_data, 
       aes(treatment, metabolised_felodipine)
       ) +
  theme_bw() +
  geom_point(data = avg_vitro, 
             aes(x = treatment, y = y, colour = treatment, fill = treatment),
             shape = 23,
             size = 3,
             show.legend = FALSE,
             na.rm = F) +
  geom_errorbar(data = avg_vitro, 
                aes(x = treatment, y = y, ymin = ymin, ymax = ymax, colour = treatment),
                width = 0.1,
                show.legend = FALSE) +
  ylim(0, max(vitro_data$metabolised_felodipine)) +
  geom_jitter(width = 0.2) +
  labs(x = "Treatment", y = "Metbolised Felodipine (ng/ml)",
       title ="Per Treatment Metabolysis of Felodipine by CYP3A4",
       subtitle = "With Diamands Representing Mean Metabolised Felodipine")

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

I like this plot but I think it could be better. I want to order the x-axis by mean from lowest to highest so that it’s easier to compare between group differences by eye. I’ll also add a geom_violin() layer so that it’s easier to understand the distribution of values within each group. However, given the small size of the data set, I think it’s important be mindful when attributing meaning to the shape of the violins.

ggplot(vitro_data, 
       aes(x = forcats::fct_reorder(treatment, metabolised_felodipine, .fun = mean), metabolised_felodipine)) +
  geom_violin(aes(colour = treatment, alpha = 0.5),
              show.legend = FALSE,
              trim = FALSE) +
    geom_jitter(width = 0.2) +
    geom_point(data = avg_vitro, 
             aes(x = treatment, y = y, colour = treatment, fill = treatment),
             shape = 23,
             size = 3,
             show.legend = FALSE) +
  geom_errorbar(data = avg_vitro, 
                aes(x = treatment, y = y, ymin = ymin, ymax = ymax, colour = treatment),
                width = 0.1,
                show.legend = FALSE) +
  theme_bw() +
  labs(x = "Treatment",
       y = "Metbolised Felodipine (ng/ml)",
       title ="Per Treatment Metabolysis of Felodipine by CYP3A4",
       subtitle = "With Diamands Representing Mean Metabolised Felodipine")

This is nice! A quick inspection seems to indicate that both treatments are able to reduce felodipine metabolism and that the combined treatment could be more effective than either treatment on their own. This is a good place to start running a proper analysis of the data and generate some helpful statistics.

lm_vitro <- lm(metabolised_felodipine ~ treatment, vitro_data)
summary(lm_vitro)

## 
## Call:
## lm(formula = metabolised_felodipine ~ treatment, data = vitro_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.560  -4.285  -0.830   3.225  18.940 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   13.760      3.362   4.093 0.000849 ***
## treatmentC    26.320      4.755   5.536 4.52e-05 ***
## treatmentG     8.300      4.755   1.746 0.100030    
## treatmentN    52.440      4.755  11.030 6.91e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.518 on 16 degrees of freedom
## Multiple R-squared:  0.8994, Adjusted R-squared:  0.8805 
## F-statistic: 47.68 on 3 and 16 DF,  p-value: 3.343e-08

par(mfrow = c(2,2))
plot(lm_vitro)

Task 1.2

In this section we are going to compare the mean extent to which CYP3A4 metabolises felodipine after exposure to coffee or grapefruit juice. Then we look at the effect of exposure to both treatment conditions at the same time. A negative control was also established by measuring felodipine metabolism in the absence of either treatment.

Although the data is presented such that there appears to be a single factor, “treatment”, with three levels; “coffee”, “grapefruit”, and “both”. We can actually separate out the coffee and grapefruit components from “both” and package the data under two factors; “grapefruit”, and “coffee”. Each with two levels; yes (“Y”), and no (“N”).

By separating the variables in this way, our analyses more accurately reflect the reality of the experiment.

split_treat <- function(data)
{
  return(data %>% mutate(grapefruit = if_else(treatment %in% c("G","B"),"Y","N"),coffee = if_else(treatment %in% c("C","B"),"Y","N")))
}
split_data <- split_treat(vitro_data)
split_data

##    treatment replicate metabolised_felodipine grapefruit coffee
## 1          N         1                   64.9          N      N
## 2          N         2                   75.0          N      N
## 3          N         3                   62.6          N      N
## 4          N         4                   59.4          N      N
## 5          N         5                   69.1          N      N
## 6          C         1                   33.0          N      Y
## 7          C         2                   44.4          N      Y
## 8          C         3                   33.3          N      Y
## 9          C         4                   46.6          N      Y
## 10         C         5                   43.1          N      Y
## 11         G         1                   20.4          Y      N
## 12         G         2                    9.5          Y      N
## 13         G         3                   17.7          Y      N
## 14         G         4                   21.7          Y      N
## 15         G         5                   41.0          Y      N
## 16         B         1                   17.6          Y      Y
## 17         B         2                   16.1          Y      Y
## 18         B         3                   10.3          Y      Y
## 19         B         4                    9.5          Y      Y
## 20         B         5                   15.3          Y      Y

Questions 1.2.1 - 1.2.4 are all concerned with whether one treatment is significantly different to another. This type of question can be approached by comparing the mean effects of treatments and running statistical tests which asks how likely it is that the observations being compared could have been drawn from the same population.

The statistical approach which comes to mind is that of hypothesis testing. The null hypothesis being that the samples being compared are not significantly different and the alternative hypothesis that the samples are indeed different. I could employ a t-test comparing sample means, or run an analysis of variance (ANOVA).

Since ANOVA runs pairwise comparisons of all the groups at once, and the fact that it is more resilient to deviations from the core assumptions, I believe that ANOVA is the best choice for these analyses.

My favourite way to run ANOVA in R is by using the lm() function. I will include the interaction effect in the model since it will relevant question 1.2.4:

lm_split <- lm(metabolised_felodipine ~ coffee * grapefruit, data = split_data)
summary(lm_split)

## 
## Call:
## lm(formula = metabolised_felodipine ~ coffee * grapefruit, data = split_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.560  -4.285  -0.830   3.225  18.940 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           66.200      3.362  19.691 1.22e-12 ***
## coffeeY              -26.120      4.755  -5.494 4.91e-05 ***
## grapefruitY          -44.140      4.755  -9.284 7.65e-08 ***
## coffeeY:grapefruitY   17.820      6.724   2.650   0.0175 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.518 on 16 degrees of freedom
## Multiple R-squared:  0.8994, Adjusted R-squared:  0.8805 
## F-statistic: 47.68 on 3 and 16 DF,  p-value: 3.343e-08

Anova(lm_split)

## Anova Table (Type II tests)
## 
## Response: metabolised_felodipine
##                   Sum Sq Df  F value    Pr(>F)    
## coffee            1480.9  1  26.2048  0.000103 ***
## grapefruit        6205.8  1 109.8108 1.427e-08 ***
## coffee:grapefruit  396.9  1   7.0238  0.017460 *  
## Residuals          904.2 16                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

1.2.1

Is there significant evidence that grapefruit juice reduced the extent to which felodipine is metabolised?

Yes, there is significant evidence that grapefruit reduced the extent to which felodipine is metabolised by CYP3A4.

The estimated effect of the grapefruit treatment is that felodipine metabolism is reduced by $44.140 ng/ml$ compared to the negative control. The associated p-value = $7.65e^{-08}$ providing strong evidence to reject the null hypothesis and accept the alternative that there a significant effect on felodipine metabolism.

1.2.2

Is there significant evidence that coffee reduces the extent to which felodipine is metabolised?

Yes, there is significant evidence that coffee reduces the extent to which felodipine is metabolised by CYP3A4.

The estimated effect of the coffee treatment is that felodipine metabolism is reduced by $26.120 ng/ml$ compared to the negative control. The associated p-value = $4.91e^{-05}$ providing strong evidence to reject the null hypothesis and accept the alternative that there a significant effect on felodipine metabolism.

1.2.3

Is there significant evidence that grapefruit juice reduces felodipine metabolism to a greater extent than does coffee?

Although the previous analysis indicated that grapefruit reduced felodipine metabolism by a greater extent ($26.120/44.140 = $) than coffee, there is no associated p-value and therefore no indication of whether the difference meets the significance threshold.

In order to determine if grapefruit juice reduces felodipine metabolism to a greater extent than coffee I will employ the aov() function so that I can run TukeyHSD(). The output from TukeyHSD() breaks down every combination of coffee and grapefruit, allowing me to directly compare the treatment effects.

aov_vitro <- aov(metabolised_felodipine ~ grapefruit * coffee, split_data)
summary(aov_vitro)

##                   Df Sum Sq Mean Sq F value   Pr(>F)    
## grapefruit         1   6206    6206 109.811 1.43e-08 ***
## coffee             1   1481    1481  26.205 0.000103 ***
## grapefruit:coffee  1    397     397   7.024 0.017460 *  
## Residuals         16    904      57                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Since the output from aov() is identical to that of Anova(lm_split), I’m happy to continue the analysis and run TukeyHSD().

vitro_tuk <- TukeyHSD(aov_vitro) %>% tidy()
vitro_tuk

## # A tibble: 8 × 7
##   term              contrast null.value estimate conf.low conf.high  adj.p.value
##   <chr>             <chr>         <dbl>    <dbl>    <dbl>     <dbl>        <dbl>
## 1 grapefruit        Y-N               0   -35.2    -42.4     -28.1  0.0000000143
## 2 coffee            Y-N               0   -17.2    -24.3     -10.1  0.000103    
## 3 grapefruit:coffee Y:N-N:N           0   -44.1    -57.7     -30.5  0.000000423 
## 4 grapefruit:coffee N:Y-N:N           0   -26.1    -39.7     -12.5  0.000259    
## 5 grapefruit:coffee Y:Y-N:N           0   -52.4    -66.0     -38.8  0.0000000385
## 6 grapefruit:coffee N:Y-Y:N           0    18.0      4.42     31.6  0.00785     
## 7 grapefruit:coffee Y:Y-Y:N           0    -8.30   -21.9       5.30 0.334       
## 8 grapefruit:coffee Y:Y-N:Y           0   -26.3    -39.9     -12.7  0.000239

The estimate associated with vitro_tuk$contrast == "N:Y-Y:N" indicates that the grapefruit only treatment reduced felodipine metabolism by $18.02ng/ml$ more than the coffee only treatment with an associated p-value = $7.85e^{-03}$. This indicates that there is significant evidence to suggest grapefruit juice reduces felodipine metabolism to a greater extent than coffee.

1.2.4

Are the effects of grapefruit juice and coffee additive; in other words, is the extent to which felodipine metabolism is reduced by their combined presence the same as the sum of their individual effects in isolation?

No, the effects of grapefruit juice and coffee are not additive.

The interaction term “coffeeY:grapefruitY” generated by lm_split tests whether the effect of the coffee or grapefruit treatment on metabolised felodipine depends on the application of the other treatment. A statistically significant p-value suggests that the combined treatment effect differs from what you would expect if treatment effects were additive.

The value of the interaction estimate indicates that felodipine metabolism by CYP3A4 is $17.82ng/ml$ less than you would expect if the treatments were additive. With an associated p-value = $3.34e^{-08}$, there is strong evidence to suggest that the combined effect is significantly different from the sum of their individual effects.

Part B: Clinical Trials

Task 1.3

For this study design and specific alternative hypothesis, what is the minimum group size N required to obtain a power of 0.8, assuming that a t-test will be used, the alternative is two sided, and the standard deviation of an individual’s systolic blood pressure is 20 mmHg.

power <- power.t.test(power = 0.8,
                      delta=10,
                      sd=20,
                      sig.level = 0.05,
                      type="two.sample",
                      alternative = "two.sided")
paste("For this study design and specific alternative hypothesis, the minimum group size required to obtain a power of 0.8, is", round(power$n, digits = 0), "participants per group.")

## [1] "For this study design and specific alternative hypothesis, the minimum group size required to obtain a power of 0.8, is 64 participants per group."

Task 1.4

The Design

500 subjects, all male with hypertension identified by “patient_index” (1-500)
- Numerical, integer data
Two treatment groups, one in which the medication was taken with 250ml of coffee, the other where the medication was taken without coffee.
- Factor (“coffee”) with two levels (“+” or “-”)
The measure is the difference in systolic blood pressure, given in units of mmHg, before taking felodipine and then 3 hours later.
- Paired observations, numerical double
If patients experienced an adverse reaction, no measurements were taken for blood pressure after treatment.
- Factor (“reaction”) with two levels (“0” or “1”)

The researchers hypothesised that consuming coffee with felodipine would lower systolic blood pressure by an average of 10mmHg more than when felodipine is taken without coffee.

Additional Thoughts

The trial isn’t blinded. It may improve the validity of the results if half of participants were randomly assigned decaf so that participants were not aware whether or not they were receiving the treatment.

trial_data <- read.csv("trial_felodipine.csv")
trial_data <- trial_data %>% 
  mutate(
    difference = systolic_before - systolic_after)
as_tibble(trial_data)

## # A tibble: 500 × 6
##    participant_index coffee systolic_before systolic_after reaction difference
##                <int> <chr>            <dbl>          <dbl>    <int>      <dbl>
##  1                 1 +                 127.           67.2        0       60  
##  2                 2 +                 146.          118.         0       28.4
##  3                 3 +                 137.          113.         0       24.7
##  4                 4 +                 159.          110          0       49.1
##  5                 5 +                 157.          128.         0       29  
##  6                 6 +                 152.          124.         0       27.4
##  7                 7 +                 155.          117.         0       37.5
##  8                 8 +                 160.          134.         0       25.8
##  9                 9 +                 150.           81.2        0       68.6
## 10                10 +                 167.          185.         0      -17.9
## # ℹ 490 more rows

trial_long <- trial_data %>%  
  pivot_longer(cols = c(systolic_before, systolic_after), 
               names_to = "timing", 
               values_to = "blood_pressure")  %>%
  group_by(timing) %>%
  mutate(mean = mean(blood_pressure, na.rm = T), sd = sd(blood_pressure, na.rm = T))
trial_long

## # A tibble: 1,000 × 8
## # Groups:   timing [2]
##    participant_index coffee reaction difference timing      blood_pressure  mean
##                <int> <chr>     <int>      <dbl> <chr>                <dbl> <dbl>
##  1                 1 +             0       60   systolic_b…          127.   161.
##  2                 1 +             0       60   systolic_a…           67.2  134.
##  3                 2 +             0       28.4 systolic_b…          146.   161.
##  4                 2 +             0       28.4 systolic_a…          118.   134.
##  5                 3 +             0       24.7 systolic_b…          137.   161.
##  6                 3 +             0       24.7 systolic_a…          113.   134.
##  7                 4 +             0       49.1 systolic_b…          159.   161.
##  8                 4 +             0       49.1 systolic_a…          110    134.
##  9                 5 +             0       29   systolic_b…          157.   161.
## 10                 5 +             0       29   systolic_a…          128.   134.
## # ℹ 990 more rows
## # ℹ 1 more variable: sd <dbl>

long_error <- trial_long  %>%
  group_by(timing) %>%
  summarise(mean_se(blood_pressure), sd = sd(blood_pressure, na.rm = T))
long_error

## # A tibble: 2 × 5
##   timing              y  ymin  ymax    sd
##   <chr>           <dbl> <dbl> <dbl> <dbl>
## 1 systolic_after   134.  133.  136.  23.7
## 2 systolic_before  161.  160.  161.  14.9

ggplot(trial_data, 
       aes(coffee, difference, colour = coffee)) +
  geom_jitter(
    width = 0.2,
    show.legend = FALSE,
    na.rm = T) +
  geom_boxplot(
    aes(fill = coffee, alpha = 0.2),
    width = 0.2,
    outlier.colour = NA,
    show.legend = FALSE,
    na.rm = T) +
  stat_summary(
    fun = "mean",
    geom = "point",
    shape = 23,
    fill = "white",
    na.rm = TRUE,
    show.legend = F) +
  coord_flip() +
    theme_bw() +
  labs(x = "Coffee Treatment",
       y = "Difference (mmHg)",
       title ="Box Plots and Jittered Points Representing the Difference is Systolic Blood \nPressure Before and After Felodipine Treatment",
       subtitle = "Mean shown as diamonds") +
  scale_y_continuous(breaks = seq(-30, 80, 10))

1. Is felodipine effective in lowering blood pressure?

Since blood pressure data are given by paired observations, I need to make sure I choose an appropriate test. I would like to run a paired t-test but first I need to check whether the data are normally distributed.

The first step is to visualise the distributions. I’ve used pivot_longer() to tidy the data so that it’s easier to produce these plots. I’m going to overlay a normal distribution generated using the mean and standard deviation parameters from the data so I can see how closely the data matches.

den_1 <- ggplot(trial_long %>%
                  filter(timing == "systolic_before"),
                aes(blood_pressure)) +
  stat_density(
    fill = "turquoise3",
    adjust = 1,
    kernel = "gaussian",
    show.legend = F,
    na.rm = T) +
  geom_function(fun = dnorm,
                args = list(mean = 160.6834, sd = 14.94013),
                na.rm = T) +
  theme_bw() +
  labs(x = "Systolic Blood Pressure (mmHg)",
       y = "Density",
       title ="Density Distribution of \nSystolic Blood Pressure (mmHg) \nBefore Treatment",
       subtitle = "Overlayed with a normal distribution")

den_2 <- ggplot(trial_long %>%
                  filter(timing == "systolic_after"),
                aes(blood_pressure)) +
  stat_density(
    fill = "indianred1",
    adjust = 1,
    kernel = "gaussian",
    show.legend = F,
    na.rm = T) +
  geom_function(
    fun = dnorm,
    args = list(mean = 134.3700, sd = 23.66987),
    na.rm = T) +
  theme_bw() +
  labs(x = "Systolic Blood Pressure (mmHg)",
       y = "Density",
       title ="Density Distribution of \nSystolic Blood Pressure (mmHg) \nAfter Treatment",
       subtitle = "Overlayed with a normal distribution")

den_1 + den_2 + plot_layout(ncol = 2)

The density plots seem to fit the assumption of normality really well. I’m also going to use shapiro.test() to generate more quantitative evidence before I decide whether I run a paired t-test.

shapiro.test(trial_data$systolic_before)

## 
##  Shapiro-Wilk normality test
## 
## data:  trial_data$systolic_before
## W = 0.99689, p-value = 0.4582

shapiro.test(trial_data$systolic_after)

## 
##  Shapiro-Wilk normality test
## 
## data:  trial_data$systolic_after
## W = 0.99503, p-value = 0.1872

The null hypothesis of a Shapiro-Wilk test is that the data is normally distributed. Since neither of our tests returned a significant result, there isn’t enough evidence to reject the null.

Now I’m feeling confident that the data meets the assumptions of a paired t-test, I can run it through the t.test() function.

t.test(Pair(systolic_before, systolic_after) ~ 1, data = trial_data)

## 
##  Paired t-test
## 
## data:  Pair(systolic_before, systolic_after)
## t = 25.965, df = 426, p-value < 2.2e-16
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  22.95966 26.72044
## sample estimates:
## mean difference 
##        24.84005

The output of our paired t-test has indicated that the blood pressure measurements taken after treatment with felodipine are significantly different.

Earlier, I used the mutate() function to calculate the difference between systolic_before and sytolic_after. A positive value reflects a decrease in blood pressure while a negative value indicates that blood pressure increased between measurements. By calculating the mean difference I can determine whether felodipine reduced or increased blood pressure.

mean(trial_data$difference, na.rm = T)

## [1] 24.84005

The mean difference is 24.84mmHg. Combining the results of all these tests provides evidence to suggest that felodipine effectively reduced systolic blood pressure. This is consistent with the findings of the in vitro study and clinical expectations.

ggplot(trial_long, 
       aes(timing, blood_pressure, colour = timing)) +
  geom_jitter(
    aes(alpha = 0.2),
    width = 0.2,
    show.legend = FALSE,
    na.rm = T) +
  geom_boxplot(
    aes(fill = timing, alpha = 0.2),
    width = 0.2,
    outlier.colour = NA,
    show.legend = FALSE,
    na.rm = T) +
  stat_summary(
    fun = "mean",
    geom = "point",
    shape = 23,
    fill = "white",
    na.rm = TRUE,
    show.legend = F) +
  coord_flip() +
    theme_bw() +
  labs(x = "Measured",
       y = "Systolic Blood Pressure (mmHg)",
       title ="Box Plots and Jittered Points Representing Systolic Blood \nPressure Before and After Felodipine Treatment",
       subtitle = "Means shown as diamonds") +
  scale_y_continuous(breaks = seq(0, 250, 20))

2. Does coffee consumption increase its effect?

To check whether coffee increases the effect of felodipine I will use lm() to generate a model comparing the observed difference in systolic blood pressure between subjects who were exposed to the coffee treatment and those who weren’t.

coffee_lm <- lm(difference ~ coffee, data = trial_data)
summary(coffee_lm)

## 
## Call:
## lm(formula = difference ~ coffee, data = trial_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -54.032 -12.875   0.268  13.697  49.025 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   21.232      1.310  16.202  < 2e-16 ***
## coffee+        7.443      1.882   3.955 8.98e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.44 on 425 degrees of freedom
##   (73 observations deleted due to missingness)
## Multiple R-squared:  0.03549,    Adjusted R-squared:  0.03322 
## F-statistic: 15.64 on 1 and 425 DF,  p-value: 8.983e-05

anova(coffee_lm)

## Analysis of Variance Table
## 
## Response: difference
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## coffee      1   5908  5908.4  15.638 8.983e-05 ***
## Residuals 425 160573   377.8                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The coefficients of the coffee+ provide an estimate of the difference between the treatment groups and the p-value associated with the estimate. According to the model, subjects who took felodipine with coffee experienced a mean reduction in systolic blood pressure 7.443mmHg greater than those who did not. The associated p-value = $8.98e^{-05}$. Therefore there is significant evidence to suggest that coffee consumption increased the effect of the felodipine treatment.

Task 1.5

The researchers were concerned that participants with more extreme hypertension had greater risk of an adverse reaction. Was there significant evidence of such a trend in the data and, if so, was the trend different between the two trial groups (coffee = “+” and coffee = “-”)? In response to these questions, be sure to describe your model/analysis and supply the corresponding R code. Provide quantitative support for your answers.

The question being asked is whether the continuous variable systolic_before can be used to predict a binary outcome “reaction” = “0” or “1”. A great way to approach this question is to generate an odds ratio for the likelihood of observing a reaction for every one unit increase in systolic_before. If the odds ratio is greater than one, then the likelihood of observing a reaction increases as systolic blood pressure increases.

trial_glm <- glm(reaction ~ systolic_before, family = binomial(), trial_data)
summary(trial_glm)

## 
## Call:
## glm(formula = reaction ~ systolic_before, family = binomial(), 
##     data = trial_data)
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -9.420625   1.521607  -6.191 5.97e-10 ***
## systolic_before  0.046617   0.009058   5.146 2.65e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 415.71  on 499  degrees of freedom
## Residual deviance: 387.01  on 498  degrees of freedom
## AIC: 391.01
## 
## Number of Fisher Scoring iterations: 5

trial_data$fitted.values <- trial_glm$fitted.values

as_tibble(trial_data)

## # A tibble: 500 × 7
##    participant_index coffee systolic_before systolic_after reaction difference
##                <int> <chr>            <dbl>          <dbl>    <int>      <dbl>
##  1                 1 +                 127.           67.2        0       60  
##  2                 2 +                 146.          118.         0       28.4
##  3                 3 +                 137.          113.         0       24.7
##  4                 4 +                 159.          110          0       49.1
##  5                 5 +                 157.          128.         0       29  
##  6                 6 +                 152.          124.         0       27.4
##  7                 7 +                 155.          117.         0       37.5
##  8                 8 +                 160.          134.         0       25.8
##  9                 9 +                 150.           81.2        0       68.6
## 10                10 +                 167.          185.         0      -17.9
## # ℹ 490 more rows
## # ℹ 1 more variable: fitted.values <dbl>

exp(coef(trial_glm))

##     (Intercept) systolic_before 
##    8.103532e-05    1.047720e+00

The output from summary(trial_glm) provides the coefficients we need to calculate an odds ratio. The calculated odds ratio was 1.047720, meaning that for every one unit increase in systolic_before, the odds of observing a reaction increased by a factor of $1.05$, or $5\%$.

ggplot(trial_data,
       aes(x = systolic_before, y = reaction)) +
  scale_x_continuous(
    breaks = seq(120, 210, 10),
    minor_breaks = NULL) +
  scale_y_continuous(
    breaks = c(0, 1),
    labels = c("No", "Yes"),
    minor_breaks = NULL) +
  geom_jitter(
    aes(colour = factor(reaction),
        fill = factor(reaction)),
    width = 0,
    height = 0.07,
    alpha = 0.5,
    show.legend = F) +
  geom_line(
    aes(y = fitted.values)) +
  theme_bw() +
  labs(x = "Systolic Blood Pressure Before Treatment (mmHg)",
       y = "Reaction",
       title ="Reaction Occurrence Plotted Against Systolic Blood Pressure \nBefore Treatment")

The second part of the task asks if the trend differs depending on whether or not the subject consumed coffee?

In statistical terms this is basically asking whether the function used to generate our model would be a better fit to the data if we included the coffee variable as a predictor. Employing a Likelihood Ratio Test (LRT) allows us to compare nested models, it’s output indicates whether there is significant evidence to suggest that a complex model is a better fit to the data than the simple model in which it is nested.

To perform my LRT I will use trial_glm as the null model and create a new model called coffee_glm to act as the alternative model. Then I can use the function anova(null, alt, test = "Chisq") to generate an analysis of deviance table. If anova() returns a significant p-value, this would indicate that the alternative model is a better fit to the data.

ggplot(trial_data,
       aes(factor(coffee), factor(reaction))) +
  geom_jitter(
    aes(colour = systolic_before),
    width = 0.2,
    height = 0.2,
    show.legend = F) +
  coord_fixed() +
  theme_bw() +
  labs(x = "Reaction",
       y = "Coffee",
       title ="Reaction Occurence Separated by \nCoffee Treatment")

ggplot(trial_data,
       aes(systolic_before, colour = factor(reaction),
           fill = factor(reaction),
           alpha = 0.5)) +
  geom_density(
    show.legend = F,
    trim = F) +
  theme_bw() +
  labs(x = "Systolic Blood Pressure Before Treatment (mmHg)",
       y = "Density",
       title ="Distribution of Systolic Blood Pressure Before Treatment \nSeparated by Reaction Occurence")

coffee_glm <- glm(reaction ~ systolic_before * coffee, family = binomial(), trial_data)
summary(coffee_glm)

## 
## Call:
## glm(formula = reaction ~ systolic_before * coffee, family = binomial(), 
##     data = trial_data)
## 
## Coefficients:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             -10.272782   2.355836  -4.361  1.3e-05 ***
## systolic_before           0.050494   0.013999   3.607  0.00031 ***
## coffee+                   1.623145   3.092238   0.525  0.59965    
## systolic_before:coffee+  -0.007455   0.018398  -0.405  0.68535    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 415.71  on 499  degrees of freedom
## Residual deviance: 384.81  on 496  degrees of freedom
## AIC: 392.81
## 
## Number of Fisher Scoring iterations: 5

anova(coffee_glm)

## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: reaction
## 
## Terms added sequentially (first to last)
## 
## 
##                        Df Deviance Resid. Df Resid. Dev Pr(>Chi)    
## NULL                                     499     415.71             
## systolic_before         1  28.7001       498     387.01 8.45e-08 ***
## coffee                  1   2.0313       497     384.98   0.1541    
## systolic_before:coffee  1   0.1646       496     384.81   0.6849    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(trial_glm, coffee_glm, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: reaction ~ systolic_before
## Model 2: reaction ~ systolic_before * coffee
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       498     387.01                     
## 2       496     384.81  2   2.1959   0.3336

With P(>Chi) = 0.3336, the results of the analysis do not provide evidence to suggest that the alternative model provides a significantly better fit to the data. This indicates that the trend does not differ depending on coffee consumption.

Problem Two

Please use the following structure to respond to each question:

Based on the question, indicate what you are looking for in the plot.
Describe what you see, and how it compares to what you are looking for. Finally, draw your conclusion in direct response to the question.

Task 2.1

Does the data show evidence of a tray effect? Carefully explain your reasoning as instructed above.

If there is no tray effect, I would expect that trays in the same row (having received the same irrigation treatment) will have equivalent distributions. Therefore if the distributions are visually different between trays of the same row, this might suggest evidence of a tray effect.

Trays 1, 6, and 8

Trays 6 and 8 appear to have equivalent distributions with maximum growth peaking at ~25. For tray 1 the distribution appears to be shifted up by ~15 units, hinting at the presence of a tray effect

Trays 2, 4, and 9

The distributions seem to shift by ~5-10 units between the trays. Tray 9 maxes out at ~15, tray 2 at ~20, and tray 4 and ~30. Again suggesting that a tray effect is present.

Trays 3, 5, and 7

The distributions for these trays are very similar. The only obvious difference being that the distribution for tray 7 is ~5 units lower than trays 3 and 5. Otherwise I’m not confident to say that this row indicates the presence of a tray effect.

Overall

The variance between trays of the same irrigation level seems to suggest that a tray effect is present.

Task 2.2

Does the data show evidence of a variety effect? Carefully explain your reasoning as instructed above.

I expect that if there is no variety effect, there would be limited separation of data points by colour. That is to say, both colours would be more or less randomly dispersed along the height axis. Because there are only four data points per treatment group, there may appear to be a variety effect for some of the treatment groups even if a variety effect doesn’t appear in the model.

In other words, knowing that pattern recognition has played such an important role in human evolution, to the extent that we tend to see patterns where none exist (search pareidolia), I think that the risk of a false positive is very high. For this reason I’m approaching this question very much biased towards the assumption that no variety effect exists.

My answer

I don’t think that there is enough separation between colours to suggest evidence of a variety effect.

Task 2.3

Does the data show evidence of a soil effect? Carefully explain your reasoning as instructed above.

If no soil effect is present, mean heights between soil level treatments would be the same. I would expect therefore, that height distributions would appear equivalent for the different levels of soil within a tray.

In the hopes of effectively conveying my thought process, I’m going to structure my analysis as follows.

If there is no soil effect, it would look like this:

Tray $x$

$S \approx R \approx N \approx B$

If the actual trays don’t broadly resemble the above example, that would suggest the presence of a soil effect.

Tray 1

$S < R \le N > B$

Tray 2

$S < R < N > B$

Tray 3

$S < R \approx N < B$

Tray 4

$S < R = N > B$

Tray 5

$S < R < N < B$

Tray 6

$S < R \approx N > B$

Tray 7

$S < R \approx N < B$

Tray 8

$S \le R \le N > B$

Tray 9

$S \le R \approx N \ge B$

Given the variation present within trays, there does appear to be evidence of a soil effect.

Task 2.4

Do you see evidence of any pairwise interaction between two of water, variety and soil? In other words, what does the plot suggest about the pairwise interactions (1) water:variety, (2) water:soil and (3) variety:soil? For each of these, carefully explain your reasoning as instructed above.

1. water:variety

If no interaction effect is present I would expect that as you move down a column between irrigation levels the order in which varieties appear along the height axis would be more or less the same.

For example, if the drought resistant variety had the greatest height in tray 1, then I would expect that to be true for trays 2 and 3 regardless of the effect of irrigation. However, as discussed when looking into variety effect, the low per treatment sample size makes visual analysis difficult unless the effect is very large.

Given the low spread of the data, I don’t have confidence that any deviations from what I expect from a lack of interaction can’t be attributed to standard deviation. For this reason I do not think there is enough evidence to suggest an interaction effect between water and variety.

2. water:soil

If an the is no water:soil interaction effect, I would expect that as you move down columns between irrigation levels, the difference in mean heights between soil treatments would be maintained.

This appears to be the case as you move from low to medium irrigation. Mean height seems to increase as you move from L/M:S -> L/M:R -> L/M:N before decreasing as you move from L/M:N -> L/M:B. When you look at the trays exposed to the high irrigation treatment, mean height increases from H:S -> H:R -> H:N, similar to the low and medium irrigation conditions. However, in contrast to those treatments, mean height increases for H:B.

This suggest to be that there is minimal water:soil interactions for low and medium irrigation treatments but a substantial interaction effect for plants exposed to the high irrigation treatment. This effect is apparent for the interaction related to water = “H” and soil = “B”.

3. variety:soil

To determine whether an interaction effect exists between variety and soil, I am looking at the change in height between soil treatments within trays, then comparing between trays to see if my observations are consistent. If there is no observable effect, I would expect that the differences in height will be roughly equivalent regardless of whether you include both varieties in the observation.

In other words, if the mean height difference between treatment “S” and treatment “R” was 10 units when you included observations from both varieties, then I would expect that to be true for when calculating the same difference for each variety independently.

Only tray 8 seems to show a difference in the effect of soil treatment between the two varieties and only for “R” and “N”. Since this effect doesn’t appear to be true for any other trays, I believe that there is insufficient evidence to suggest that a variety:soil interaction exists.

Problem Three

The Design

These data originate from an observational study, not an experiment.

200 subjects aged 90yrs old, each given a patient ID (1-200)
- Numeric, integer
Outcome variable - Alzheimer’s diagnosis (Yes or No)
- Factor with two levels
For subjects with Alzheimer’s, age of onset is given, (Years of age)
- Numeric, integer
Fat consumption in units of grams per week (g/wk)
- Numeric, integer
Physical activity measured in units of hours per week (hr/wk)
- Numeric, double
Fat consumption and physical activity data relies on self reporting 50-70 years after the fact.

The research questions:

Do those two lifestyle choices have significant effect on the risk of getting the disease?
Does the data support the hypothesis that people who exercise more and eat less saturated fats develop Alzheimer’s disease later in life?

Additional Thoughts

I’m not confident in the reliability of the data. Self reported data, especially in relation to diet and exercise, is highly prone to biased or outright false responses. Additionally, the researchers are asking subjects who are 90 years old, some of whom have Alzheimer’s, to recall how they ate and exercised between 50-70 years in the past. If you can accurately recall how you were eating and exercising this time last month I would be shocked (unless you’re a body builder like Marcin and Eric appear to be). Furthermore, saturated fat consumption is a very specific metric that would need to be estimated based on the types of food being eaten. If someone presented me with the data in real life, I would tell them to come back when they had something real.

Task 3.1

alz_data <- read.csv("alzheimer.csv", stringsAsFactors = TRUE)
alz_data <- alz_data %>%
  mutate(exercise_adj = exercise_hours + 1)
as_tibble(alz_data)

## # A tibble: 200 × 6
##    patient_id disease age_onset fat_consumption exercise_hours exercise_adj
##         <int>   <int>     <int>           <int>          <dbl>        <dbl>
##  1          1       1        55              90            1.8          2.8
##  2          2       0        NA              80            5.6          6.6
##  3          3       0        NA              83            6.7          7.7
##  4          4       0        NA              68            0.4          1.4
##  5          5       1        72              87            2.6          3.6
##  6          6       1        78              86            4.5          5.5
##  7          7       1        61              74            1.9          2.9
##  8          8       0        NA              53            2.9          3.9
##  9          9       1        79              94            0.7          1.7
## 10         10       1        74              91            2.6          3.6
## # ℹ 190 more rows

alz_stat <- alz_data %>%
  group_by(disease) %>%
  summarise(
    fat_mean = mean(fat_consumption),
    fat_sd = sd(fat_consumption),
    ex_mean = mean(exercise_hours),
    ex_sd = sd(exercise_hours))
as_tibble(alz_stat)

## # A tibble: 2 × 5
##   disease fat_mean fat_sd ex_mean ex_sd
##     <int>    <dbl>  <dbl>   <dbl> <dbl>
## 1       0     71.1   15.8    2.73  1.90
## 2       1     88.8   17.3    1.68  1.83

Let’s begin by visualising the distributions of the parameters. I want to see the following:

The frequency distribution of fat_consumption
The frequency distribution of exercise_hours
The frequency distribution of age_onset
The differences in the distributions of fat_consumption and exercise_hours between those with without the disease

Density Distributions for fat_consumption

den_3 <- ggplot(alz_data %>%
                     filter(disease == "0"),
                   aes(fat_consumption)) +
  stat_density(
    fill = "indianred1",
    adjust = 1,
    kernel = "gaussian",
    show.legend = F) +
  geom_function(fun = dnorm, args = list(mean = 71.05607, sd = 15.77904)) +
   theme_bw() +
  labs(x = "Saturated Fat Consumption (g/wk)",
       y = "Density",
       title ="Distribution of Saturated Fat Consumption \nfor Subjects without Alzheimer's",
       subtitle = "Normal Distribution Overlayed")

den_4 <- ggplot(alz_data %>%
                  filter(disease == "1"),
                aes(fat_consumption)) +
  stat_density(
    fill = "turquoise3",
    adjust = 1,
    kernel = "gaussian",
    show.legend = F) +
  geom_function(fun = dnorm, args = list(mean = 88.76344, sd = 17.28684)) +
   theme_bw() +
  labs(x = "Saturated Fat Consumption (g/wk)",
       y = "Density",
       title ="Distribution of Saturated Fat Consumption \nfor Subjects with Alzheimer's",
       subtitle = "Normal Distribution Overlayed")

den_5 <- ggplot(alz_data,
                aes(fat_consumption)) +
  stat_density(
    fill = "goldenrod",
    adjust = 1,
    kernel = "gaussian",
    show.legend = F) +
  geom_function(fun = dnorm, args = list(mean = mean(alz_data$fat_consumption), sd = sd(alz_data$fat_consumption))) +
   theme_bw() +
  labs(x = "Saturated Fat Consumption (g/wk)",
       y = "Density",
       title ="Overall Distribution of Saturated \nFat Consumption",
       subtitle = "Normal Distribution Overlayed") 

patchwork::wrap_plots(den_3, den_4, den_5, ncol = 2)

Density Distributions for `exercise_hours’

den_6 <- ggplot(alz_data %>%
                     filter(disease == "0"),
                   aes(exercise_hours)) +
  stat_density(
    fill = "indianred1",
    adjust = 1,
    kernel = "gaussian",
    show.legend = F) +
  geom_function(fun = dnorm, args = list(mean = 2.727103, sd = 1.826946)) +
   theme_bw() +
  labs(x = "Weekly Exercise Hours",
       y = "Density",
       title ="Distribution of Exercise Hours \nfor Subjects without Alzheimer's",
       subtitle = "Normal Distribution Overlayed") 

den_7 <- ggplot(alz_data %>%
                  filter(disease == "1"),
                aes(exercise_hours)) +
  stat_density(
    fill = "turquoise3",
    adjust = 1,
    kernel = "gaussian",
    show.legend = F) +
  geom_function(fun = dnorm, args = list(mean = 1.679570, sd = 1.826946)) +
  theme_bw() +
  labs(x = "Weekly Exercise Hours",
       y = "Density",
       title ="Distribution of Exercise Hours \nfor Subjects with Alzheimer's",
       subtitle = "Normal Distribution Overlayed") 

den_8 <- ggplot(alz_data,
                aes(exercise_hours)) +
  stat_density(
    fill = "goldenrod",
    adjust = 1,
    kernel = "gaussian",
    show.legend = F) +
  geom_function(fun = dnorm, args = list(mean = mean(alz_data$exercise_hours), sd = sd(alz_data$exercise_hours))) +
  theme_bw() +
  labs(x = "Weekly Exercise Hours",
       y = "Density",
       title ="Overal Distribution of Exercise Hours",
       subtitle = "Normal Distribution Overlayed") 

wrap_plots(den_6, den_7, den_8, nrow = 2)

A qualitative assessment of the data clearly shows that fat_consumption and age_onset are not normally distributed. While fat_consumption appears closer to a normal distribution, I’d like to run a Shapiro-Wilk normality test so that I have a quantitative assessment.

shapiro.test(alz_data$fat_consumption)

## 
##  Shapiro-Wilk normality test
## 
## data:  alz_data$fat_consumption
## W = 0.99517, p-value = 0.774

shapiro.test(alz_data$exercise_hours)

## 
##  Shapiro-Wilk normality test
## 
## data:  alz_data$exercise_hours
## W = 0.92196, p-value = 8.125e-09

The output above aligns with my qualitative assessment that fat_consumption data are normally distributed but exercise_hours are not. This might result in some issues when running statistical analysis, so transforming the exercise_hours data could be worthwhile.

Since exercise_hours is skewed heavily to the left, performing a log transformation might result in a distribution which more closely resembles normality.

alz_log <- log(alz_data$exercise_hours + 1)
alz_data$exercise_log <- alz_log

shapiro.test(alz_data$exercise_log)

## 
##  Shapiro-Wilk normality test
## 
## data:  alz_data$exercise_log
## W = 0.93132, p-value = 4.306e-08

Comparing the p-values of normality test on the transformed and original data suggests that although neither are normally distributed, there is a slight improvement in normality after the transformation. However, I don’t think that this minuscule improvement is enough to justify employing the transformed data in favour of the original.

ggplot(alz_data,
       aes(fat_consumption)) +
  geom_histogram(binwidth = 1) +
  theme_bw() +
  labs(x = "Saturated Fat Consumption (g/wk)",
       y = "Count",
       title ="Histogram Showing Saturated Fat Consumption (g/wk)")

ggplot(alz_data,
       aes(exercise_hours)) +
  geom_histogram(binwidth = 0.1) +
  theme_bw() +
  labs(x = "Physical Activity (hr/wk)",
       y = "Count",
       title ="Histogram Showing Hours of Physical Activity Per Week")

ggplot(alz_data,
       aes(age_onset)) +
  geom_histogram(
    binwidth = 1,
    na.rm = TRUE) +
  theme_bw() +
  labs(x = "Age Onset",
       y = "Count",
       title ="Histogram Showing Age of Alzheimer's Onset")

fat_density <- ggplot(alz_data,
                      aes(fat_consumption, fill = factor(disease), colour = factor(disease), alpha = 0.5)) +
  geom_density(show.legend = F) +
  theme_bw() +
  labs(x = "Saturated Fat Consumption (g/wk)",
       y = "Density",
       title ="Density Distribution For Saturated Fat Consumption (g/wk)")

phys_density <- ggplot(alz_data,
                       aes(exercise_hours, fill = factor(disease), colour = factor(disease), alpha = 0.5)) +
  geom_density(show.legend = F) +
  theme_bw() +
  labs(x = "Physical Activity (hr/wk)",
       y = "Density",
       title ="Density Distribution For Hours of Physical Activity Per Week")

onset_density <- ggplot(alz_data,
                       aes(age_onset)) +
  geom_density(
    show.legend = F,
    fill = "goldenrod",
    colour = "goldenrod",
    na.rm = T) +
  theme_bw() +
  labs(x = "Age Onset",
       y = "Density",
       title ="Density Distribution For Age of Alzheimer's Onset")

fat_density

phys_density

onset_density

fat_glm <- glm(disease ~ fat_consumption, family = binomial(), alz_data)
fat_glm

## 
## Call:  glm(formula = disease ~ fat_consumption, family = binomial(), 
##     data = alz_data)
## 
## Coefficients:
##     (Intercept)  fat_consumption  
##         -5.3002           0.0647  
## 
## Degrees of Freedom: 199 Total (i.e. Null);  198 Residual
## Null Deviance:       276.3 
## Residual Deviance: 226   AIC: 230

alz_data$fat.fit <- fat_glm$fitted.values

ggplot(alz_data,
       aes(fat_consumption, disease)) +
  scale_x_continuous(breaks = seq(0, 130, 10),
                     minor_breaks = NULL) +
  scale_y_continuous(breaks = c(0, 1),
                     labels = c("No", "Yes"), minor_breaks = NULL) +
  geom_jitter(
    aes(colour = factor(disease), fill = factor(disease)),
    width = 0,
    height = 0.1, 
    alpha = 0.5,
    show.legend = F) +
  geom_line(aes(y = fat.fit)) +
  theme_bw() +
  labs(x = "Saturated Fat Consumption (g/wk)",
       y = "Alzheimer's",
       title ="Alzheimer's Diagnosis Plotted Against Saturated Fat Consumption")

phys_glm <- glm(disease ~ exercise_hours, family = binomial(), alz_data)
phys_glm

## 
## Call:  glm(formula = disease ~ exercise_hours, family = binomial(), 
##     data = alz_data)
## 
## Coefficients:
##    (Intercept)  exercise_hours  
##         0.5248         -0.3068  
## 
## Degrees of Freedom: 199 Total (i.e. Null);  198 Residual
## Null Deviance:       276.3 
## Residual Deviance: 260.9     AIC: 264.9

alz_data$phys.fit <- phys_glm$fitted.values

ggplot(alz_data,
       aes(exercise_hours, disease)) +
  scale_x_continuous(breaks = seq(0, 9, 1),
                     minor_breaks = NULL) +
  scale_y_continuous(breaks = c(0, 1),
                     labels = c("No", "Yes"),
                     minor_breaks = NULL) +
  geom_jitter(
    aes(colour = factor(disease),
        fill = factor(disease)),
    width = 0,
    height = 0.1, 
    alpha = 0.5,
    show.legend = F) +
  geom_line(
    aes(y = phys.fit)) +
    theme_bw() +
  labs(x = "Exercise (hr/wk)",
       y = "Alzheimer's",
       title ="Alzheimer's Diagnosis Plotted Against Weekly Exercise")

Task 3.2

Using the Likelihood Ratio Test (LRT), answer the main research question: Do those two lifestyle choices have significant effect on the risk of getting the disease?

Performing an LRT requires the generation of a null model and an alternative model. In this case the null model is quite simple to produce.

The null model (H₀)

This model will not account for exercise hours, nor saturated fat eaten. The model implies that Alzheimer’s diagnoses have no predictors.

mod_null <- glm(disease ~ 1, data = alz_data, family = binomial)
summary(mod_null)

## 
## Call:
## glm(formula = disease ~ 1, family = binomial, data = alz_data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -0.1402     0.1418  -0.989    0.323
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 276.28  on 199  degrees of freedom
## Residual deviance: 276.28  on 199  degrees of freedom
## AIC: 278.28
## 
## Number of Fisher Scoring iterations: 3

The alternative model (H₁)

There are many different options for how I construct the alternative model. Thankfully I don’t need to choose. Instead, I can make every combination of the alternative model and put them all in the same anova() along with the null model.

mod_alt <- glm(disease ~ fat_consumption + exercise_hours, data = alz_data, family = binomial)
summary(mod_alt)

## 
## Call:
## glm(formula = disease ~ fat_consumption + exercise_hours, family = binomial, 
##     data = alz_data)
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -4.51496    0.96435  -4.682 2.84e-06 ***
## fat_consumption  0.05919    0.01105   5.356 8.52e-08 ***
## exercise_hours  -0.15687    0.09265  -1.693   0.0904 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 276.28  on 199  degrees of freedom
## Residual deviance: 223.06  on 197  degrees of freedom
## AIC: 229.06
## 
## Number of Fisher Scoring iterations: 4

mod_alt_int <- glm(disease ~ fat_consumption * exercise_hours, data = alz_data, family = binomial)
summary(mod_alt_int)

## 
## Call:
## glm(formula = disease ~ fat_consumption * exercise_hours, family = binomial, 
##     data = alz_data)
## 
## Coefficients:
##                                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    -4.734091   1.357067  -3.488 0.000486 ***
## fat_consumption                 0.061893   0.016151   3.832 0.000127 ***
## exercise_hours                 -0.053718   0.449859  -0.119 0.904949    
## fat_consumption:exercise_hours -0.001321   0.005640  -0.234 0.814833    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 276.28  on 199  degrees of freedom
## Residual deviance: 223.01  on 196  degrees of freedom
## AIC: 231.01
## 
## Number of Fisher Scoring iterations: 4

mod_alt_fat <- glm(disease ~ fat_consumption, data = alz_data, family = binomial)
summary(mod_alt_fat)

## 
## Call:
## glm(formula = disease ~ fat_consumption, family = binomial, data = alz_data)
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -5.3002     0.8785  -6.033 1.61e-09 ***
## fat_consumption   0.0647     0.0108   5.992 2.07e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 276.28  on 199  degrees of freedom
## Residual deviance: 225.98  on 198  degrees of freedom
## AIC: 229.98
## 
## Number of Fisher Scoring iterations: 4

mod_alt_ex <- glm(disease ~ exercise_hours, data = alz_data, family = binomial)
summary(mod_alt_ex)

## 
## Call:
## glm(formula = disease ~ exercise_hours, family = binomial, data = alz_data)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     0.52480    0.22635   2.319 0.020421 *  
## exercise_hours -0.30678    0.08287  -3.702 0.000214 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 276.28  on 199  degrees of freedom
## Residual deviance: 260.94  on 198  degrees of freedom
## AIC: 264.94
## 
## Number of Fisher Scoring iterations: 4

anova(mod_null, mod_alt_fat, mod_alt_ex, mod_alt, mod_alt_int, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: disease ~ 1
## Model 2: disease ~ fat_consumption
## Model 3: disease ~ exercise_hours
## Model 4: disease ~ fat_consumption + exercise_hours
## Model 5: disease ~ fat_consumption * exercise_hours
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1       199     276.28                          
## 2       198     225.98  1   50.296 1.322e-12 ***
## 3       198     260.94  0  -34.955              
## 4       197     223.06  1   37.875 7.541e-10 ***
## 5       196     223.01  1    0.054    0.8158    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The results of the analysis of deviance indicate that there are two models which are a better fit to the data that the null, mod_alt_fat and mod_alt. mod_alt_fat returned a deviance of 50.296 with a p-value = $1.322e^{-12}$ and mod_alt returned a deviance of $37.875$ with a p-value = $7.541e^{-10}$.

With a lower deviance and AIC, the analysis suggest that mod_alt provides the best fit to the data.

summary(mod_alt)

## 
## Call:
## glm(formula = disease ~ fat_consumption + exercise_hours, family = binomial, 
##     data = alz_data)
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -4.51496    0.96435  -4.682 2.84e-06 ***
## fat_consumption  0.05919    0.01105   5.356 8.52e-08 ***
## exercise_hours  -0.15687    0.09265  -1.693   0.0904 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 276.28  on 199  degrees of freedom
## Residual deviance: 223.06  on 197  degrees of freedom
## AIC: 229.06
## 
## Number of Fisher Scoring iterations: 4

So, do the two lifestyle choices have significant effect on the risk of getting the disease?

The likelihood ratio tests has provided us with strong evidence to suggest that fat consumption and exercise hours do have significant effect on the risk of developing Alzheimer’s. When we examine the coefficients of the model above, we can spot something very interesting. There is strong evidence (p-value = $8.52e^{-08}$) to suggest that increased fat consumption leads to increased Alzheimer’s risk, however a significant effect is not observed for exercise hours… We may need to investigate this further.

The final thing I would like to investigate is the how the model predicts the likelihood of observing a positive diagnosis for every one unit increase of the predictors. In other words, I want to generate the odds ratios.

exp(coef(mod_alt))

##     (Intercept) fat_consumption  exercise_hours 
##      0.01094407      1.06097945      0.85481627

According to the calculated odds ratios, for every one unit increase in fat_consumption, the likelihood of observing a positive diagnosis increases by $6\%$. In contrast, for every one unit increase in exercise hours the likelihood of diagnosis decreases by $15\%$.

I think that these results are really interesting because even though the associated by value is insufficient to suggest a exercise hours confer a significant effect, the odds ratio generated indicates that just an extra 6 minutes of exercise a week reduces the likelihood of Alzheimer’s diagnosis by $15\%$! With Alzheimer’s and similar degenerative neurological disorders quickly becoming the leading cause of death, I’m tempted to say that even weak evidence for an effect of this magnitude warrants further investigation.

Task 3.3

Using the LRT answer: Do we observe a statistically significant effect of the amount of exercise on the risk of the disease? Perform statistical test and write your conclusion in 2 - 5 sentences.

No, we do not observe a statistically significant effect of the amount of exercise on the risk of Alzheimer’s.

In regards to the LRT performed in Task 3.2, I believe that my previous analysis has sufficiently addressed the question. Additionally, using the output of mod_alt_ex to perform a new analysis would not be wise given that we previously determined it is not a better fit to the data than the null model.

Task 3.4

Do we observe a statistically significant effect of the amount of exercise on the age of onset of the disease? Choose and perform statistical test or tests, and write your conclusion in 2 - 5 sentences.

I had a really difficult time deciding on an approach to this question. Neither distribution is normally distributed and the outcome is not binary so I can’t apply an of the statistical approaches from previous tasks.

I ended up deciding to employ a correlation test and after some online investigation decided that Kendall’s Tau would be most suited to these data.

ggplot(alz_data,
       aes(age_onset, mean(exercise_hours))) +
  geom_col(na.rm = T,
           show.legend = F) +
  theme_bw() +
  labs(x = "Age Onset",
       y = "Mean Exercise (hr/wk)",
       title ="Histogram Displaying Mean Exercise by Age of Alzheimer's Onset")

cor.test(alz_data$exercise_hours, alz_data$age_onset, method = "kendall")

## 
##  Kendall's rank correlation tau
## 
## data:  alz_data$exercise_hours and alz_data$age_onset
## z = -1.2067, p-value = 0.2276
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
##         tau 
## -0.08895736

Kendall’s rank correlation tau did not produce evidence to suggest that a significant correlation exist between exercise hours and age of onset. The output generated provided an insignificant p-value = $0.227$ for an extremely weak negative correlation ($\tau = -0.08895736$).

Task 3.5

Does the data support the hypothesis that people who exercise more tend to eat less saturated fats? You have freedom of choice of the statistical method. Write your conclusion in 2 - 5 sentences.

Similar to the previous task, I will use cor.test(..., method = "kendall") to generate a correlation coefficient.

ggplot(alz_data, aes(exercise_hours, fat_consumption, colour = disease)) +
  geom_jitter(width = 0.1, show.legend = FALSE) +
  geom_smooth(method = "lm", aes(colour = NULL), show.legend = FALSE)

## `geom_smooth()` using formula = 'y ~ x'

cor.test(alz_data$exercise_hours, alz_data$fat_consumption, method = "kendall")

## 
##  Kendall's rank correlation tau
## 
## data:  alz_data$exercise_hours and alz_data$fat_consumption
## z = -5.2524, p-value = 1.502e-07
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
##        tau 
## -0.2562325

With a p-value = $1.502e^{-07}$, there is evidence to suggest a significant negative correlation between exercise hours and saturated fat consumption ($\tau = -0.2562325$). Therefore I believe that there is significant evidence to suggest that those who exercise more tend to eat less saturated fats.

BIOL2202: Assignment 02

Matthew Nakkan

2024-10-12