biod <- read.csv("http://www.stat.ufl.edu/~winner/data/biodiesel_transest.csv") attach(biod); names(biod) head(biod) pairs(biod[,c(2:5,7,9)]) n <- length(run_id) Y <- cbind(meth_pct,eth_pct,prop1_pct) Z <- cbind(rep(1,n),b_temp,b_pres,b_time) r <- ncol(Z) - 1 m <- ncol(Y) (Beta_hat <- solve(t(Z) %*% Z) %*% t(Z) %*% Y) cor(Y) Yhat <- Z %*% Beta_hat Eps <- Y - Yhat (SSCP_E <- t(Eps) %*% Eps) (S_E <- (1/(n-(r+1))) * SSCP_E) (Sigma_hat <- (1/n) * SSCP_E) (V.Beta_hat1 <- S_E[1,1] * solve(t(Z) %*% Z)) (V.Beta_hat2 <- S_E[2,2] * solve(t(Z) %*% Z)) (V.Beta_hat3 <- S_E[3,3] * solve(t(Z) %*% Z)) (COV.B1.B2 <- S_E[1,2] * solve(t(Z) %*% Z)) (COV.B1.B3 <- S_E[1,3] * solve(t(Z) %*% Z)) (COV.B2.B3 <- S_E[2,3] * solve(t(Z) %*% Z)) z0 <- c(1, 325, 9, 37.5) beta.1 <- Beta_hat[,1] se.beta.1 <- sqrt(diag(V.Beta_hat1)) t.beta.1 <- beta.1 / se.beta.1 pv.beta.1 <- 2*(1-pt(abs(t.beta.1),n-(r+1))) out.beta1 <- cbind(beta.1,se.beta.1,t.beta.1,pv.beta.1) colnames(out.beta1) <- c("Estimate","Std Err","t", "P-value") round(out.beta1,4) predict.z0 <- t(t(z0) %*% Beta_hat) se.mean.z0 <- matrix(rep(0,3),ncol=1) se.pred.z0 <- matrix(rep(0,3),ncol=1) for (i in 1:3) { se.mean.z0[i] <- sqrt(S_E[i,i] * (t(z0) %*% solve(t(Z) %*% Z) %*% z0)) se.pred.z0[i] <- sqrt(S_E[i,i] * (1 + t(z0) %*% solve(t(Z) %*% Z) %*% z0)) } ci.crit <- sqrt((m*(n-r-1)/(n-r-m))*qf(.975,m,n-r-m)) mean.ci.lo <- predict.z0-ci.crit*se.mean.z0 mean.ci.hi <- predict.z0+ci.crit*se.mean.z0 pred.pi.lo <- predict.z0-ci.crit*se.pred.z0 pred.pi.hi <- predict.z0+ci.crit*se.pred.z0 pred.out <- cbind(predict.z0,mean.ci.lo,mean.ci.hi,pred.pi.lo,pred.pi.hi) colnames(pred.out) <- c("Predicted","Mean LB","Mean UB","Pred LB","Pred UB") round(pred.out,2) ### Test H0: Beta_Time = 0 Z1 <- Z[,1:3] q <- ncol(Z1) - 1 (Beta_hat1 <- solve(t(Z1) %*% Z1) %*% t(Z1) %*% Y) Yhat1 <- Z1 %*% Beta_hat1 Eps1 <- Y - Yhat1 SSCP_E1 <- t(Eps1) %*% Eps1 (Sigma_hat1 <- (1/n) * SSCP_E1) ## Chi-Square Test Lambda.2n <- (det(Sigma_hat) / det(Sigma_hat1)) X2.TS <- -(n-r-1-0.5*(m-r+q+1))*log(Lambda.2n) X2.CV <- qchisq(.95,m*(r-q)) X2.PV <- 1-pchisq(X2.TS,m*(r-q)) X2.out <- cbind(X2.TS, X2.CV, X2.PV) colnames(X2.out) <- c("Test Stat", "X2(.05)", "P-value") round(X2.out, 3) ## F-Tests L <- matrix(c(0,0,0,1),ncol=4) # L <- matrix(c(0,0,1,0),ncol=4) # L <- matrix(c(0,1,0,0),ncol=4) H <- t(L%*%Beta_hat) %*% solve(L %*% solve(t(Z)%*%Z) %*% t(L)) %*% (L%*%Beta_hat) E <- SSCP_E H1 <- L %*% solve(t(Z)%*%Z) %*% t(L) Lambda <- det(E)/det(H+E) ### Wilk's Lambda V <- sum(diag(H %*% solve(H+E))) ### Pillai's Trace q1 <- qr(H+E)$rank q2 <- 1 # rank(H1)$rank nu <- n-(r+1) s <- min(q1,q2) m1 <- 0.5*(abs(q1-q2)-1) m2 <- 0.5*(nu-q1-1) t.W <- ifelse(q1^2+q2^2-5>0,sqrt((q1^2*q2^2-4)/(q1^2+q2^2-5)),1) lambda <- nu - (q1-q2+1)/2 u <- (q1*q2-2) / 4 F.W <- ((1-Lambda^(1/t.W))/Lambda^(1/t.W)) * ((lambda*t.W - 2*u)/(q1*q2)) df1.W <- q1*q2 df2.W <- lambda*t.W - 2*u F.W.CV <- qf(.95,df1.W,df2.W) F.W.PV <- 1-pf(F.W,df1.W,df2.W) F.P <- (V/(s-V)) * ((2*m2+s+1)/(2*m1+s+1)) df1.P <- 2*m1+s+1 df2.P <- 2*m2+s+1 F.P.CV <- qf(.95,df1.P,df2.P) F.P.PV <- 1-pf(F.P,df1.P,df2.P) wilks.out <- cbind(F.W,df1.W,df2.W,F.W.CV,F.W.PV) colnames(wilks.out) <- c("Wilks F","df1","df2","F(.05)","P-value") round(wilks.out,4) pillais.out <- cbind(F.P,df1.P,df2.P,F.P.CV,F.P.PV) colnames(pillais.out) <- c("Pillais F","df1","df2","F(.05)","P-value") round(pillais.out,4) ### Regression based on lm function biod.mod1 <- lm(Y ~ b_temp + b_pres + b_time) summary(biod.mod1) summary(manova(biod.mod1),test="Wilks") biod.mod2 <- lm(Y ~ b_temp + b_pres) summary(biod.mod2) summary(manova(biod.mod2),test="Wilks") biod.mod3 <- lm(Y ~ b_temp + b_time) summary(biod.mod3) summary(manova(biod.mod3),test="Wilks") biod.mod4 <- lm(Y ~ b_pres + b_time) summary(biod.mod4) summary(manova(biod.mod4),test="Wilks") anova(biod.mod2, biod.mod1, test="Wilks") anova(biod.mod3, biod.mod1, test="Wilks") anova(biod.mod4, biod.mod1, test="Wilks")