This function provides a convenient interface for carrying out manova using the scores from PCA and the factors (groups) stored in a Spectra object. The function will do anova as well, if you only provide one vector of scores, though this is probably of limited use. A Spectra object contains group information stored in its spectra$groups element, but you can also use splitSpectraGroups to generate additional groups/factors that might be more useful than the original.

hypTestScores(spectra, pca, pcs = 1:3, fac = NULL, ...)

Arguments

spectra

An object of S3 class Spectra().

pca

An object of class prcomp.

pcs

An integer vector giving the PCA scores to use as the response in the manova analysis.

fac

A character vector giving the factors to be used in the manova. They will be searched for within the Spectra object.

...

Additional arguments to be passed downstream, in this case to aov. Untested.

Value

The results of the analysis print to the console unless assigned. If assigned, the object class is one of several described in aov depending upon the data passed to it.

Details

This function is an extraordinarily thin wrapper which helps the user to avoid writing a very tedious formula specification.

See also

splitSpectraGroups which can be used to create additional factor elements in the Spectra object, which can then be used with this function. Additional documentation at https://bryanhanson.github.io/ChemoSpec/

Author

Bryan A. Hanson (DePauw University).

Examples


data(metMUD2)

# Original factor encoding:
levels(metMUD2$groups)
#> [1] "AbcD" "AbCD" "ABcD" "ABCD"

# Split those original levels into 2 new ones (re-code them)
new.grps <- list(geneBb = c("B", "b"), geneCc = c("C", "c"))
mM3 <- splitSpectraGroups(metMUD2, new.grps)
#> 	Additional data was found: geneBb
#> 	Additional data was found: geneCc

# Now do the PCA and anova, with 3 ways to see the results
pca <- c_pcaSpectra(mM3)
#> 	Additional data was found: geneBb
#> 	Additional data was found: geneCc
res <- hypTestScores(mM3, pca, fac = c("geneBb", "geneCc"))
#> 	Additional data was found: geneBb
#> 	Additional data was found: geneCc
#> [1] "geneBb" "geneCc"
#> [1] "geneBb*geneCc"
res
#> Call:
#>    manova(formula = form)
#> 
#> Terms:
#>                   geneBb   geneCc geneBb:geneCc Residuals
#> PC1             476.0737   0.0814        2.1953    3.6544
#> PC2               0.0457  19.8592        3.6427    4.0036
#> PC3               0.0000   0.0749        0.1187    3.7585
#> Deg. of Freedom        1        1             1        16
#> 
#> Residual standard errors: 0.4779147 0.5002239 0.4846729
#> Estimated effects may be unbalanced
summary(res)
#>               Df  Pillai approx F num Df den Df    Pr(>F)    
#> geneBb         1 0.99270   634.87      3     14 3.450e-15 ***
#> geneCc         1 0.83976    24.46      3     14 7.902e-06 ***
#> geneBb:geneCc  1 0.57216     6.24      3     14  0.006517 ** 
#> Residuals     16                                             
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
summary.aov(res)
#>  Response PC1 :
#>               Df Sum Sq Mean Sq   F value    Pr(>F)    
#> geneBb         1 476.07  476.07 2084.3630 < 2.2e-16 ***
#> geneCc         1   0.08    0.08    0.3565  0.558820    
#> geneBb:geneCc  1   2.20    2.20    9.6117  0.006876 ** 
#> Residuals     16   3.65    0.23                        
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#>  Response PC2 :
#>               Df  Sum Sq Mean Sq F value    Pr(>F)    
#> geneBb         1  0.0457  0.0457  0.1826  0.674845    
#> geneCc         1 19.8592 19.8592 79.3659 1.337e-07 ***
#> geneBb:geneCc  1  3.6427  3.6427 14.5576  0.001522 ** 
#> Residuals     16  4.0036  0.2502                      
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#>  Response PC3 :
#>               Df Sum Sq  Mean Sq F value Pr(>F)
#> geneBb         1 0.0000 0.000001  0.0000 0.9981
#> geneCc         1 0.0749 0.074882  0.3188 0.5802
#> geneBb:geneCc  1 0.1187 0.118667  0.5052 0.4875
#> Residuals     16 3.7585 0.234908               
#> 

# You can also call this function on the existing groups:
res <- hypTestScores(metMUD2, pca, fac = "groups")
#> [1] "groups"
#> [1] "groups"