[GRASS-user] i.pca vs. r.covar/m.eigensystem/r.mapcalc

GRASS GRASS user
1

I wanted to check/verify some i.pca results so I calculated 2 principal
components using r.covar, m.eigensystem and r.mapcalc. The results are
not the same with i.pca's output.

(Maybe this post http://www.nabble.com/Eigen-vectors-tt15686112.html was
about the same issue)

# i.pca performed on (3) modis surface reflectance bands (2, 6 and 7 -
500m spatial resolution)
r.info pca.2007.242.500.b267.1 -h
[...]
Eigen vectors:
   ( -0.64 -0.65 -0.42 )
   ( -0.71 0.29 0.64 )
   ( 0.29 -0.70 0.65 )
[...]

# calculating "by hand"
## (a) based on covariance
r.covar
map=MOD2007_242_500_sur_refl_b02,MOD2007_242_500_sur_refl_b06,MOD2007_242_500_sur_refl_b07 > covariance_b67

# eigenvalues and eigenvectors
(echo 3; cat covariance_b67) | m.eigensystem
E 193252.2225840657 0.0000000000 19.81
V 0.8480879794 0.0000000000
V -0.0973632655 0.0000000000
V -0.5621384304 0.0000000000
N 0.8297330080 0.0000000000
N -0.0952560549 0.0000000000
N -0.5499721988 0.0000000000
W 364.7544643092 0.0000000000
W -41.8750018841 0.0000000000
W -241.7703198653 0.0000000000

E 770286.8983235050 0.0000000000 78.97
V 0.5123041585 0.0000000000
V 0.8425333994 0.0000000000
V 0.6269758062 0.0000000000
N 0.4384249584 0.0000000000
N 0.7210319581 0.0000000000
N 0.5365598486 0.0000000000
W 384.7880047100 0.0000000000
W 632.8208355692 0.0000000000
W 470.9170625750 0.0000000000

E 11834.0908624288 0.0000000000 1.21
V 0.2892235465 0.0000000000
V -0.5746367180 0.0000000000
V 0.5358742690 0.0000000000
N 0.3454369570 0.0000000000
N -0.6863229556 0.0000000000
N 0.6400266474 0.0000000000
W 37.5782238381 0.0000000000
W -74.6613734450 0.0000000000
W 69.6250477264 0.0000000000

## (b) based on correlation
r.covar -r
map=MOD2007_242_500_sur_refl_b02,MOD2007_242_500_sur_refl_b06,MOD2007_242_500_sur_refl_b07 > correlation_b267

# eigenvalues and eigenvectors
(echo 3; cat correlation_b267) | m.eigensystem
E 2.2858261608 0.0000000000 76.19
V -0.5735822577 0.0000000000
V -0.7657536915 0.0000000000
V -0.6794534712 0.0000000000
N -0.4887914931 0.0000000000
N -0.6525548605 0.0000000000
N -0.5790121159 0.0000000000
W -0.7390013610 0.0000000000
W -0.9865943596 0.0000000000
W -0.8754054597 0.0000000000

E 0.6796466529 0.0000000000 22.65
V 0.8645087622 0.0000000000
V -0.1090761065 0.0000000000
V -0.6068722494 0.0000000000
N 0.8141383928 0.0000000000
N -0.1027208167 0.0000000000
N -0.5715130018 0.0000000000
W 0.6711812672 0.0000000000
W -0.0846837448 0.0000000000
W -0.4711592330 0.0000000000

E 0.0345271863 0.0000000000 1.15
V 0.2515640507 0.0000000000
V -0.6024904816 0.0000000000
V 0.4666495170 0.0000000000
N 0.3134669897 0.0000000000
N -0.7507466869 0.0000000000
N 0.5814790267 0.0000000000
W 0.0582468451 0.0000000000
W -0.1394999392 0.0000000000
W 0.1080474816 0.0000000000

# range for i.pca's pc 1
r.info pca.2007.242.500.b267.1 -r
min=-7713.335449
max=-21.046232

# range for i.pca's pc 2
r.info pca.2007.242.500.b267.2 -r
min=-2313.784484
max=1838.939406

# calculating PC 1 and 2 "by hand"
# pc 1
r.mapcalc testpca1="(-0.64)*MOD2007_242_500_sur_refl_b02 +
(-0.65)*MOD2007_242_500_sur_refl_b06 +
(-0.42)*MOD2007_242_500_sur_refl_b07"
# range of testpca1
min=-7743.680000
max=-87.370000

# pc 2
r.mapcalc testpca2="(-0.71)*MOD2007_242_500_sur_refl_b02 +
0.29*MOD2007_242_500_sur_refl_b06 + 0.64*MOD2007_242_500_sur_refl_b07"
# range of testpca2
min=-2283.720000
max=1866.570000

Can someone explain which values (V or N or W) correspond to the
eigenvectors reported by i.pca? I don't think they coincide. But, as far
as I understand the PCA concept, they should or they should be close
enough.

Kind regards, Nikos

2

My apologies for BUMPing :slight_smile:

On Sun, 2008-10-12 at 00:20 +0200, Nikos Alexandris wrote:

I wanted to check/verify some i.pca results so I calculated 2 principal
components using r.covar, m.eigensystem and r.mapcalc. The results are
not the same with i.pca's output.

(Maybe this post http://www.nabble.com/Eigen-vectors-tt15686112.html was
about the same issue)

# i.pca performed on (3) modis surface reflectance bands (2, 6 and 7 -
500m spatial resolution)

*** and region resolution set to 500m of course***

r.info pca.2007.242.500.b267.1 -h
[...]
Eigen vectors:
   ( -0.64 -0.65 -0.42 )
   ( -0.71 0.29 0.64 )
   ( 0.29 -0.70 0.65 )
[...]

I repeated the same i.pca but setting this time a different region
resolution (250m) than before (500m) and I get slightly different
eigenvectors:

*** with region resolution=250m ***
Eigen values:
( -0.64 -0.65 -0.42 )
( -0.71 0.28 0.64 )
( -0.30 0.71 -0.64 )

I understand that the region "resolution" influences the results of
i.pca. OK, then how would someone re-produce the same PCs doing the math
in R and not in GRASS?

Could someone explain or give any hints to learn more about i.pca?

Thank you, Nikos

3

On Monday 27 October 2008, Nikos Alexandris wrote:

My apologies for BUMPing :slight_smile:

On Sun, 2008-10-12 at 00:20 +0200, Nikos Alexandris wrote:
> I wanted to check/verify some i.pca results so I calculated 2 principal
> components using r.covar, m.eigensystem and r.mapcalc. The results are
> not the same with i.pca's output.
>
> (Maybe this post http://www.nabble.com/Eigen-vectors-tt15686112.html was
> about the same issue)
>
> # i.pca performed on (3) modis surface reflectance bands (2, 6 and 7 -
> 500m spatial resolution)

*** and region resolution set to 500m of course***

> r.info pca.2007.242.500.b267.1 -h
> [...]
> Eigen vectors:
> ( -0.64 -0.65 -0.42 )
> ( -0.71 0.29 0.64 )
> ( 0.29 -0.70 0.65 )
> [...]

I repeated the same i.pca but setting this time a different region
resolution (250m) than before (500m) and I get slightly different
eigenvectors:

*** with region resolution=250m ***
Eigen values:
( -0.64 -0.65 -0.42 )
( -0.71 0.28 0.64 )
( -0.30 0.71 -0.64 )

I understand that the region "resolution" influences the results of
i.pca. OK, then how would someone re-produce the same PCs doing the math
in R and not in GRASS?

Could someone explain or give any hints to learn more about i.pca?

Thank you, Nikos

Hi,

Remember that the region settings may cause re-sampling of your image due to
changes in resolution or grid-cell alignment. Here is a re-enactment of your
example above in GRASS, then in R- using some color imagery. Although I am
not an expert in PCA, I think that the methods in R/GRASS here are
comparable. In the output from R, the row and column names of the returned
matrix are printed, and it looks like the rotation vectors (eigen vectors)
for the 3rd PC are the least stable. This seems to make sense as (in this
case) the 3rd PC accounts for the least amount of variance-- and is thus
probably just picking up noise.

Not sure about applying the rotation manually and how that compares with
i.pca ...

Cheers,

Dylan

------------------- example ----------------------

# intial test at 1m res
g.region res=1 -a
i.pca in=naip.red,naip.green,naip.blue out=pca --o
d.rgb r=pca.2 g=pca.1 b=pca.3

Eigen values:
( -0.74 -0.55 -0.38 )
( 0.67 -0.67 -0.31 )
( 0.08 0.49 -0.87 )

# now try 10m res
g.region res=10 -a
i.pca in=naip.red,naip.green,naip.blue out=pca --o
d.rgb r=pca.2 g=pca.1 b=pca.3

Eigen values:
( -0.73 -0.56 -0.38 )
( 0.68 -0.65 -0.34 )
( 0.06 0.51 -0.86 )

## now try it in R

# init libraries
library(spgrass6)

# read in same data used in GRASS PCA calc
system('g.region res=1 -a')
x <- readRAST6(c('naip.red','naip.green', 'naip.blue'))

# use prcomp, as it directly returns rotation matrix
# same as Eigen values?
x.pca <- prcomp(x@data)

# transpose result for printing in same order as GRASS
signif(t(x.pca$rotation), 2)

    naip.red naip.green naip.blue
PC1 -0.730 -0.56 -0.38
PC2 0.680 -0.62 -0.39
PC3 0.023 0.54 -0.84

# do again, at reduces res:
system('g.region res=10 -a')

x <- readRAST6(c('naip.red','naip.green', 'naip.blue'))

# use prcomp, as it directly returns rotation matrix
# same as Eigen values- yes, see the manual page for prcomp
x.pca <- prcomp(x@data)

# transpose result for printing in same order as GRASS
signif(t(x.pca$rotation), 2)

PC1 -0.730 -0.56 -0.38
PC2 0.680 -0.64 -0.36
PC3 0.041 0.53 -0.85

--
Dylan Beaudette
Soil Resource Laboratory
http://casoilresource.lawr.ucdavis.edu/
University of California at Davis
530.754.7341

4

On Tue, 2008-10-28 at 09:39 -0700, Dylan Beaudette wrote:

On Monday 27 October 2008, Nikos Alexandris wrote:
> My apologies for BUMPing :slight_smile:
>
> On Sun, 2008-10-12 at 00:20 +0200, Nikos Alexandris wrote:
> > I wanted to check/verify some i.pca results so I calculated 2 principal
> > components using r.covar, m.eigensystem and r.mapcalc. The results are
> > not the same with i.pca's output.
> >
> > (Maybe this post http://www.nabble.com/Eigen-vectors-tt15686112.html was
> > about the same issue)
> >
> > # i.pca performed on (3) modis surface reflectance bands (2, 6 and 7 -
> > 500m spatial resolution)
>
> *** and region resolution set to 500m of course***
>
> > r.info pca.2007.242.500.b267.1 -h
> > [...]
> > Eigen vectors:
> > ( -0.64 -0.65 -0.42 )
> > ( -0.71 0.29 0.64 )
> > ( 0.29 -0.70 0.65 )
> > [...]
>
> I repeated the same i.pca but setting this time a different region
> resolution (250m) than before (500m) and I get slightly different
> eigenvectors:
>
> *** with region resolution=250m ***
> Eigen values:
> ( -0.64 -0.65 -0.42 )
> ( -0.71 0.28 0.64 )
> ( -0.30 0.71 -0.64 )
>
>
> I understand that the region "resolution" influences the results of
> i.pca. OK, then how would someone re-produce the same PCs doing the math
> in R and not in GRASS?
>
> Could someone explain or give any hints to learn more about i.pca?
>
> Thank you, Nikos
>

Hi,

Remember that the region settings may cause re-sampling of your image due to
changes in resolution or grid-cell alignment. Here is a re-enactment of your
example above in GRASS, then in R- using some color imagery. Although I am
not an expert in PCA, I think that the methods in R/GRASS here are
comparable. In the output from R, the row and column names of the returned
matrix are printed, and it looks like the rotation vectors (eigen vectors)
for the 3rd PC are the least stable. This seems to make sense as (in this
case) the 3rd PC accounts for the least amount of variance-- and is thus
probably just picking up noise.

Not sure about applying the rotation manually and how that compares with
i.pca ...

Cheers,

Dylan

[...]

Hi Dylan!

Sorry for replying so late and thank you for your detailed post. I now
see what I was doing wrong when trying to produce pc's with R. Here is
how I try to replicate your example using my data:

#load library(spgrass6) and data
library(spgrass6)
x <-
readRAST6(c('MOD2007_242_500_sur_refl_b02','MOD2007_242_500_sur_refl_b06','MOD2007_242_500_sur_refl_b07'))

# prcomp (gives me an error) ## maybe because of NULL's ??
x.pca <- prcomp(x@data)
Error in svd(x, nu = 0) : infinite or missing values in 'x'

----

# Trying with a subset with no NULL's ## resolution is set to 250m
i.pca in=b2,b6,b7 out=test_pca_b267

Eigen values:
( 0.79 0.61 0.06 )
( 0.44 -0.50 -0.74 )
( 0.42 -0.61 0.67 )

# set resolution to 500m and repeat i.pca
g.region res=500 -pa
i.pca in=b2,b6,b7 out=test_pca_b267_500m

Eigen values:
( 0.79 0.61 0.06 )
( 0.44 -0.50 -0.74 )
( 0.42 -0.61 0.67 )

## Well, this subset is rather small, thus no differences between 250m
and 500m (?)
g.region -p

[...]
cells: 7455
[...]

# trying with R
x <- readRAST6(c('b2','b6','b7'))
x.pca <- prcomp(x@data)

Error in svd(x, nu = 0) : infinite or missing values in 'x'

## again the same error!

I'll try another time again.

Kind regards, Nikos

5

On Wed, Nov 5, 2008 at 5:39 AM, Nikos Alexandris
<nikos.alexandris@felis.uni-freiburg.de> wrote:

On Tue, 2008-10-28 at 09:39 -0700, Dylan Beaudette wrote:

On Monday 27 October 2008, Nikos Alexandris wrote:
> My apologies for BUMPing :slight_smile:
>
> On Sun, 2008-10-12 at 00:20 +0200, Nikos Alexandris wrote:
> > I wanted to check/verify some i.pca results so I calculated 2 principal
> > components using r.covar, m.eigensystem and r.mapcalc. The results are
> > not the same with i.pca's output.
> >
> > (Maybe this post http://www.nabble.com/Eigen-vectors-tt15686112.html was
> > about the same issue)
> >
> > # i.pca performed on (3) modis surface reflectance bands (2, 6 and 7 -
> > 500m spatial resolution)
>
> *** and region resolution set to 500m of course***
>
> > r.info pca.2007.242.500.b267.1 -h
> > [...]
> > Eigen vectors:
> > ( -0.64 -0.65 -0.42 )
> > ( -0.71 0.29 0.64 )
> > ( 0.29 -0.70 0.65 )
> > [...]
>
> I repeated the same i.pca but setting this time a different region
> resolution (250m) than before (500m) and I get slightly different
> eigenvectors:
>
> *** with region resolution=250m ***
> Eigen values:
> ( -0.64 -0.65 -0.42 )
> ( -0.71 0.28 0.64 )
> ( -0.30 0.71 -0.64 )
>
>
> I understand that the region "resolution" influences the results of
> i.pca. OK, then how would someone re-produce the same PCs doing the math
> in R and not in GRASS?
>
> Could someone explain or give any hints to learn more about i.pca?
>
> Thank you, Nikos
>

Hi,

Remember that the region settings may cause re-sampling of your image due to
changes in resolution or grid-cell alignment. Here is a re-enactment of your
example above in GRASS, then in R- using some color imagery. Although I am
not an expert in PCA, I think that the methods in R/GRASS here are
comparable. In the output from R, the row and column names of the returned
matrix are printed, and it looks like the rotation vectors (eigen vectors)
for the 3rd PC are the least stable. This seems to make sense as (in this
case) the 3rd PC accounts for the least amount of variance-- and is thus
probably just picking up noise.

Not sure about applying the rotation manually and how that compares with
i.pca ...

Cheers,

Dylan

[...]

Hi Dylan!

Sorry for replying so late and thank you for your detailed post. I now
see what I was doing wrong when trying to produce pc's with R. Here is
how I try to replicate your example using my data:

#load library(spgrass6) and data
library(spgrass6)
x <-
readRAST6(c('MOD2007_242_500_sur_refl_b02','MOD2007_242_500_sur_refl_b06','MOD2007_242_500_sur_refl_b07'))

# prcomp (gives me an error) ## maybe because of NULL's ??
x.pca <- prcomp(x@data)
Error in svd(x, nu = 0) : infinite or missing values in 'x'

----

# Trying with a subset with no NULL's ## resolution is set to 250m
i.pca in=b2,b6,b7 out=test_pca_b267

Eigen values:
( 0.79 0.61 0.06 )
( 0.44 -0.50 -0.74 )
( 0.42 -0.61 0.67 )

# set resolution to 500m and repeat i.pca
g.region res=500 -pa
i.pca in=b2,b6,b7 out=test_pca_b267_500m

Eigen values:
( 0.79 0.61 0.06 )
( 0.44 -0.50 -0.74 )
( 0.42 -0.61 0.67 )

## Well, this subset is rather small, thus no differences between 250m
and 500m (?)
g.region -p

[...]
cells: 7455
[...]

# trying with R
x <- readRAST6(c('b2','b6','b7'))
x.pca <- prcomp(x@data)

Error in svd(x, nu = 0) : infinite or missing values in 'x'

## again the same error!

I'll try another time again.

Kind regards, Nikos

Sounds like NAs are getting in the way. This is a limitation of
several R functions, that becomes particularly annoying when working
with GIS data. Here are some tips on dealing with NA in GIS data:

http://casoilresource.lawr.ucdavis.edu/drupal/node/664

Cheers,

Dylan

6

On Wed, 2008-11-05 at 07:06 -0800, Dylan Beaudette wrote:

> # trying with R
> x <- readRAST6(c('b2','b6','b7'))
> x.pca <- prcomp(x@data)
>
> Error in svd(x, nu = 0) : infinite or missing values in 'x'
>
>
> ## again the same error!
>
> I'll try another time again.
>
> Kind regards, Nikos
>

Sounds like NAs are getting in the way. This is a limitation of
several R functions, that becomes particularly annoying when working
with GIS data. Here are some tips on dealing with NA in GIS data:

http://casoilresource.lawr.ucdavis.edu/drupal/node/664

Cheers,

Dylan

Dylan, great info!!!

I've bombed in another wall :frowning:

x.nas <- which(is.na(x@data@MOD2007_242_500_sur_refl_b02) &

is.na(x@data@MOD2007_242_500_sur_refl_b06) & is.na(x@data
$MOD2007_242_500_sur_refl_b07))
Error in which(is.na(x@data@MOD2007_242_500_sur_refl_b02) &
is.na(x@data@MOD2007_242_500_sur_refl_b06) & :
  trying to get slot "MOD2007_242_500_sur_refl_b02" from an object
(class "data.frame") that is not an S4 object

I am working on it. Maybe something relevant with this error in [1]

[1] https://stat.ethz.ch/pipermail/r-devel/2007-December/047743.html

7

On Friday 07 November 2008, Nikos Alexandris wrote:

On Wed, 2008-11-05 at 07:06 -0800, Dylan Beaudette wrote:
> > # trying with R
> > x <- readRAST6(c('b2','b6','b7'))
> > x.pca <- prcomp(x@data)
> >
> > Error in svd(x, nu = 0) : infinite or missing values in 'x'
> >
> >
> > ## again the same error!
> >
> > I'll try another time again.
> >
> > Kind regards, Nikos
>
> Sounds like NAs are getting in the way. This is a limitation of
> several R functions, that becomes particularly annoying when working
> with GIS data. Here are some tips on dealing with NA in GIS data:
>
> http://casoilresource.lawr.ucdavis.edu/drupal/node/664
>
> Cheers,
>
> Dylan

Dylan, great info!!!

I've bombed in another wall :frowning:

> x.nas <- which(is.na(x@data@MOD2007_242_500_sur_refl_b02) &

is.na(x@data@MOD2007_242_500_sur_refl_b06) & is.na(x@data
$MOD2007_242_500_sur_refl_b07))
Error in which(is.na(x@data@MOD2007_242_500_sur_refl_b02) &
is.na(x@data@MOD2007_242_500_sur_refl_b06) & :
  trying to get slot "MOD2007_242_500_sur_refl_b02" from an object
(class "data.frame") that is not an S4 object

I am working on it. Maybe something relevant with this error in [1]

[1] https://stat.ethz.ch/pipermail/r-devel/2007-December/047743.html

You have a syntax error here:
x@data@MOD2007_242_500_sur_refl_b06
          ^^^^

should be:
x@data$MOD2007_242_500_sur_refl_b06

subtle but important.

$ is used to access columns from a dataframe
@ is used to access slots

Cheers,

Dylan

--
Dylan Beaudette
Soil Resource Laboratory
http://casoilresource.lawr.ucdavis.edu/
University of California at Davis
530.754.7341

8

On Fri, 2008-11-07 at 09:54 -0800, Dylan Beaudette wrote:

On Friday 07 November 2008, Nikos Alexandris wrote:
> On Wed, 2008-11-05 at 07:06 -0800, Dylan Beaudette wrote:

[...]

> I've bombed in another wall :frowning:
>
> > x.nas <- which(is.na(x@data@MOD2007_242_500_sur_refl_b02) &
>
> is.na(x@data@MOD2007_242_500_sur_refl_b06) & is.na(x@data
> $MOD2007_242_500_sur_refl_b07))
> Error in which(is.na(x@data@MOD2007_242_500_sur_refl_b02) &
> is.na(x@data@MOD2007_242_500_sur_refl_b06) & :
> trying to get slot "MOD2007_242_500_sur_refl_b02" from an object
> (class "data.frame") that is not an S4 object

[...]

You have a syntax error here:
x@data@MOD2007_242_500_sur_refl_b06
          ^^^^

should be:
x@data$MOD2007_242_500_sur_refl_b06

subtle but important.

$ is used to access columns from a dataframe
@ is used to access slots

Cheers,
Dylan

Dylan, thank you. It works now... sort of!? Sorry for another long post.
Perhaps I should continue this thread in another ML.

#read raster maps
x <-
readRAST6(c('MOD2007_242_500_sur_refl_b02','MOD2007_242_500_sur_refl_b06','MOD2007_242_500_sur_refl_b07'))

#extract NAs
x.nas <- which(is.na(x@data$MOD2007_242_500_sur_refl_b02) & is.na(x@data
$MOD2007_242_500_sur_refl_b06) & is.na(x@data
$MOD2007_242_500_sur_refl_b07))

#vectors pointing to no-NAs
x.vals <- which( !is.na(x@data$MOD2007_242_500_sur_refl_b02) & !
is.na(x@data$MOD2007_242_500_sur_refl_b06) & !is.na(x@data
$MOD2007_242_500_sur_refl_b07))

#get no-NAs
x.no.na <- x@data[x.vals, ]

#pca
x.pca <- prcomp(x.no.na)

# output
x.pca

Standard deviations:
[1] 881.8701 438.7193 108.9755

Rotation:
                                   PC1 PC2 PC3
MOD2007_242_500_sur_refl_b02 0.4371624 0.83100302 -0.3439811
MOD2007_242_500_sur_refl_b06 0.7210191 -0.09520195 0.6863440
MOD2007_242_500_sur_refl_b07 0.5376063 -0.54806074 -0.6407877

# probably the "correct" function to use (and compare) is the princomp
and not prcomp
## read [1] [2] for details

# princomp on x.no.na ## now I am stuck here :slight_smile:
x.pca_corr <- princomp(x.no.na, cor = TRUE)

# x.pca_corr
Call:
princomp(x = x.no.na, cor = TRUE)

Standard deviations:
   Comp.1 Comp.2 Comp.3
1.5137012 0.8211861 0.1853699

3 variables and 88101 observations.

# or "cor = FALSE"
x.pca_cova <- princomp(x.no.na, cor = FALSE)

# x.pca_cova
Call:
princomp(x = x.no.na, cor = FALSE)

Standard deviations:
  Comp.1 Comp.2 Comp.3
881.8651 438.7169 108.9748

3 variables and 88101 observations.

# I would like to compare R's results with i.pca's but I don't know why
or how I can control the "print-out" of the eigenvalues

## If I get the eigenvalues I would...
#compare with i.pca reported eigenvalues## read first mails on this
thread :slight_smile:

Eigen values:
( -0.64 -0.65 -0.42 )
( -0.71 0.28 0.64 )
( -0.30 0.71 -0.64 )

My understanding is that i.pca reports in rows, so in mt exampl it is:
pc1 ( -0.64 -0.65 -0.42 )
pc2 ( -0.71 0.28 0.64 )
pc3 ( -0.30 0.71 -0.64 )

Some interesting stuff to check:

1. prcomp and princomp accept an argument "na.action"

a function which indicates what should happen when the data contain
'NA's.
The default is set by the 'na.action' setting of 'options', and is
'na.fail' if that
is unset. The 'factory-fresh' default is 'na.omit'.

## question: (something like na.action = na.omit would help?)

---

[1] taken from "?prcomp":
The calculation is done by a singular value decomposition of the
(centered and possibly scaled) data matrix, not by using 'eigen'
on the covariance matrix. This is generally the preferred method
for numerical accuracy. The 'print' method for these objects
prints the results in a nice format and the 'plot' method produces
a scree plot.

[2] "?princomp":
The calculation is done using 'eigen' on the correlation or
covariance matrix, as determined by 'cor'. This is done for
compatibility with the S-PLUS result. A preferred method of
calculation is to use 'svd' on 'x', as is done in 'prcomp'.