Finally all is clear _but_ one thing: the only “unknown” variable (to [...]]]>

Testing i.pca ~ prcomp(), m.eigensystem ~ princomp()

1. Briefly about PCA
2. The modules/functions that implement PCA in GRASS & R
3. My claims (Entitled Comments)
4. Evidence (=the numbers derived from i.pca, prcomp, princomp,
m.eigensystem using some MODIS surface reflectance bands).

Finally all is clear _but_ one thing: the only “unknown” variable (to
me) is still the Eigenvalue provided by i.pca. I can’t nail this one:
*** how is it calculated? *** Looks like it is some _weighted_
variance… !??

The only thing I noticed with some testing (not thoroughly though) is
that the Eigenvalue (=Variance) reported by “i.pca” is ~0.58 of the
variance reported by “prcomp()” (sdev^2) !? That is:

sqrt(Eigenvalue(i.pca)) / sdev(prcomp) = 0.58xxx

Hamish, if you have the time (or anyone , could you please “translate”
the code in “grass6_dev/imagery/i.pca/main.c”, concerning the “eigval”
variable, in some pseudocode. It is really the last thing to clarify ( I
think ).

Kindest regards, Nikos

————————
# 1. PCA ### ### ### ###

*. PCA is performed via:

– method 1 (eigenvectors) : by determining the eigenvalues and
eigenvectors of a given covariance or correlation matrix.

– covariance matrix -> non-standardised
– correlation matrix -> standardised

– method 2 (SVD): based on by a singular value decomposition of the
data matrix. This is a more general solution to PCA than the
“eigenvectors”.

…The difference between the two methods is explained in [1]…
…SVD is used for numerical accuracy… (R Documentation)
…Results between the two methos can differ a bit…

…Data centering reduces the square mean error of
*approximating* the input data
…Data scaling

*. PCA returns:

– eigenvalues (=variance _or_ sdev^2)
– eigenvectors (=loadings, rotation, weighting coefficient)

——————————————-
# 2. Implementations of PCA ### ### ### ###

* R’s implementation of:

– method 1 (eigen) is the _princomp()_ function

-*applies*:
\data centering
\data scaling

-returns:
\sdev  –  standard deviation of principal components (that is:
the sqrt(Variance)
\loadings  –  eigenvectors

– method 2 (SVD) is the _prcomp()_ function

-options:
\data centering (center= TRUE | FALSE) [default is TRUE]
\data scaling (scale=TRUE | FALSE) [default is FALSE]

-returns:
\sdev (as above)
\rotation (same as loadings)

* GRASS’ implementation of:

– method 1 (eigen) is the _m.eigensystem_ module

-returns:
E -> eigenvalue -> Variance(E) -> sdev(E)^2
V -> eigenvectors associated with E
N -> normalized eigenvectors V –> This is _data centering_.
W -> N vector multiplied by the square root of the magnitude of
the eigenvalue (E), in other words: W = N * sdev(E)

– method 2 (SVD) is the   _i.pca_ module

-returns:
\Eigenvalues (latest version by Hamish)
\Eigenvectors

—————————–
# 3. Comments ### ### ### ###

# 1: it is (?) meaningless to compare the numerical results
(eigenvectors) of “i.pca” with the results of “m.eigensystem” since they
implement two different algebraic solutions to PCA with different
options (data centering is NOT performed by i.pca). The different
algebraic solutions don’t give big differences. Data centering does.

# 2: “i.pca” corresponds to “prcomp( ,center=FALSE,scale=FALSE)” and
“m.eigensystem” corresponds to “princomp()”

# 3: The standard deviations (sdev) reported by prcomp(,center=TRUE,
scale=FALSE) _match_ the variances (=eigenvalues = sdev^2) reported by
“m.eigensystem”.

# 4: Obviously, “prcomp()” with center=TRUE and scaling=FALSE by default
gives (almost) same results as “princomp()”.

——————————————
# 4. Comparing PCA results ### ### ### ###

# prcomp() #############################################################

### Details for prcomp {stats} from R Documentation
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.
## more details in [2]
#Default settings of prcomp(): _center=TRUE_ and _scale=FALSE_

### with center=FALSE, scale=FALSE
prcomp(mod07, center=FALSE, scale=FALSE) <<== this corresponds to i.pca

Standard deviations:
[1] 4288.3788  476.8904  114.3971

Rotation:
PC1        PC2        PC3
MOD2007_242_500_sur_refl_b02 -0.6353238  0.7124070 -0.2980602
MOD2007_242_500_sur_refl_b06 -0.6485551 -0.2826985  0.7067234
MOD2007_242_500_sur_refl_b07 -0.4192135 -0.6423066 -0.6416403

### with center=TRUE, scale=FALSE
prcomp(mod07, center=TRUE, scale=FALSE)
Standard deviations:
[1] 857.5749 436.0928 108.5085

Rotation:
PC1         PC2        PC3
MOD2007_242_500_sur_refl_b02 0.4184468  0.83881578 -0.3482677
MOD2007_242_500_sur_refl_b06 0.7249935 -0.07751872  0.6843794
MOD2007_242_500_sur_refl_b07 0.5470710 -0.53886820 -0.6405735

### with center=TRUE, scale=TRUE
prcomp(mod07, center=TRUE, scale=TRUE)
Standard deviations:
[1] 1.5030740 0.8397807 0.1885121

Rotation:
PC1        PC2        PC3
MOD2007_242_500_sur_refl_b02 0.4807060  0.8202848 -0.3099267
MOD2007_242_500_sur_refl_b06 0.6561893 -0.1020540  0.7476634
MOD2007_242_500_sur_refl_b07 0.5816677 -0.5627768 -0.5873201

### with center=FALSE, scale=TRUE
prcomp(mod07, center=FALSE, scale=TRUE)
Standard deviations:
[1] 1.71889117 0.20787871 0.04689961

Rotation:
PC1         PC2        PC3
MOD2007_242_500_sur_refl_b02 -0.5747640  0.74015136 -0.3490305
MOD2007_242_500_sur_refl_b06 -0.5812896 -0.06907305  0.8107597
MOD2007_242_500_sur_refl_b07 -0.5759763 -0.66888331 -0.4699430

# i.pca ##############################################################

# perform pca
i.pca input=MOD07_b02,MOD07_b06,MOD07_b07 output=TEST

Eigen values, (vectors), and [percent importance]:
PC1 6307563.04 ( -0.64 -0.65 -0.42 ) [ 98.71% ]
PC2 78023.63 ( -0.71 0.28 0.64 ) [ 1.22% ]
PC3 4504.60 ( -0.30 0.71 -0.64 ) [ 0.07% ]

### Comment: Comparing with prcomp()’s results it is obvious that
“i.pca” implements the SVD method _without_ centering and _without_
scaling the data prior to the analysis. ###

# princomp() – cov #####################################################

# Details for princomp {stats} from R Documentation
The calculation is done using eigen on the correlation or covariance
matrix

### using _covariance_ matrix
# prints only standard deviations
princomp(mod07)

Call:
princomp(x = mod07)

Standard deviations:
Comp.1   Comp.2   Comp.3
857.5737 436.0922 108.5083

3  variables and  350596 observations.

# get loadings (=rotation = eigenvectors)
(princomp(mod07))$loadings

Loadings:
Comp.1 Comp.2 Comp.3
MOD2007_242_500_sur_refl_b02 -0.418  0.839  0.348
MOD2007_242_500_sur_refl_b06 -0.725        -0.684
MOD2007_242_500_sur_refl_b07 -0.547 -0.539  0.641

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

# m.eigensystem – cov #################################################

# Determine eigen values and vectors of _covariance_ matrix
## Marked with “#” are the _N_’s that correspond to the eigenvectors
derived from R’s princomp() ##
(echo 3; r.covar MOD07_b02,MOD07_b06,MOD07_b07)|m.eigensystem

r.covar: complete …
100%
E    778244.0258462029          .0000000000    79.20
V          .5006581842          .0000000000
V          .8256483300          .0000000000
V          .6155834548          .0000000000
N      #   .4372107421   #      .0000000000
N      #   .7210155161   #      .0000000000
N      #   .5375717557   #      .0000000000
W       385.6991853500          .0000000000
W       636.0664787886          .0000000000
W       474.2358050886          .0000000000
E    192494.5769628266          .0000000000    19.59
V         -.8689798010          .0000000000
V          .0996340298          .0000000000
V          .5731134848          .0000000000
N      #  -.8309940700   #      .0000000000
N      #   .0952787255   #      .0000000000
N      #   .5480609638   #      .0000000000
W      -364.5920328433          .0000000000
W        41.8027823088          .0000000000
W       240.4573848757          .0000000000
E     11876.4548199713          .0000000000     1.21
V          .2872248982          .0000000000
V         -.5731591248          .0000000000
V          .5351449518          .0000000000
N      #   .3439413070   #      .0000000000
N      #  -.6863370819   #      .0000000000
N      #   .6408165005   #      .0000000000
W        37.4824307850          .0000000000
W       -74.7964308085          .0000000000
W        69.8356366100          .0000000000

# princomp() – cov #####################################################

### using _correlation_ matrix
princomp(mod07, cor=TRUE)

Call:
princomp(x = mod07, cor = TRUE)

Standard deviations:
Comp.1    Comp.2    Comp.3
1.5030740 0.8397807 0.1885121

3  variables and  350596 observations.

# get loadings
(princomp(mod07, cor=TRUE))$loadings

Loadings:
Comp.1 Comp.2 Comp.3
MOD2007_242_500_sur_refl_b02 -0.481  0.820  0.310
MOD2007_242_500_sur_refl_b06 -0.656 -0.102 -0.748
MOD2007_242_500_sur_refl_b07 -0.582 -0.563  0.587

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

# m.eigensystem – corr #################################################

# Determine eigen values and vectors of _correlation_ matrix
## Marked with “#” are the _N_’s that correspond to the eigenvectors
derived from R’s princomp() ##
(echo 3; r.covar -r MOD07_b02,MOD07_b06,MOD07_b07)|m.eigensystem

r.covar: complete …
100%
E         2.2915877718          .0000000000    76.39
V         -.5755655569          .0000000000
V         -.7660355041          .0000000000
V         -.6809380186          .0000000000
N     #   -.4896413269   #      .0000000000
N     #   -.6516766616   #      .0000000000
N     #   -.5792830912   #      .0000000000
W         -.7412186091          .0000000000
W         -.9865075560          .0000000000
W         -.8769182329          .0000000000
E          .6740687010          .0000000000    22.47
V          .8667178982          .0000000000
V         -.1116525720          .0000000000
V         -.6069908335          .0000000000
N      #   .8145815825   #      .0000000000
N      #  -.1049362531   #      .0000000000
N      #  -.5704780699   #      .0000000000
W          .6687852213          .0000000000
W         -.0861544341          .0000000000
W         -.4683721194          .0000000000
E          .0343435272          .0000000000     1.14
V          .2486404469          .0000000000
V         -.6006166822          .0000000000
V          .4655120098          .0000000000
N      #   .3109794470   #      .0000000000
N      #  -.7512029762   #      .0000000000
N      #   .5822249325   #      .0000000000
W          .0576307320          .0000000000
W         -.1392129859          .0000000000
W          .1078979635          .0000000000

[1] http://www.snl.salk.edu/~shlens/pub/notes/pca.pdf

[2] Copy-pasted from R Documentation

### prcomp() returns the following:

1. sdev
the standard deviations of the principal components (i.e.,
the square roots of the eigenvalues of the cov./cor. matrix,
though the calculation is actually done with the singular
values of the data matrix).

2. rotation
the matrix of variable loadings (i.e., a matrix whose
columns contain the eigenvectors). The function princomp
returns this in the element loadings.

* Note

The signs of the columns of the rotation matrix are
arbitrary, and so may differ between different programs
for PCA, and even between different builds of R.

转载：http://osgeo-org.1560.n6.nabble.com/Testing-i-pca-prcomp-m-eigensystem-princomp-td4049804.html

PCA ############################### install.packages(“vegan”) library(vegan)

> STpcoa<-read.table(file=”bactera_16s_final.subsample.phylip.tre1.weighted.phylip.pcoa.axes”, header=T,row.names=1) > STpcoa axis1 axis2 axis3 axis4 Cellulose -0.020878 -0.234601 0.167454 0 Foodwaste -0.234592 0.221741 0.085802 0 Sludge 0.368882 0.100725 -0.010570 0 Xylan -0.113413 -0.087865 -0.242686 0 >pl.STpcoa<-princomp(STpcoa) > summary(pl.STpcoa) Importance of components: Comp.1 Comp.2 Comp.3 Comp.4 Standard deviation 0.2260563 0.1746944 0.1536319 0 Proportion of Variance 0.4856521 0.2900347 [...]]]> ###############################

PCA
###############################
install.packages(“vegan”)
library(vegan)

> STpcoa<-read.table(file=”bactera_16s_final.subsample.phylip.tre1.weighted.phylip.pcoa.axes”, header=T,row.names=1)
> STpcoa
axis1 axis2 axis3 axis4
Cellulose -0.020878 -0.234601 0.167454 0
Foodwaste -0.234592 0.221741 0.085802 0
Sludge 0.368882 0.100725 -0.010570 0
Xylan -0.113413 -0.087865 -0.242686 0
>pl.STpcoa<-princomp(STpcoa)
> summary(pl.STpcoa)
Importance of components:
Comp.1 Comp.2 Comp.3 Comp.4
Standard deviation 0.2260563 0.1746944 0.1536319 0
Proportion of Variance 0.4856521 0.2900347 0.2243133 0
Cumulative Proportion 0.4856521 0.7756867 1.0000000 1

> ls(pl.STpcoa)
[1] “call” “center” “loadings” “n.obs” “scale” “scores” “sdev”
> class(pl.STpcoa)
[1] “princomp”
> nmds.col<-c(rep(“green”, 1), rep(“blue”, 1), rep(“black”,1), rep(“red”,1))
> plot(pl.STpcoa$scores, col=nmds.col, pch=20)
> legend(x=0.12, y=0.25, c(“Cellulose”,”Foodwaste”,”Sludge”,”Xylan”),c(“green”,”blue”,”black”,”red”),bty=”n”)

> biplot(pl.STpcoa)

##############################
> pl2.STpcoa<-prcomp(STpcoa)
> class(pl2.STpcoa)
[1] “prcomp”
> pl2.STpcoa$sd^2
[1] 0.06813525 0.04069083 0.03147035 0.00000000
> summary(pl2.STpcoa)
Importance of components:
PC1 PC2 PC3 PC4
Standard deviation 0.2610 0.2017 0.1774 0
Proportion of Variance 0.4857 0.2900 0.2243 0
Cumulative Proportion 0.4857 0.7757 1.0000 1
> pl2.STpcoa$x
PC1 PC2 PC3 PC4
Cellulose -0.02087762 -0.2346017 0.16745306 0
Foodwaste -0.23459165 0.2217407 0.08580311 0
Sludge 0.36888225 0.1007250 -0.01056992 0
Xylan -0.11341297 -0.0878640 -0.24268626 0
> pl2.STpcoa$rotation
PC1 PC2 PC3 PC4
axis1 1.000000e+00 -9.352971e-08 -8.835157e-07 0
axis2 9.353330e-08 1.000000e+00 4.064575e-06 0
axis3 8.835153e-07 -4.064576e-06 1.000000e+00 0
axis4 0.000000e+00 0.000000e+00 0.000000e+00 1
> plot(pl2.STpcoa$x, col=nmds.col, pch=20)
> legend(x=0.12, y=0.25, c(“Cellulose”,”Foodwaste”,”Sludge”,”Xylan”),c(“green”,”blue”,”black”,”red”),bty=”n”)

> screeplot(pl2.STpcoa,type=”lines”,main=”Scree Plot”)

bactera_16s_final.subsample.phylip.tre1.weighted.phylip.pcoa.axes

‘princomp’ can only be used with more units than variables”

prcomp是基于SVD分解（svd()函数，princomp是基于特征向量eigen()函数)

Good video source:

http://www.youtube.com/watch?v=oZ2nfIPdvjY

http://www.youtube.com/watch?v=I5GxNzKLIoU&feature=relmfu

http://www.planta.cn/forum/viewtopic.php?t=16754&highlight=%D3%EF%D1%D4

###########################################

以下所有代码包括练习数据，都可在R平台上直接运行。

#主成分分析和主成分回归 主成分分析的思想是Pearson 1901年提出的，Hotelling 1933进一步发展 在R中，进行主成分分析用到princomp() 函数

用法 princomp(x, cor = FALSE, scores = TRUE, covmat = NULL, subset = rep(TRUE, nrow(as.matrix(x))), …)

# 分析用数据 # cor 是否用样本的协方差矩阵作主成分分析 prcomp() 二 summary()函数 三 [...]]]> Difference between prcomp and princomp:

‘princomp’ can only be used with more units than variables”

prcomp是基于SVD分解（svd()函数，princomp是基于特征向量eigen()函数)

Good video source:

http://www.youtube.com/watch?v=oZ2nfIPdvjY

http://www.youtube.com/watch?v=I5GxNzKLIoU&feature=relmfu

http://www.planta.cn/forum/viewtopic.php?t=16754&highlight=%D3%EF%D1%D4

###########################################

以下所有代码包括练习数据，都可在R平台上直接运行。

#主成分分析和主成分回归
主成分分析的思想是Pearson 1901年提出的，Hotelling 1933进一步发展
在R中，进行主成分分析用到princomp() 函数

用法
princomp(x, cor = FALSE, scores = TRUE, covmat = NULL,
subset = rep(TRUE, nrow(as.matrix(x))), …)

# 分析用数据
# cor 是否用样本的协方差矩阵作主成分分析
prcomp()
二 summary()函数
三 loadings()函数
四 predict() 函数
五 screeplot() 函数
六 biplot() 函数

实例
某中学随机抽取某年级30名学生，测量其身高，体重，胸围，坐高，针对这30名中学生身体四项指标数据做主成分分析。
student<-data.frame(
X1=c(148, 139, 160, 149, 159, 142, 153, 150, 151, 139,
140, 161, 158, 140, 137, 152, 149, 145, 160, 156,
151, 147, 157, 147, 157, 151, 144, 141, 139, 148 ),
X2=c(41, 34, 49, 36, 45, 31, 43, 43, 42, 31,
29, 47, 49, 33, 31, 35, 47, 35, 47, 44,
42, 38, 39, 30, 48, 36, 36, 30, 32, 38 ),
X3=c(72, 71, 77, 67, 80, 66, 76, 77, 77, 68,
64, 78, 78, 67, 66, 73, 82, 70, 74, 78,
73, 73, 68, 65, 80, 74, 68, 67, 68, 70),
X4=c(78, 76, 86, 79, 86, 76, 83, 79, 80, 74,
74, 84, 83, 77, 73, 79, 79, 77, 87, 85,
82, 78, 80, 75, 88, 80, 76, 76, 73, 78 )
)
#主成分分析
student.pr <- princomp(student, cor = TRUE)
#显示结果
summary(student.pr, loadings=TRUE)
#预测，显示各样本主成分的值
pre<-predict(student.pr)
#显示碎石图
screeplot(student.pr,type=”lines”)
# 主成分分析散点图
biplot(student.pr)

例二
对128个成年男子的身材进行测量，每人测得16项指标，身高，坐高，胸围，头高，裤长，下档，手长，领围，前胸，后背，肩厚，肩宽，袖长，肋围，腰围，腿肚，分别用X1-X16表示。16项指标的相关矩阵R。从相关矩阵出发进行主成分分析，随16项指标进行分类。
命令
x<-c(
1.00,
0.79, 1.00,
0.36, 0.31, 1.00,
0.96, 0.74, 0.38, 1.00,
0.89, 0.58, 0.31, 0.90, 1.00,
0.79, 0.58, 0.30, 0.78, 0.79, 1.00,
0.76, 0.55, 0.35, 0.75, 0.74, 0.73, 1.00,
0.26, 0.19, 0.58, 0.25, 0.25, 0.18, 0.24, 1.00,
0.21, 0.07, 0.28, 0.20, 0.18, 0.18, 0.29,-0.04, 1.00,
0.26, 0.16, 0.33, 0.22, 0.23, 0.23, 0.25, 0.49,-0.34, 1.00,
0.07, 0.21, 0.38, 0.08,-0.02, 0.00, 0.10, 0.44,-0.16, 0.23, 1.00,
0.52, 0.41, 0.35, 0.53, 0.48, 0.38, 0.44, 0.30,-0.05, 0.50, 0.24, 1.00,
0.77, 0.47, 0.41, 0.79, 0.79, 0.69, 0.67, 0.32, 0.23, 0.31, 0.10, 0.62, 1.00,
0.25, 0.17, 0.64, 0.27, 0.27, 0.14, 0.16, 0.51, 0.21, 0.15, 0.31, 0.17, 0.26, 1.00,
0.51, 0.35, 0.58, 0.57, 0.51, 0.26, 0.38, 0.51, 0.15, 0.29, 0.28, 0.41, 0.50, 0.63, 1.00,
0.21, 0.16, 0.51, 0.26, 0.23, 0.00, 0.12, 0.38, 0.18, 0.14, 0.31, 0.18, 0.24, 0.50, 0.65, 1.00
)
names<-c(“X1″, “X2″, “X3″, “X4″, “X5″, “X6″, “X7″, “X8″, “X9″,
“X10″, “X11″, “X12″, “X13″, “X14″, “X15″, “X16″)
R<-matrix(0, nrow=16, ncol=16, dimnames=list(names, names))
for (i in 1:16){
for (j in 1:i){
R<-x[(i-1)*i/2+j]; R[j,i]<-R
}
}
#主成分分析
pr<-princomp(covmat=R)
load<-loadings(pr)

#
plot(load[,1:2])
text(load[,1], load[,2], adj=c(-0.4, 0.3))

主成分回归
考虑进口总额Y与三个自变量：国内总产值，存储量，总消费量之间的关系。现收集了1949-1959共11年的数据，试做线性回归和主成分回归分析。
conomy<-data.frame(
x1=c(149.3, 161.2, 171.5, 175.5, 180.8, 190.7, 202.1, 212.4, 226.1, 231.9, 239.0),
x2=c(4.2, 4.1, 3.1, 3.1, 1.1, 2.2, 2.1, 5.6, 5.0, 5.1, 0.7),
x3=c(108.1, 114.8, 123.2, 126.9, 132.1, 137.7, 146.0, 154.1, 162.3, 164.3, 167.6),
y=c(15.9, 16.4, 19.0, 19.1, 18.8, 20.4, 22.7, 26.5, 28.1, 27.6, 26.3)
)

线性回归
lm.sol<-lm(y~x1+x2+x3, data=conomy)
summary(lm.sol)
主成分回归

# 主成分分析
conomy.pr<-princomp(~x1+x2+x3, data=conomy, cor=T)
summary(conomy.pr, loadings=TRUE)
pre<-predict(conomy.pr)
conomy$z1<-pre[,1]; conomy$z2<-pre[,2]
lm.sol<-lm(y~z1+z2, data=conomy)
summary(lm.sol)

> gtsdata_test=read.table(“gtsdata.txt”, header=T) > gtsenv=read.table(“gtsenv.txt”, header=T) > gtsdata_data_t<-t(gtsdata_data) > decorana(gtsdata_data_t)

Call: decorana(veg = gtsdata_data_t)

Detrended correspondence analysis with 26 segments. Rescaling of axes with 4 iterations.

DCA1 DCA2 DCA3 DCA4 Eigenvalues 0.8634 0.4834 0.23788 0 Decorana values 0.8721 0.3793 0.07223 0 Axis lengths 5.3292 2.1115 1.80907 0

> gts.ca=cca(gtsdata_data_t) > gts.ca Call: cca(X = [...]]]>  

> gtsdata_test=read.table(“gtsdata.txt”, header=T)
> gtsenv=read.table(“gtsenv.txt”, header=T)
> gtsdata_data_t<-t(gtsdata_data)
> decorana(gtsdata_data_t)

Call:
decorana(veg = gtsdata_data_t)

Detrended correspondence analysis with 26 segments.
Rescaling of axes with 4 iterations.

DCA1 DCA2 DCA3 DCA4
Eigenvalues 0.8634 0.4834 0.23788 0
Decorana values 0.8721 0.3793 0.07223 0
Axis lengths 5.3292 2.1115 1.80907 0

> gts.ca=cca(gtsdata_data_t)
> gts.ca
Call: cca(X = gtsdata_data_t)

Inertia Rank
Total 1.653
Unconstrained 1.653 3
Inertia is mean squared contingency coefficient

Eigenvalues for unconstrained axes:
CA1 CA2 CA3
0.8721 0.5037 0.2776

> plot(gts.ca,scaling=3)

> gtsdata_data_t_del<-gtsdata_data_t[1:3,]
> gtsdata_data_t_del

> gts.rda=rda(gtsdata_data_t_del,gtsenv)
> gts.rda
Call: rda(X = gtsdata_data_t_del, Y = gtsenv)

Inertia Proportion Rank
Total 101790 1
Constrained 101790 1 2
Unconstrained 0 0 0
Inertia is variance
Some constraints were aliased because they were collinear (redundant)

Eigenvalues for constrained axes:
RDA1 RDA2
81240 20549

plot(gts.rda,display=c(“sp”,”bp”,”si”),scaling=3)

gtsenv.txt

gtsdata.txt

a.str # t1 t2 t3 # g1 “1” “3” NA # g2 “2” “5” “6”

a.num <- apply(a.str, c(1,2), as.numeric)

a.num # t1 t2 t3 # g1 1 3 NA # g2 2 5 6

Note: 第一行，第一列位置要为空！！

gtsdata

gtsenv.txt

赖江山.pdf

>bac_4sampledata=read.csv(“/home/R_heatmap/4sample_R_cluster.csv”, sep=”\t”) > row.names(bac_4sampledata)<-bac_4sampledata$Group > bac_4sample_Datamatrix<-data.matrix(bac_4sampledata[,2:5]) > heatmap.2(bac_4sample_Datamatrix, distfun=dist,col=greenred(256), scale=”row”, key=TRUE, symkey=FALSE, density.info=”none”, trace=”none”, cexRow=0.5, cexCol=0.7,margin=c(7,30), keysize=1.5);

4sample_R_cluster_stdtop100

The reason your attempt at passing distfun = dist(method = 'euclidean') didn’t work is thatdistfun (and hclustfun) are supposed to simply be name of functions. So if you want [...]]]> Glancing at the code for heatmap.2 I’m fairly sure that the default is to use dist, and it’s default is in turn to use euclidean distances.

The reason your attempt at passing distfun = dist(method = 'euclidean') didn’t work is thatdistfun (and hclustfun) are supposed to simply be name of functions. So if you want to alter defaults and pass arguments you need to write a wrapper function like this:

heatmap.2(...,hclustfun = function(x) hclust(x,method = 'centroid'),...) 

As I mentioned, I’m fairly certain that heatmap.2 is using euclidean distances by default, but a similar solution can be used to alter the distance function used:

heatmap.2(...,distfun = function(x) dist(x,method = 'euclidean'),...)

library("gplots")
library("RColorBrewer")

test <- matrix(c(79,38.6,30.2,10.8,22,
81,37.7,28.4,9.7,19.9,
82,36.2,26.8,9.8,20.9,
74,29.9,17.2,6.1,13.9,
81,37.4,20.5,6.7,14.6),ncol=5,byrow=TRUE)
colnames(test) <- c("18:0","18:1","18:2","18:3","20:0")
rownames(test) <- c("Sample 1","Sample 2","Sample 3", "Sample 4","Sample 5")
test <- as.table(test)
mat=data.matrix(test)

heatmap.2(mat,
dendrogram="row",
Rowv=TRUE,
Colv=NULL,
distfun = dist,
hclustfun = hclust,
xlab = "Lipid Species",
ylab = NULL,
colsep=c(1),
sepcolor="black",
key=TRUE,
keysize=1,
trace="none",
density.info=c("none"),
margins=c(8, 12),
col=bluered
)