#
# quadratic-reg-var-R.doc
#

data()  # list data sets in R

data(Loblolly) 

plot(height~age, data=Loblolly)

lob.lmfit <- lm(height~ age + I(age^2), data=Loblolly)

plot(predict(lob.lmfit),resid(lob.lmfit))
abline(h=0)

qqnorm(resid(lob.lmfit))

# get VAR and N for "height" at each AGE

var.ht<- unlist(lapply(split(Loblolly$height,Loblolly$age),var))
var.ht
#        3         5        10        15        20        25 
#0.1628951 0.6651654 2.3650951 3.8057940 4.8921824 5.1476071 

rep.ht<- unlist(lapply(split(Loblolly$height,Loblolly$age),length))
#rep.ht
# 3  5 10 15 20 25 
#14 14 14 14 14 14 

Loblolly$wts<- rep(1/var.ht, rep.ht)

# wt=1/var(Ht) at each age
lob.wfit <- lm(height~ age + I(age^2), data=Loblolly, weights=wts) 

# wt=1/(predicted HT) with predictions from OLS
lob.wfit2 <- lm(height~ age + I(age^2), data=Loblolly, weights=1/predict(lob.lmfit))


# compare the SE from the different fits

summary(lob.lmfit)
summary(lob.wfit)
summary(lob.wfit2)


plot(height~age, data=Loblolly)




# 
#  output excerpts
#

Call:
lm(formula = height ~ age + I(age^2), data = Loblolly)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.7902 -1.1496  0.0183  0.9401  4.1815 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -7.607232   0.608017  -12.51   <2e-16 ***
age          3.959044   0.109236   36.24   <2e-16 ***
I(age^2)    -0.049838   0.003884  -12.83   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 1.702 on 81 degrees of freedom
Multiple R-squared: 0.9934,	Adjusted R-squared: 0.9932 
F-statistic:  6079 on 2 and 81 DF,  p-value: < 2.2e-16