# To get basically the same data set as described on pages 10-13 of textbook: library(MASS) set.seed(1203) mymat1 <- mvrnorm(200, mu=c(0,0), Sigma=matrix(c(1,0.5,0.5,1),ncol=2)) mymat2 <- log(abs(mymat1)) + 12 # mymat2 can be thought of as the observed data matrix which is not multivariate normal. chisplot(mymat2) boxcox(lm(mymat2[,1]~1),lambda=seq(-10,10,1/5)) boxcox(lm(mymat2[,2]~1),lambda=seq(-10,10,1/5)) # To see the actual lambda values: box.x <- boxcox(lm(mymat2[,1]~1),lambda=seq(-10,10,1/5),plotit=F)\$x box.y <- boxcox(lm(mymat2[,1]~1),lambda=seq(-10,10,1/5),plotit=F)\$y cbind(box.x, box.y) box.x <- boxcox(lm(mymat2[,2]~1),lambda=seq(-10,10,1/5),plotit=F)\$x box.y <- boxcox(lm(mymat2[,2]~1),lambda=seq(-10,10,1/5),plotit=F)\$y cbind(box.x, box.y) # Good choices for the lambdas: lambda=c(5.5, 6) mymat3 <- cbind( ((mymat2[,1])^5.5 - 1)/5.5, ((mymat2[,2])^6 - 1)/6 ) # mymat3 is the transformed data matrix #### A more efficient way to do the transformations of each column of mymat2. #### This assumes all the lambda_j's are in a vector called "lambda", as above. newmat <- NULL for (j in 1:length(lambda)){ if (lambda[j]!=0){ newcol<- (airpol.data.sub[,j]^lambda[j] - 1)/lambda[j] } else { newcol<-log(airpol.data.sub[,j]) } newmat <- cbind(newmat,newcol) } mymat3 <- newmat colnames(mymat3) <-colnames(mymat2) # mymat3 is the transformed data matrix #### #### chisplot(mymat3)