pdf("C:\\Rmisc\\graphs\\cholest.pdf") dLDL <- c(-35.8,-45.0,-52.7,-49.7,-58.2,-66.0) DOSE <- c(1,2.5,5,10,20,40) r <- c(15,17,12,14,18,13) lnDOSE <- log(DOSE) cholest <- data.frame(dLDL,lnDOSE,r) attach(cholest) cholest.wls <- lm(dLDL ~ lnDOSE, weights=r) summary(cholest.wls) anova(cholest.wls) x <- seq(0,4,0.01) y <- predict(cholest.wls,list(lnDOSE=x)) plot(lnDOSE,dLDL,main="Weighted Least Squares Regression of dLDL on ln(DOSE)", xlab="ln(DOSE)",ylab="dLDL",pch=16) lines(x,y) y <- as.matrix(dLDL) x0 <- matrix(rep(1,nrow(y)),ncol=1) x1 <- as.matrix(lnDOSE) x <- cbind(x0,x1) w<- diag(sqrt(r)) xs <- w %*% x; ys <- w %*% y; betahatw <- solve(t(xs) %*% xs) %*% t(xs) %*% ys es <- ys - xs %*% betahatw ssew <- t(es) %*% es msew <- ssew/(6-2) vbetahatw <- diag(msew[1,1] * solve(t(xs) %*% xs)) sebetahatw <- sqrt(vbetahatw) yhatw <- x %*% betahatw ew <- y-yhatw print(cbind(betahatw,sebetahatw)) print(cbind(yhatw,ew)) dev.off()