0

After conducting a PCA on a stack of rasters (similar to this & in the 2014 Raster Package documentation), I'd like to review my eigenvalues, eigenvectors, and loadings...

Typical calls for princomp which return the scree plot, proportion of variation, cumulative proportion - summary(), print(), hist() - don't seem to be pulling the information from my RasterBrick output. Here's example code: 

#bring in rasterBrick
logo <- brick(system.file("external/rlogo.grd", package="raster")) 
#select random samples & do PCA
sr <- sampleRandom(logo, 100) 
pca <- princomp(sr) 
# to visualize pcs as rasters
x <- predict(logo, pca, index=1:3)
plot(x)

##ANSWERED QUESTION
summary(pca) # importance of components
plot (pca) # scree plot
loadings (pca) #eigens

summary() returns what appears to be summary statistics on the raster layer not values from the analysis; print() shows min and max values of the raster, etc.

Thanks for your thoughts, especially any clarity on how to find the eignvalues & eigenvectors relevant to the PCA.

Community
  • 1
  • 1
lorena
  • 25
  • 1
  • 5

1 Answers1

2

The eigenvectors/loadings are stored in the loadings element of the model object returned by princomp. See the Value section of the help for princomp (run ?princomp). Here's the key section:

Value

princomp returns a list with class "princomp" containing the following components:

loadings the matrix of variable loadings (i.e., a matrix whose columns contain the eigenvectors). This is of class "loadings": see loadings for its print method.

You can access the loadings with loadings(pca). The first matrix below contains the eigenvector of each principal component.

loadings(pca)

Loadings:
      Comp.1 Comp.2 Comp.3
red    0.588 -0.505  0.631
green  0.584 -0.274 -0.764
blue   0.559  0.818  0.134

               Comp.1 Comp.2 Comp.3
SS loadings     1.000  1.000  1.000
Proportion Var  0.333  0.333  0.333
Cumulative Var  0.333  0.667  1.000

The summary function gives you the proportion of variance explained by each PC:

summary(pca)

Importance of components:
                            Comp.1      Comp.2       Comp.3
Standard deviation     136.9251939 16.85462507 1.4842831706
Proportion of Variance   0.9849601  0.01492417 0.0001157405
Cumulative Proportion    0.9849601  0.99988426 1.0000000000

Another thing you can always do with any R object is run str, which will tell you what the object contains. For example, see below for what a princomp model object contains, and note that one of the elements is loadings. 

str(pca)

List of 7
 $ sdev    : Named num [1:3] 136.5 17.63 1.43
  ..- attr(*, "names")= chr [1:3] "Comp.1" "Comp.2" "Comp.3"
 $ loadings: loadings [1:3, 1:3] 0.587 0.583 0.562 -0.515 -0.267 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:3] "red" "green" "blue"
  .. ..$ : chr [1:3] "Comp.1" "Comp.2" "Comp.3"
 $ center  : Named num [1:3] 162 165 173
  ..- attr(*, "names")= chr [1:3] "red" "green" "blue"
 $ scale   : Named num [1:3] 1 1 1
  ..- attr(*, "names")= chr [1:3] "red" "green" "blue"
 $ n.obs   : int 100
 $ scores  : num [1:100, 1:3] 85.4 110.4 151.3 149 22.8 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:3] "Comp.1" "Comp.2" "Comp.3"
 $ call    : language princomp(x = sr)
 - attr(*, "class")= chr "princomp"
 - attr(*, "class")= chr "princomp"
eipi10
  • 91,525
  • 24
  • 209
  • 285
  • thanks, though I'm running the pca on a rasterbrick, so many of these relationships seem to be breaking down. for example, a normal pca summary call returns standard deviation, proportion of variance, cumulative proportion. however, the rasterbrick pca hasn't given me this information... yet! here's a peak at example code `logo <- brick(system.file("external/rlogo.grd", package="raster")) names(logo) sr <- sampleRandom(logo, 100) pca <- princomp(sr) x <- predict(logo, pca, index=1:3)` – lorena Oct 15 '14 at 17:09
  • What do you mean by "...hasn't give me this information..."? Are you getting an error or are the output values just different from what you expected? I ran your code and it seems to work. It would be helpful if you updated your question to add the code you've run, the output from the code, and a description about what, specifically, is not what you wanted or expected. – eipi10 Oct 15 '14 at 17:31
  • Hopefully the updated question clarifies it. The code runs; however, the outputs aren't aligned with normal outputs for pca. It seems that I'm just getting information on the PC rasters. I might be misinterpreting the output, but that's why I'm asking the wider community for insight... – lorena Oct 15 '14 at 18:01
  • See the updates to my answer. Are you running `summary` on the raster or on the PCA model object? – eipi10 Oct 15 '14 at 18:13
  • I've been running `summary` and `print` on the raster object "x", which I now see was incorrect... Thanks for working through this with me, analyzing raster data in R has been a learning experience! – lorena Oct 15 '14 at 18:41
  • For future reference, please post your code, a sample of your data (if needed) and the output you expect in your initial question. This will make it *much* easier to help you. In this case, it would have been obvious immediately what was wrong if you had posted the exact code you actually ran. For more details, see this post on how to create a [reproducible example](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example). – eipi10 Oct 15 '14 at 18:44