# David Wells, BIFX 553, 2/16/17 homework

# For dat1, I decided that the model works best using untransformed values
# and witout the x1*x2 interaction. Several models with transformed values
# and without the x1*x2 interaction looked good at first but turned out to be
# likely heteroscedastic. The x1*x2 interactions seemded to make the q-q
# plots less linear.


# Load the packages and data
library(broom)
library(car)
library(tidyverse)
library(splines)

load('~/BIFX553/Feb16hw.RData')

dat1.lm <- lm(y ~ x1 + x2, data=dat1)

# Model Summary
tidy(dat1.lm)

# The model isn't improved by included by including x1*x2.
dat1_xinter.lm <- lm(y ~ x1 + x2 + x1*x2, data=dat1)
anova(dat1.lm, dat1_xinter.lm)

# The model doesn't show significant autocorrelation.
durbinWatsonTest(dat1.lm)

# The model doesn't show significant heteroscedasticity.
ncvTest(dat1.lm)

# Q-q plot is linear.
qqp(dat1.lm)

# This model shows only one outlier.
outlierTest(dat1.lm)

# The variance influence factors are lower than other models'
vif(dat1.lm)


# For dat2 I took the log10 of the y values and chose the model log10(y) ~ x1.
# I tried several models but they were heteroscedastic and their C+R plots weren't
# linear like the model I chose. My model has only one outlier.

dat2 <- mutate(dat2, ly = log10(y))
dat2.lm <- lm(ly ~ x1,data=dat2)
tidy(dat2.lm)
crPlots(dat2.lm)
durbinWatsonTest(dat2.lm)
ncvTest(dat2.lm)
outlierTest(dat2.lm)


# For dat3 I went with untransformed values but the C+R plot suggests that
# a non-paremetric regression line will fit the data. I used spline to
# plot the line. The model satisfies our other assumptions.

dat3.lm <- lm(y~x1,data=dat3)

tidy(dat3.lm)
crPlots(dat3.lm)
durbinWatsonTest(dat3.lm)
ncvTest(dat3.lm)
outlierTest(dat3.lm)

# Plot the line.

dat3_sp <- ggplot(dat3.lm, aes(x = x1, y = y)) +
  geom_point() +
  geom_smooth(method = "loess")

dat3_sp + geom_smooth(method = 'lm', se = FALSE, color = 'orange3',
  linetype = 2, formula = y ~ bs(x, knots = 0, degree = 1))