# Linear Models
#--------------

# Simulated Example
#------------------

y <- rep(1:4, each=25)     # note this is numeric, not a factor
z <- 5 + 2*y + rnorm(100)  # first response var
a <- rep(1:5, times=20)
za <- 4 + 3*y +2*a +rnorm(100) # a second response var

table(y,a)  # balanced simulation

ds <- data.frame(y,a,z,za)

# Some model specification
#-------------------------
(reg1 <- lm(z ~ 1 + y + a, data=ds)) # Explicit intercept
lm(z ~ y + a, data=ds)               # Implicit intercept
lm(z ~ a + 1 + y, data=ds)           # Somewhat arbitrary order

(reg2 <- lm(z ~ 0 + y + a, data=ds)) # Zero intercept

lm(z ~ 5 + y + a, data=ds) # ERROR, 5 is not a valid term

ds$int <- 5  # recycling!
lm(z ~ offset(int) + y + a, data=ds)     # Offset the model
lm(z ~ offset(int) + 0 + y + a, data=ds) # Offset the zero dose model


lm(z ~ 1 + y + a - 1, data=ds)           # Remove a term

lm(z ~ (1 + y) + a, data=ds)  # Group terms, visually

lm(z ~ I(1 + y) + a, data=ds) # A calculated term, "protected"
lm(z ~ I(y^2), data=ds)       # another
lm(log(z) ~ y + a, data=ds)   # another

lm(z ~ 1 + (y + a), data=ds)  # Another grouping
lm(z ~ 1 + I(y + a), data=ds) # Another calculated term

lm(z ~ 1 + 7*(y + a), data=ds)    # ERROR, confuses crossing and multiplication
lm(z ~ 1 + I(7*(y + a)), data=ds) # Calculated term
lm(z ~ y*a, data=ds)              # crossed term
lm(z ~ y + a + y:a, data=ds)      # explicit main and interaction terms

m <- formula(z~y)      # formulas can be objects
str(m)
lm(m) 
plot(m)

# Extracting information from a model fit
#----------------------------------------
reg1 <- lm(z ~ y, data=ds)
str(reg1)
summary(reg1)
anova(reg1)

model1 <- lm(za ~ y, data=ds)
model2 <- lm(za ~ a + y, data=ds)

anova(model1, model2)
 
coef(model1)
plot(za ~ y)
abline(coef(model1))

confint(model2)

####################
# Example Regression
####################
library(lattice)
xyplot(dist ~ speed, data=cars, type=c("p", "r"))

regcar <- lm(dist ~ speed, data=cars)
summary(regcar)

# Simulated ANOVA
#----------------
yf <- as.factor(y)

bwplot(z~yf)
aov1 <- lm(z ~ yf)
summary(aov1)
anova(aov1)
anova(reg1)

xyplot(residuals(reg1)~predict(reg1))
xyplot(residuals(aov1)~predict(aov1))

#############
# Example AOV
#############
aov2 <- lm(weight ~ feed, chickwts)
summary(aov2)
plot(chickwts)
xyplot(weight ~ feed, data=chickwts)
xyplot(residuals(aov2) ~ predict(aov2), group=chickwts$feed, 
    auto.key=TRUE, xlab="Predicted values", ylab="Residuals")