Example R programs and commands 12. Multivariate scatterplots and MANOVAs # All lines preceded by the "#" character are my comments. # All other left-justified lines are my input. # All other indented lines are the R program output. # Generate some normally distributed random points: x <- rnorm(50,0,1) # 50 copies of N(0,1) y <- rnorm(50,0,10) # 50 copies of N(0,10) xy <- cbind(x,y) # 50 coordinate pairs (x,y) plot(xy) # Rotation of coordinates: rot <- matrix(c(1,-1,1,1),nrow=2,ncol=2) xy0 <- xy %*% rot # rotate by 45 degrees plot(xy0) # Repeat to get 50 more rotated and shifted bivariate normals: x <- rnorm(50,0,1) # 50 copies of N(0,1) y <- rnorm(50,0,10) # 50 copies of N(0,10) xy <- cbind(x,y) xy1 <- xy %*% rot + 1 # rotate by 45 degrees and shift by (1,1) plot(xy1) # Combine the two 50-point data sets and label them: xys<-rbind(xy0,xy1) hilo <- factor(gl(2,50,len=2*50)) # label "1" for xy0, "2" for xy1 blackred<-c(gl(2,50)) plot(xys, col=blackred) # Black for xy0, red for xy1 # Plot the means and ranges of the individual coordinates: plot(xys[,1], col=blackred) # first coordinate values plot(xys[,1]~hilo) # first coordinate box plot plot(xys[,2], col=blackred) # second coordinate values plot(xys[,2]~hilo) # second coordinate box plot # Accumulate the sums of squares: fit <- manova(xys~hilo) # ANOVA on the individual coordinates: summary.aov(fit) # MANOVA on the pairs of coordinates: summary.manova(fit) ################################################### # Single-factor MANOVA on tabulated data: # # Suppose we are given the following centipede # haemolymph amino acid data: # # Male centipedes: # # Alanine Aspartic Acid Tyrosine # ------- ------------- -------- # 7.0 17.0 19.7 # 7.3 17.2 20.3 # 8.0 19.3 22.6 # 8.1 19.8 23.7 # 7.9 18.4 22.0 # # Female centipedes: # # Alanine Aspartic Acid Tyrosine # ------- ------------- -------- # 7.3 17.4 22.5 # 7.7 19.8 24.9 # 8.2 20.2 26.1 # 8.3 22.6 27.5 # 6.4 23.4 28.1 # # Read the data from the table using "scan()": data=scan() 7.0 17.0 19.7 7.3 17.2 20.3 8.0 19.3 22.6 8.1 19.8 23.7 7.9 18.4 22.0 7.3 17.4 22.5 7.7 19.8 24.9 8.2 20.2 26.1 8.3 22.6 27.5 6.4 23.4 28.1 # # Each measurement is a row of 3 values, so make a matrix: # <> "ncol=3" since there are three coordinates per measurement # <> "byrow=TRUE" to fill the matrix by rows # <> "dimnames" is optional but it puts nice labels atop the columns; # the first NULL means just default labels on the rows. aacon<-matrix(data,ncol=3,byrow=TRUE,dimnames=list(NULL,c("Ala","Asp","Tyr"))) aacon Ala Asp Tyr [1,] 7.0 17.0 19.7 [2,] 7.3 17.2 20.3 [3,] 8.0 19.3 22.6 [4,] 8.1 19.8 23.7 [5,] 7.9 18.4 22.0 [6,] 7.3 17.4 22.5 [7,] 7.7 19.8 24.9 [8,] 8.2 20.2 26.1 [9,] 8.3 22.6 27.5 [10,] 6.4 23.4 28.1 # # Make factor labels for the single factor: # <> the first 5 triples are for males, the next 5 for females # <> encapsulating "gl()" within "factor()" is optional gender<-factor(gl(2,5,10,labels=c("Male","Female"))) gender [1] Male Male Male Male Male Female Female Female Female Female Levels: Male Female # # Function "manova()" computes the sums of squares and crossed # products and labels them for use in ANOVAs and MANOVAs: manova(aacon~gender) Call: manova(aacon ~ gender) Terms: gender Residuals Ala 0.016 3.320 Asp 13.689 29.000 Tyr 43.264 30.820 Deg. of Freedom 1 8 Residual standard error: 0.644205 1.903943 1.962779 Estimated effects may be unbalanced # Get a p-value for the effect of gender on the combination of 3 # measurements with the "summary()" function. If applied to the # output of "manova()" it uses Pillai's statistic by default. # summary(manova(aacon~gender)) Df Pillai approx F num Df den Df Pr(>F) gender 1 0.9146 21.4232 3 6 0.001317 ** Residuals 8 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 # Use Wilks', Hotelling-Lawley's, or Roy's statistics instead: summary(manova(aacon~gender),test="W") summary(manova(aacon~gender),test="H") summary(manova(aacon~gender),test="R") # Perform 3 ANOVAs on the individual responses: summary(aov(aacon~gender)) Response Ala : Df Sum Sq Mean Sq F value Pr(>F) gender 1 0.016 0.016 0.0386 0.8492 Residuals 8 3.320 0.415 Response Asp : Df Sum Sq Mean Sq F value Pr(>F) gender 1 13.689 13.689 3.7763 0.0879 . Residuals 8 29.000 3.625 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 Response Tyr : Df Sum Sq Mean Sq F value Pr(>F) gender 1 43.264 43.264 11.23 0.01006 * Residuals 8 30.820 3.853 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 # Note: alternative calls are possible: summary.aov(manova(aacon~gender)) summary.manova(aov(aacon~gender)) summary.aov(aov(aacon~gender)) summary.manova(manova(aacon~gender)) ################################################### # Multifactorial MANOVA # # Suppose we are given the following bird plasma calcium and water # loss data: # # No Hormone Treatment: # -------------------------------------------- # Female Male # -------------------- -------------------- # Plasma Ca H20 loss Plasma Ca H20 loss # --------- -------- --------- -------- # 16.5 76 14.5 80 # 18.4 71 11.0 72 # 12.7 64 10.8 77 # # Hormone Treatment: # -------------------------------------------- # Female Male # -------------------- -------------------- # Plasma Ca H20 loss Plasma Ca H20 loss # --------- -------- --------- -------- # 39.1 71 32.0 65 # 26.2 70 23.8 69 # 21.3 63 28.8 67 # # Read the data: data<-scan() 16.5 76 14.5 80 18.4 71 11.0 72 12.7 64 10.8 77 39.1 71 32.0 65 26.2 70 23.8 69 21.3 63 28.8 67 # # Put the data into a matrix: # <> "ncol=2" since there are two coordinates per measurement # <> "byrow=TRUE" to fill the matrix by rows # <> "dimnames" is optional but it puts nice labels atop the columns; # the first NULL means just default labels on the rows. # bird<-matrix(data,ncol=2,byrow=TRUE,dimnames=list(NULL,c("PLCA","H2OL"))) bird PLCA H2OL [1,] 16.5 76 [2,] 14.5 80 [3,] 18.4 71 [4,] 11.0 72 [5,] 12.7 64 [6,] 10.8 77 [7,] 39.1 71 [8,] 32.0 65 [9,] 26.2 70 [10,] 23.8 69 [11,] 21.3 63 [12,] 28.8 67 # # Make factor labels for the two factors: # <> the first 6 rows are for "no hormone treatment", the next 6 for # "hormone treatment". # <> rows are alternating "male" and "female". # <> there are 12 measurements in all. This is the default for # "hormone" but must be specified for "gender". # <> encapsulating "gl()" within "factor()" is optional. # hormone<-factor(gl(2,6,labels=c("No hormone treatment","Hormone treatment"))) gender<-factor(gl(2,1,12,labels=c("Female","Male"))) # # # Function "manova()" computes the sums of squares and crossed # products and labels them for use in ANOVAs and MANOVAs: # # All factors and interactions: # manova(bird~hormone*gender) Call: manova(bird ~ hormone * gender) Terms: hormone gender hormone:gender Residuals PLCA 635.1075 14.7408 7.2075 228.7533 H2OL 102.0833 18.7500 36.7500 151.3333 Deg. of Freedom 1 1 1 8 Residual standard error: 5.347351 4.349329 Estimated effects may be unbalanced # Get a p-value for the effect of gender on the combination of 3 # measurements with the "summary()" function. If applied to the # output of "manova()" it uses Pillai's statistic by default. # summary(manova(bird~hormone*gender)) Df Pillai approx F num Df den Df Pr(>F) hormone 1 0.8557 20.7624 2 7 0.001140 ** gender 1 0.2541 1.1922 2 7 0.358439 hormone:gender 1 0.3089 1.5647 2 7 0.274351 Residuals 8 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 # Use Wilks', Hotelling-Lawley's, or Roy's statistics instead: summary(manova(bird~hormone*gender),test="W") summary(manova(bird~hormone*gender),test="H") summary(manova(bird~hormone*gender),test="R") # Perform ANOVAs on the individual responses: summary(aov(bird~hormone*gender)) Response PLCA : Df Sum Sq Mean Sq F value Pr(>F) hormone 1 635.11 635.11 22.2111 0.001516 ** gender 1 14.74 14.74 0.5155 0.493191 hormone:gender 1 7.21 7.21 0.2521 0.629151 Residuals 8 228.75 28.59 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 Response H2OL : Df Sum Sq Mean Sq F value Pr(>F) hormone 1 102.083 102.083 5.3965 0.04869 * gender 1 18.750 18.750 0.9912 0.34861 hormone:gender 1 36.750 36.750 1.9427 0.20087 Residuals 8 151.333 18.917 --- Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 # Note: alternative calls are possible: summary.aov(manova(bird~hormone*gender)) summary.manova(aov(bird~hormone*gender)) summary.aov(aov(bird~hormone*gender)) summary.manova(manova(bird~hormone*gender))