Post date: Nov 15, 2019 3:42:59 AM
## Look at the relationship between the B and S pairs: are they close to each other on the PCA?
data_4 <- cbind(size,pairs,data)
data.S <- data_4[data_4[,1]=="S",]
data.B <- data_4[data_4[,1]=="B",]
res.pca.S <- prcomp(data.S[,-c(1,2,3,4)], center = TRUE)
res.pca.B <- prcomp(data.B[,-c(1,2,3,4)], center = TRUE)
# basic pca plots
fviz_pca_ind(res.pca.S,
axes=c(1,2),
label = "none", # hide individual labels
habillage =data.S[,4], # color by groups
#palette = c("#af8dc3", "#E7B800", "#5ab4ac"),
point.size = 3,
title = "Big ramets"
#addEllipses = TRUE # Concentration ellipses
) +
theme(legend.position = "none")
# Now compare PCA1 for S and B
PCA1.S <- cbind(data.S[,1:4], res.pca.S$x[,1])
PCA1.B <- cbind(data.B[,1:4], res.pca.B$x[,1])
# Compare S to B and keep the indices of S that are found in B
ind <- is.element(PCA1.S[,2],PCA1.B[,2])
# sim = similar
PCA1.S.sim <- PCA1.S[ind,]
# Store unique values for S
x <- unique(PCA1.S.sim[,2])
# Compare B to the new S, and keep of B the values found in S
ind <- is.element(PCA1.B[,2],PCA1.S.sim[,2])
PCA1.B.sim <- PCA1.B[ind,]
# Store unique values for B
y <- unique(PCA1.B.sim[,2])
# Compare new B and new S
setequal(x,y)
# Calculate the mean of duplicates
# Small trees
keep <- c()
mean.S <- rep(NA,each = length(y)[1])
count <- 1
j <- 1
for (i in c(1:(dim(PCA1.S.sim)[1]-1)))
{
x1 <- PCA1.S.sim[i,5]
y1 <- PCA1.S.sim[i,2]
y2 <- PCA1.S.sim[i+1,2]
if (y1 == y2)
{
count = count + 1
}
if (y1 != y2 && count == 1)
{
mean.S[j] <- x1
count <- 1
j = j + 1
} else if (y1 != y2 && count > 1)
{
mean.S[j] <- mean(PCA1.S.sim[i-count:i,5])
count <- 1
j = j + 1
}
}
# Big trees
keep <- c()
mean.B <- rep(NA,each = length(y)[1])
count <- 1
j <- 1
for (i in c(1:(dim(PCA1.B.sim)[1]-1)))
{
x1 <- PCA1.B.sim[i,5]
y1 <- PCA1.B.sim[i,2]
y2 <- PCA1.B.sim[i+1,2]
if (y1 == y2)
{
count = count + 1
}
if (y1 != y2 && count == 1)
{
mean.B[j] <- x1
count <- 1
j = j + 1
} else if (y1 != y2 && count > 1)
{
mean.B[j] <- mean(PCA1.B.sim[i-count:i,5])
count <- 1
j = j + 1
}
}
plot(mean.S,mean.B)
see<-cbind(mean.S,mean.B)
I realize I should compare the PCA values after the PCA is calculated with all samples together.