## Known weights example
Data <- read.delim("http://web.as.uky.edu/statistics/users/pbreheny/760/data/statin.txt")
summary(fit <- lm(LDL~Statin, Data))
summary(fit.w <- lm(LDL~Statin, Data, weights=n))
sqrt(summary(fit)$r.squared)
sqrt(summary(fit.w)$r.squared)

## Hand speed
Data <- read.delim("http://web.as.uky.edu/statistics/users/pbreheny/760/data/handspeed.txt")
xyplot(Time~Age, Data, col="slateblue", alpha=.7, pch=19)

## Justifying the mean-variance relationship
M <- with(Data, by(Time, Age, mean))
SD <- with(Data, by(Time, Age, sd))
plot(M, SD, pch=19, xlab="Mean", ylab="SD") ## Roughly linear

## Linear fit
fit.w <- lm(Time~Age, Data)
for (i in 1:20) {
  w <- 1/fit.w$fitted.values^2
  fit.w <- lm(Time~Age, Data, weights=w)
}
fit <- lm(Time~Age, Data)
xx <- seq(min(Data$Age), max(Data$Age), len=99)
New <- data.frame(Age=xx)
with(Data, plot(Time~Age, pch=19, col="gray"))
col <- c("#FF4E37FF", "#008DFFFF")
lines(xx, predict(fit, New), type="l", col=col[1], lwd=3)
lines(xx, predict(fit.w, New), type="l", col=col[2], lwd=3)
legend("topleft", c("OLS", "WLS"), col=col, lwd=3)
summary(fit)
summary(fit.w)

## Bad idea
fit.w <- lm(Time~Age, Data)
for (i in 1:20) {
  w <- fit.w$fitted.values^2
  fit.w <- lm(Time~Age, Data, weights=w)
}
summary(fit.w)

## Quadratic fit
fit.w <- lm(Time~Age+I(Age^2), Data)
for (i in 1:20) {
  w <- 1/fit.w$fitted.values^2
  fit.w <- lm(Time~Age+I(Age^2), Data, weights=w)
}
fit <- lm(Time~Age+I(Age^2), Data)
with(Data, plot(Time~Age, pch=19, col="gray"))
lines(xx, predict(fit, New), type="l", col=col[1], lwd=3)
lines(xx, predict(fit.w, New), type="l", col=col[2], lwd=3)
legend("topleft", c("OLS", "WLS"), col=col, lwd=3)
summary(fit)
summary(fit.w)

## Changepoint at 60
f <- function(x, t) (x-t)*(x>t)
fit <- fit.w <- lm(Time~Age+f(Age,60), Data)
for (i in 1:20) {
  w <- 1/fit.w$fitted.values^2
  fit.w <- lm(Time~Age+f(Age,60), Data, weights=w)
}
pdf("handspeed3.pdf",height=5,width=10)
with(Data, plot(Time~Age, pch=19, col="gray"))
lines(xx, predict(fit, New), type="l", col=pal(2)[1], lwd=3)
lines(xx, predict(fit.w, New), type="l", col=pal(2)[2], lwd=3)
toplegend(c("OLS", "WLS"), col=pal(2), lwd=3, ncol=2)
dev.off()
require(visreg)
visreg(fit)
visreg(fit.w)
summary(fit)
summary(fit.w)