## Paul Johnson
## 2013-04-01
## regression-doublelog-1.R
## Question. What do you get when you "log both sides" of a
## regression model?
## In economics and biology, the double-log model is very
## common. It describes an interactive, multiplicative process.
##
## With just one predictor, the theory is
## y = b0 * x^b1 * exp(e)
## where I wish I could add subscripts on y, x, and e.
## Using the rules of logs, that simplifies to the linear
## model. If there were predictors x1,x2, and so forth,
## they would be added as additional multiplicative terms
## with exponents.
## Recall the log laws
## 1. log(x^b1) = b1 * log(x)
## 2. log(c*d) = log(c) + log(d)
set.seed(123123)
dat <- data.frame( x = runif(1000))
## The 3 parameters are b0, b1, and stde
b0 <- 0.1
b1 <- 1.2
stde <- 2
dat$y <- b0 * dat$x^(b1) * exp(rnorm(1000, m = 0, s = stde))
plot(y ~ x, data = dat, main="No Apparent Relationship?")
## Interesting. Looks like nothing.
m0 <- lm(y ~ x, data = dat)
summary(m0)
## One of the usual transformations we try is the log, which
## is justified either by idea of changing a "skewed" variable to
## a more normal one, or by the desire to fit the interactive
## model above.
## Test out plot's built-in antilogger
plot(y ~ x, data = dat, log = "xy")
## Better fix the labels, urgently
plot(y ~ x, data = dat, log = "xy", xlab = "log(x)", ylab = "log(y)")
## Previous same as this, this is more usual to the way I would do this:
plot(log(y) ~ log(x), data = dat)
m1 <- lm(log(y) ~ log(x), data = dat)
summary(m1)
abline(m1, col = "red")
## The challenge of interpretation is to get predicted
## values in a scale that interests us.
## I notice termplot does the same thing that rockchalk::plotCurves
## will do
termplot(m1, partial = TRUE, se = TRUE)
library(rockchalk)
## needs version 1.7.2 or newer
m1pc <- plotCurves(m1, plotx = "x", interval = "conf")
## if your rockchalk is not on the development path, can still
## construct that manually. Here's how
nd <- data.frame(x = plotSeq(dat$x, 40))
pdat <- predict(m1, newdata = nd, interval = "conf")
head(pdat)
nd <- cbind(nd, pdat)
nd
plot(log(y) ~ x, data = dat, col = gray(.7))
lines(fit ~ x, data = nd, col = "red")
lines(lwr ~ x, data = nd, col = "red", lty = 2)
lines(upr ~ x, data = nd, col = "red", lty = 2)
## I really want to plot this in the (x,y) plane
## not the log(x) or log(y). Here's an easy way
## using the results from before.
nd$ypred <- exp(nd$fit)
head(nd)
titl <- "double-log regression in x,y space"
plot(y ~ x, data = dat, col = gray(.8), main = titl)
lines(ypred ~ x, data = nd, col = "red", lwd = 2)
## The student can go back to the beginning and
## change the coefficients to see how these plots
## might change.