sg_spread <- read.csv("http://users.stat.ufl.edu/~winner/sta4210/mydata/shotgun_spread.csv", header=T) attach(sg_spread); names(sg_spread) sg_spread1 <- sg_spread[cartridge==1,] detach(sg_spread); attach(sg_spread1) mean.spread <- as.numeric(tapply(sqrtA, range, mean)) sd.spread <- as.numeric(tapply(sqrtA, range, sd)) range1 <- seq(10,50,10) par(mfrow=c(2,1)) plot(range1,mean.spread) plot(range1,sd.spread) par(mfrow=c(1,1)) plot(range, sqrtA) ### Matrix Approach X <- cbind(rep(1,length(sqrtA)),range,range^2) W <- diag(1/SD_range^2) Y <- sqrtA (beta.ols <- solve(t(X) %*% X) %*% t(X) %*% Y) e.ols <- Y - X %*% beta.ols (sse.ols <- t(e.ols) %*% e.ols) (s2.ols <- sse.ols/(length(Y)-ncol(X))) (v.beta.ols <- s2.ols[1,1] * solve(t(X) %*% X)) se.beta.ols <- sqrt(diag(v.beta.ols)) t.beta.ols <- beta.ols/se.beta.ols (beta.wls <- solve(t(X) %*% W %*% X) %*% t(X) %*% W %*% Y) e.wls <- Y - X %*% beta.wls (sse.wls <- t(e.wls) %*% W %*% e.wls) (s2.wls <- sse.wls/(length(Y)-ncol(X))) (v.beta.wls <- s2.wls[1,1] * solve(t(X) %*% W %*% X)) se.beta.wls <- sqrt(diag(v.beta.wls)) t.beta.wls <- beta.wls/se.beta.wls print(round(cbind(beta.ols,se.beta.ols,t.beta.ols, beta.wls,se.beta.wls,t.beta.wls),5)) regweight <- 1/(SD_range^2) ### Ordinary Least Squares sg.mod1 <- lm(sqrtA ~ range + I(range^2)) summary(sg.mod1) #### Weighted Least Squares sg.mod2 <- lm(sqrtA ~ range + I(range^2), weight=regweight) summary(sg.mod2) vcov(sg.mod2) deviance(sg.mod2) anova(sg.mod2) ### Robust Standard Errors (beta.ols <- solve(t(X) %*% X) %*% t(X) %*% Y) e.ols <- Y - X %*% beta.ols E2 <- e.ols %*% t(e.ols) * diag(length(Y)) (v.beta.ols.robust<- solve(t(X) %*% X) %*% t(X) %*% E2 %*% X %*% solve(t(X) %*% X)) se.beta.ols.robust <- sqrt(diag(v.beta.ols.robust)) t.beta.ols.robust <- beta.ols/se.beta.ols.robust print(round(cbind(beta.ols,se.beta.ols.robust,t.beta.ols.robust),5))