lm.r

m1 <- lm(I(log(price)) ~ I(carat^(1/3)), data = diamonds)
m2 <- update(m1, ~ . + carat)
m3 <- update(m2, ~ . + cut)
m4 <- update(m3, ~ . + color)
m5 <- update(m4, ~ . + clarity)
mtable(m1, m2, m3, m4, m5)


diamondsbig$logprice <- log(diamondsbig$price)
m1 <- lm(logprice ~ I(carat^(1/3)), diamondsbig[diamondsbig$price < 10000 & diamondsbig$cert == "GIA", ])
m2 <- update(m1, ~ . + carat)
m3 <- update(m2, ~ . + cut)
m4 <- update(m3, ~ . + color)
m5 <- update(m4, ~ . + clarity)
mtable(m1, m2, m3, m4, m5, sdigits = 3)


# Predictions
thisDiamond = data.frame(carat = 1.00, cut = "V.Good", color = "I", clarity="VS1")
modelEstimate = predict(m5, newdata = thisDiamond, interval="prediction", level = .95)


dat = data.frame(m4$model, m4$residuals) 

with(dat, sd(m4.residuals))

with(subset(dat, carat > .9 & carat < 1.1), sd(m4.residuals)) 

dat$resid <- as.numeric(dat$m4.residuals)
ggplot(aes(y = resid, x = round(carat, 2)), data = dat) + 
  geom_line(stat = "summary", fun.y = sd) 

lm_dplyr.r

# http://stackoverflow.com/a/33870137/3149679

library(dplyr)

d <- data.frame(
  state=rep(c('NY', 'CA'), c(10, 10)),
  year=rep(1:10, 2),
  response=c(rnorm(10), rnorm(10))
)

fitted_models = d %>% group_by(state) %>% do(model = lm(response ~ year, data = .))
fitted_models$model


library(broom)
fitted_models %>% tidy(model)
fitted_models %>% glance(model)
fitted_models %>% augment(model)