Large-sample tests based on quadratic forms T2 <- function(X,A,h) # Test H0: A mu = h # X is n by p data matrix, n iid vectors, length p # F for multivariate normal, chisq for n large { h <- cbind(h) ; d <- length(h) ; n <- nrow(X) if(d != nrow(A))stop("h and A incompatible") if(ncol(X) != ncol(A))stop("X and A incompatible") xbar <- cbind(apply(X,2,mean)) ; S <- var(X) Tsq <- n * ( t(A%*%xbar - h) %*% solve(A %*% S %*% t(A) ) %*% (A%*%xbar - h) ) F <- Tsq * (n-d)/(n*d-d) ; pF <- 1-pf(F,d,n-d) pChi <- 1-pchisq(Tsq,d) T2 <- cbind(numeric(4)) rownames(T2) <- c("T-squared","F","P-value for F", "P-value for Chiquare") T2[1] <- Tsq ; T2[2] <- F T2[3] <- pF ; T2[4] <- pChi T2 } # End of function T2 > sweat <- read.table("sweat.dat") > T2(sweat,diag(3),c(4,50,10)) [,1] T-squared 9.73877286 F 2.90454629 P-value for F 0.06492834 P-value for Chiquare 0.02092230 > > noise1 <- read.table("noise1.dat") > noise <- as.matrix(noise1[,2:6]) > contr <- rbind(c(1, -1, 0, 0, 0), + c(0, 1, -1, 0, 0), + c(0, 0, 1, -1, 0), + c(0, 0 , 0, 1, -1)) > T2(noise1,contr,c(0,0,0)) # Pretend to make a mistake Error in T2(noise1, contr, c(0, 0, 0)) : h and A incompatible > T2(noise,contr,c(0,0,0,0)) [,1] T-squared 6.453032e+01 F 1.531228e+01 P-value for F 1.614600e-08 P-value for Chiquare 3.231859e-13 > ------- Now a version in which you supply the covariance matrix and mean. It applies to indepenent groups too. T3 <- function(xbar,Sn,A,h) # Test H0: A mu = h # xbar is vector of sample means, Sn estimated # covariance matrix of xbar { xbar <- cbind(xbar) ; h <- cbind(h) ; d <- length(h) if(d != nrow(A))stop("h and A incompatible") if(ncol(A) != ncol(Sn))stop("Sn and A incompatible") Tsq <- ( t(A%*%xbar - h) %*% solve(A %*% Sn %*% t(A) ) %*% (A%*%xbar - h) ) pChi <- 1-pchisq(Tsq,d) T3 <- cbind(numeric(3)) rownames(T3) <- c("Chi-squared","df","P-value") T3[1] <- Tsq ; T3[2] <- d ; T3[3] <- 1-pchisq(Tsq,d) T3 } # End of function T3 noise <- as.matrix(noise1[,2:6]) contr <- rbind(c(1, -1, 0, 0, 0), c(0, 1, -1, 0, 0), c(0, 0, 1, -1, 0), c(0, 0 , 0, 1, -1)) Sn <- var(noise)/n nmean <- apply(noise,2,mean) > Sn <- var(noise)/n > nmean <- apply(noise,2,mean) > nmean V2 V3 V4 V5 V6 39.82167 36.83000 35.30167 34.38500 31.44500 T3(nmean,Sn,contr,c(0,0,0,0)) > T3(nmean,Sn,contr,c(0,0,0,0)) [,1] Chi-squared 6.453032e+01 P-value 3.231859e-13 ------- Simulated market research data for 5 companies. For each company, assess: Unaided & Aided Brand & Advertising awareness. That's 20 vars. Two waves. UNAIDED AD AWARENESS AIDED AD AWARENESS UNAIDED BRAND AWARENESS UNAIDED BRAND AWARENESS We have two big files. X1 has n1 rows and 20 columns, X2 has n2 rows and 20 columns. Here are the results in tabular form. Wave1 Wave2 UA 26.466667 21.6438356 AA 43.733333 39.5890411 UB 23.066667 18.2876712 AB 57.133333 53.6986301 UA 35.133333 30.0684932 AA 44.066667 38.8356164 UB 52.466667 45.2739726 AB 57.733333 53.2191781 UA 13.600000 4.3835616 AA 31.933333 27.9452055 UB 20.800000 18.6986301 AB 44.000000 0.4109589 UA 25.066667 3.9726027 AA 42.066667 33.6986301 UB 45.666667 36.5753425 AB 57.133333 52.9452055 UA 5.000000 2.8767123 AA 5.333333 4.2465753 UB 8.400000 5.8219178 AB 12.666667 7.7397260 xbar1 <- apply(X1,2,mean) ; S1 <- var(X1) xbar2 <- apply(X2,2,mean) ; S2 <- var(X2) > n1 ; n2 [1] 1500 [1] 1460 > cbind(xbar1,xbar2)*100 xbar1 xbar2 [1,] 26.466667 21.6438356 [2,] 43.733333 39.5890411 [3,] 23.066667 18.2876712 [4,] 57.133333 53.6986301 [5,] 35.133333 30.0684932 [6,] 44.066667 38.8356164 [7,] 52.466667 45.2739726 [8,] 57.733333 53.2191781 [9,] 13.600000 4.3835616 [10,] 31.933333 27.9452055 [11,] 20.800000 18.6986301 [12,] 44.000000 0.4109589 [13,] 25.066667 3.9726027 [14,] 42.066667 33.6986301 [15,] 45.666667 36.5753425 [16,] 57.133333 52.9452055 [17,] 5.000000 2.8767123 [18,] 5.333333 4.2465753 [19,] 8.400000 5.8219178 [20,] 12.666667 7.7397260 Sn <- diag(40) # 40 by 40 Identity Sn[1:20,1:20] <- S1/n1 ; Sn[21:40,21:40] <- S2/n2 # Yielding block diagonal structure save(n1,xbar1,S1,n2,xbar2,S2,Sn,file="marketresearch.rsave") # retrieve with load Now test if any change between wave 1 and 2 ------------------------------------------- XBAR <- c(xbar1,xbar2) cont <- numeric(20*40) ; dim(cont) <- c(20,40) cont[1:20,1:20] <- diag(20) ; cont[1:20,21:40] <- -diag(20) T3(XBAR,Sn,cont,numeric(20)) > T3(XBAR,Sn,cont,numeric(20)) [,1] Chi-squared 1560.319 df 20.000 P-value 0.000 ########################################################################## # Delta Method: Is % change in advertising awareness different from % # change in brand awareness? Test aided and unaided # awareness simultaneously. Just Brand A. ########################################################################## gxbar <- c((XBAR[21]-XBAR[1])/XBAR[1],(XBAR[22]-XBAR[2])/XBAR[2], (XBAR[23]-XBAR[3])/XBAR[3],(XBAR[24]-XBAR[4])/XBAR[4]) names(gxbar) <- c("Unaided Ad","Aided Ad","Unaided Brand","Aided Brand") # Now make gdot(muhat) block1 <- diag(c(-XBAR[21]/XBAR[1]^2,-XBAR[22]/XBAR[2]^2, -XBAR[23]/XBAR[3]^2,-XBAR[24]/XBAR[4]^2)) block2 <- diag(c(1/XBAR[1],1/XBAR[2],1/XBAR[3],1/XBAR[4])) z16 <- numeric(4*16) ; dim(z16) <- c(4,16) gdot <- cbind(block1,z16,block2,z16) Vn <- gdot %*% Sn %*% t(gdot) gcontrast <- rbind(c(1,0,-1,0),c(0,1,0,-1)) gxbar T3(gxbar,Vn,gcontrast,c(0,0)) > gxbar Unaided Ad Aided Ad Unaided Brand Aided Brand -0.18222284 -0.09476278 -0.20718188 -0.06011733 > > T3(gxbar,Vn,gcontrast,c(0,0)) [,1] Chi-squared 1.0618616 df 2.0000000 P-value 0.5880574