Introduction

In this project, I aimed to answer various questions through data analysis in R. The main goal of the analysis was to discover the linear model used to produce the synthetic data used in this project. This data is based on Influence of Life Stress on Depression a paper by Caspi et al. that tells a finding of gene-environment interaction. It focuses on a genetic factor with two variants relating to serotonin. Someone who has the shorter allele variant was found to have a stronger relationship between stressful events and depression than someone with the long allele variant. While the data itself is synthetic, it provides a situation of analysis quite similar to the one in that paper. To answer the question of ‘what model best fits my data?’ I needed to ask simpler questions first.

Methods

After importing my data into Rstudio, I tested to make sure there was no missing data. After I ensured everything was there, I began to plan how I was going to find the linear regression model. I started by plotting scatter plots of Y with each of the four environmental variables. Using the ggplot2 package I graphed each of these and found that E1 and E4 had a very weak relationship while E2 and E3 showed a positive linear relationship.

I then followed this by creating separate linear models with each environmental variable and found again that E2 and E3 appeared to have a stronger relationship with Y. The summary function on each of these smaller models showed that E2 had a t-value of 14.22 with a very low p-value and E3 had a t-value of 20.35 with a very low p-value showing that both of these variables are significant. E1 and E4 both had t-values near or less than 1 and p-values higher than 0.1 making them appear not significant.

Generalized Additive Model

From here I decided to use the Generalized Additive Model to see if my variables might be non-linear. A summary of the model showed again that E1 and E2 were not significant but also revealed that E2 and especially E3 may have a nonlinear relationship. This was shown with their edf values being above one and E3 being above two.

To confirm I had polynomial terms I then made models of each variable and raised them by a degree to see if they would be better fits. Unfortunately, my results for this showed only worsening t-values and adjusted R2 values for each of the models, so I began to question if they were nonlinear. To continue my search for significant variables I decided to calculate the strongest correlations to Y in my variables. I used dplyr to find and then sort these correlations from strongest to weakest. Doing this showed that my most significant variables were E3, E2, and G12 with correlations all above 0.1 with G14 being the next highest near 0.07.

##         E3         E2        G12        G14 
## 0.48576742 0.36202713 0.14581856 0.06818992

Early Linear Models

I finally created my first linear model using all of the environmental variables. This model had an adjusted R2 of 0.3769 and upon viewing the residual plot there was a slight pattern. In addition, the Normal Q-Q plot showed that it was not normally distributed.

From here I introduced the rest of my variables and raised them to the second order to include any interaction variables. This model had an adjusted R2 of 0.402, this was a slight improvement over the previous model. Upon plotting the model with all variables, this residual plot showed even more of a pattern but the Q-Q plot was closer to seeming normal. This was a tradeoff and a slight improvement over the environmental-only model.

As it seemed that I had acquired a substantial amount of information on my independent variables, I looked at my dependent variable. I decided to use the ‘box cox’ method from the MASS package. The method gave me a lambda of 1.9 which I then used to transform my Y variable to Y1.9.

This new model (Y1.9 = (E1+E2+E3+E4+G1+…G20)2 gave me an adjusted R2 of 0.4014. While this was slightly worse, my residual plot improved along with the Q-Q plot, showing that the model was more normal and that there was less heteroscedasticity.

Lasso and Stepwise Regression

I decided to use the ‘LASSO’ method to narrow down my significant variables from here. I reconfirmed that E3 was significant and was able to learn that, according to this method, G12:G14 had the highest impact on Y as its coefficient was the highest. E2 was not shown on its own, only listed with E3 as an interaction variable.

##     variable coefficient
## 266  G12:G14   4.4221005
## 3         E3   2.2399873
## 171   G4:G11   1.4988470
## 229   G8:G15  -1.1335375
## 48     E2:E3   0.5565101
## 273  G13:G14   0.5321984

I wanted to include another method to further confirm this analysis, so I implemented stepwise regression. I found that after imputing my transformed model, the stepwise regression model summary showed me that E2:E3, E3, and G12:G14 were my strongest predictors. Interestingly, E2 was again not listed on its own. I decided to not include any variables beyond this as the adjusted R2 did not increase at a high enough rate to warrant it and the Bayesian Information Criterion (BIC) did not decrease at a substantial rate after E3 was included.

Model Summary
model adjR2 BIC
(Intercept) + E2:E3 0.368500782390762 -603.917263550716
(Intercept) + E2:E3 + G12:G14 0.397029060423632 -659.800268956789
(Intercept) + E3 + E2:E3 + G12:G14 0.407085508273982 -676.187882510353
(Intercept) + E3 + E2:E3 + G4:G11 + G12:G14 0.410925153232056 -678.713977292598
(Intercept) + E3 + E2:E3 + G4:G11 + G9:G14 + G12:G14 0.413506981159028 -678.414535578878

I then measured them on their own, by splitting any interaction terms, to ensure they were all significant and to see if I had overlooked any terms. They all proved to be significant with t-values above 4.5 and p-values near 0. E2 and E3 were both significant on their own with t-values above 17 despite E2 not being listed on its own in either of the variable selection tests. Before putting my variables back together and testing, I cross-checked my strongest variables from the stepwise regression with the ones from my ‘LASSO’ method. I found that both contained variables E2:E3, E3, and G12:G14 together. I decided to only put the G variables back together as their t-values together were higher than separated, and they were listed as interaction terms in both variable selection tests. I decided to keep E2 and E3 separated as when I split them, their t-values were significant and their p-values approached 0. Putting my variables back together where I saw fit, I found the t-values to all be above 8 with E3 proving to be the highest at 23.98.

Estimate Std. Error t value Pr(>|t|)
(Intercept) 46.561903 3.0772104 15.131205 0
E2 4.947677 0.2771651 17.851007 0
E3 6.697147 0.2792249 23.984780 0
G12:G14 6.931125 0.8358687 8.292122 0

Results

After much trial and error, my final model was Y1.9 = 46.562 + 4.948E2 + 6.697E3 + 6.931G12:G14 + \(\epsilon\). This final model had the highest adjusted R2 at 0.4084 and also was the most normally distributed based on analysis from the Q-Q plot. In addition, the residual plot had no pattern.

I also used an ANOVA table to show these variables were significant. They all had large sum of squares showing they explained the variance. They also had high F-values all above 65 with p-values near 0.

ANOVA Table

## Analysis of Variance Table
## 
## Response: I(Y^1.9)
##             Df Sum Sq Mean Sq F value    Pr(>F)    
## E2           1  63866   63866 298.837 < 2.2e-16 ***
## E3           1 120045  120045 561.707 < 2.2e-16 ***
## G12:G14      1  14695   14695  68.759 2.678e-16 ***
## Residuals 1339 286164     214                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

As I conducted various tests at the same time, I accounted for multiple comparisons using the Holm-Bonferroni method. This is a stepwise procedure that is less conservative then the Bonferonni method, reducing risk for type 2 errors while still controlling for type 1. This adjusted my p-values and showed again that E2, E3, and my interaction variable G12:G14 had adjusted p-values close to 0, showing they are all statistically significant. I also looked for significant predictors with the threshold being a p-value less than 0.05 to keep the odds of a type 1 error low. The second output below shows the chosen significant predictors that fall in that threshold.

##   (Intercept)            E2            E3       G12:G14 
##  1.443693e-47  1.258880e-63 1.902618e-105  2.677973e-16
## [1] "(Intercept)" "E2"          "E3"          "G12:G14"

Limitations

Despite coming to all of these conclusions and answering my main question to some extent, I had some limitations in my analysis. I struggled to understand the relationship of some of my significant independent variables to the dependent variable. I kept finding that E2 and E3 had some kind of nonlinear relationship, yet every attempt at showing that in a model only worsened the fit. I was able to increase my adjusted R2 using more interaction terms and combing the E2 and E3 variables, but once adjusted p-values were introduced, no predictors were meeting the minimum threshold except for the predictors in my final model. Although I had some limitations, I found a model that fit better than the Environmental variable model and the model containing all variables. Overall this model may be missing some nonlinear relationships, but the significant variables are all accounted for.

Appendix

Code

library(dplyr)
library(ggplot2)
library(tidyverse)
read_csv("P2_849114.csv")
data <- read.csv("P2_849114.csv")
head(data)
glimpse(data)
summary(data)

#Scatter plot of Y vs Ei
ggplot(data, aes(x = E1, y = Y)) + geom_point() + ggtitle("Y vs E1") + geom_smooth(method = "loess")
ggplot(data, aes(x = E2, y = Y)) + geom_point() + ggtitle("Y vs E2") + geom_smooth(method = "loess")
ggplot(data, aes(x = E3, y = Y)) + geom_point() + ggtitle("Y vs E3") + geom_smooth(method = "loess")
ggplot(data, aes(x = E4, y = Y)) + geom_point() + ggtitle("Y vs E4") + geom_smooth(method = "loess")


lm_E1 <- lm(Y ~ (E1), data)
lm_E2 <- lm(Y ~ (E2), data)
lm_E3 <- lm(Y ~ (E3), data)
lm_E4 <- lm(Y ~ (E4), data)
summary(lm_E1)
summary(lm_E2)
summary(lm_E3)
summary(lm_E4)



library(mgcv)
gam_model <- gam(Y ~ s(E1) + s(E2) + s(E3) + s(E4), data = data)
summary(gam_model)

# Checking for polynomial terms
lm_E1_poly <- lm(Y ~ I(E1^2), data = data)
lm_E2_poly <- lm(Y ~ I(E2^2), data = data)
lm_E3_poly <- lm(Y ~ I(E3^2), data = data)
lm_E4_poly <- lm(Y ~ I(E4^2), data = data)

summary(lm_E1_poly)
summary(lm_E2_poly)
summary(lm_E3_poly)
summary(lm_E4_poly)

#Correlation of DV with IVs
correlations <- data %>%
  dplyr::select(Y, starts_with("E"), starts_with("G")) %>%
  cor(use = "pairwise.complete.obs")
View(correlations)
sorted_correlations <- sort(abs(correlations["Y", -1]), decreasing = TRUE)
View(sorted_correlations)
head(sorted_correlations, 
     n = 4)


#Linear Model with only Environmental Variables
lm_env <- lm(Y ~ E1+E2+E3+E4, data)
summary(lm_env)
plot(lm_env)

#Linear Model with all Variables (Environmental and Genetic)
lm_all <- lm(Y ~ (E1+E2+E3+E4+G1+G2+G3+G4+G5+G6+G7+G8+G9+G10+
               G11+G12+G13+G14+G15+G16+G17+G18+G19+G20)^2, data)
summary(lm_all)
plot(lm_all) 

library(MASS)
boxcox(lm_all, lambda = seq(-4, 4, by = 0.1))

lm_all_trans <- lm(I(Y^1.9) ~ (E1+E2+E3+E4+G1+G2+G3+G4+G5+G6+G7+G8+G9+G10+
                                  G11+G12+G13+G14+G15+G16+G17+G18+G19+G20)^2, data)
summary(lm_all_trans)
plot(lm_all_trans)


#Lasso Method 
library(glmnet)
x <- model.matrix(Y ~ (E1+E2+E3+E4+G1+G2+G3+G4+G5+G6+G7+G8+G9+G10+
               G11+G12+G13+G14+G15+G16+G17+G18+G19+G20)^2, data)[,-1]
y <- I(data$Y^1.9)
set.seed(315)
cv.lasso <- cv.glmnet(x, y, alpha = 1)
lasso.best <- glmnet(x, y, alpha = 1, lambda = cv.lasso$lambda.min)
coef(lasso.best)

coefs <- coef(lasso.best, s = cv.lasso$lambda.min)
coefs_vector <- as.numeric(coefs[-1, ])
variable_names <- rownames(coefs)[-1]

coefs_df <- data.frame(variable = variable_names, coefficient = coefs_vector)
sorted_coefs <- coefs_df[order(-abs(coefs_df$coefficient)), ]

print(sorted_coefs)

#Stepwise Regression
library(leaps)
model <- regsubsets( model.matrix(lm_all_trans)[,-1], I((data$Y)^1.9),
                     nbest = 1 , nvmax=5, 
                     method = 'forward', intercept = TRUE )
temp <- summary(model)

library(knitr)
predictors <- colnames(model.matrix(lm_all_trans))
lm_select <- apply(temp$which,
                   1, function(x) paste0(predictors[x], collapse = " + "))
kable(data.frame(cbind( model = lm_select, adjR2 = temp$adjr2, BIC = temp$bic)),
      caption='Model Summary')

#Splitting interaction terms
lm_main_effects <- lm( I(Y^1.9) ~ E1+E2+E3+E4+G1+G2+G3+G4+G5+G6+G7+G8+G9+G10+G11+G12+G13+G14+G15+G16+G17+G18+G19+G20, data)
temp_me <- summary(lm_main_effects)
kable(temp_me$coefficients[ abs(temp_me$coefficients[,4]) <= 0.001, ], caption='Sig Coefficients')

#Putting them back together
model_2 <- lm( I(Y^1.9) ~ (E2+E3+G12:G14), data)
temp_m2 <- summary(model_2) 
kable(temp_m2$coefficients[ abs(temp_m2$coefficients[,"t value"]) >= 4, ])

#Final Model
post_model<- lm( I(Y^1.9) ~ E2+E3+G12:G14, data)
summary(post_model)
plot(post_model)

lm_poly_E2 <- lm(I(Y^1.9) ~ poly(E2, 2, raw = TRUE) + E3 + G12:G14, data = data)
lm_poly_E3 <- lm(I(Y^1.9) ~ E2 + poly(E3, 2, raw = TRUE) + G12:G14, data = data)
lm_poly_E2E3 <- lm(I(Y^1.9) ~ poly(E2, 2, raw = TRUE) + poly(E3, 2, raw = TRUE) + G12:G14, data = data)

summary(lm_poly_E2)
summary(lm_poly_E3)
summary(lm_poly_E2E3)

plot(lm_poly_E2)
plot(lm_poly_E3)
plot(lm_poly_E2E3)

anova(post_model)

poly_check_model<- lm( I(Y^1.9) ~ E2+E3+I(E2^2)+I(E3^2)+G12:G14, data)
summary(poly_check_model)
plot(poly_check_model)

library(stats)

mc_model <- lm(I(Y^1.9)~ E2:E3+G12:G14 , data )
summary(mc_model)
p_values <- summary(mc_model)$coefficients[, "Pr(>|t|)"]

#Bonferroni method
p_adjusted <- p.adjust(p_values, method = "holm")

print(p_adjusted)


threshold <- 0.05  
significant_predictors <- names(p_adjusted)[p_adjusted < threshold]

print(significant_predictors)

mc_model_final <- lm(I(Y^1.9)~ E2+E3+G12:G14 , data)

p_values_final <- summary(mc_model_final)$coefficients[, "Pr(>|t|)"]

p_adjusted_final <- p.adjust(p_values_final, method = "holm")
threshold <- 0.05  
significant_predictors_final <- names(p_adjusted_final)[p_adjusted_final < threshold]

print(significant_predictors_final)