practice_3.Rmd

---
title: "R Notebook"
output: html_notebook
---
```{r}
#install_github("vqv/ggbiplot")
library(ggbiplot)
data(iris)
pca.obj <- prcomp(iris[,1:4], center=TRUE, scale.=TRUE)
str(pca.obj)
P <- ggbiplot(pca.obj,
         obs.scale = 1, 
         var.scale=1,
         ellipse=T,
         circle=F,
         varname.size=3,
         var.axes=T,
         groups=iris$Species, #no need for coloring, I'm making the points invisible
         alpha=0) #invisible points, I add them below
P$layers <- c(geom_point(aes(color=iris$Species), cex=5), P$layers) #add geom_point in a layer underneath (only way I have to change the size of the points in ggbiplot)
png(filename="test.png", height=600, width=600)
print(#or ggsave()
    P
)
dev.off()
```
```{r}
library(magrittr)
```

```{r}
#introduce NAs
iris$Sepal.Length[sample(1:150, 5)] <- NA
iris$Sepal.Width[sample(1:150, 5)] <- NA
iris$Petal.Length[sample(1:150, 5)] <- NA
iris$Petal.Width[sample(1:150, 5)] <- NA
#pca.obj2 <- prcomp(iris[,1:4], center=TRUE, scale.=TRUE) #cannot use prcomp with NAs
#source("https://bioconductor.org/biocLite.R")
#biocLite("pcaMethods")
library(pcaMethods)
pca.obj2 <- pcaMethods::pca(iris[,1:4], method="nipals", nPcs=3, center=TRUE, scale.=TRUE)
class(pca.obj2)  #class pcaRes pcaMethods
pca.obj2
str(pca.obj2)
ggbiplot(pca.obj2)
pca.obj2 %>% na.omit() %>%  ggbiplot()
```
```{r}
library(pcaMethods)
```
```{r}
data("metaboliteDataComplete")
mdc <- scale(metaboliteDataComplete, center=TRUE, scale=FALSE)
length(mdc)
class(mdc)
cond <- runif(length(mdc)) < 0.05
mdcOut <- mdc 
mdcOut[cond] <- 10
```
```{r}
resSvd <- pca(mdc, method = "svd", nPcs = 5, center = FALSE)
resSvdOut <- pca(mdcOut, method = "svd", nPcs = 5, center = FALSE)
resRobSvd <- pca(mdcOut, method = "robustPca", nPcs = 5, center = FALSE)
```

```{r}
par(mfrow=c(2,2))
plot(loadings(resSvd)[,1], loadings(resSvdOut)[,1], 
xlab = "Loading 1 SVD", ylab = "Loading 1 SVD with outliers")
plot(loadings(resSvd)[,1], loadings(resRobSvd)[,1],
xlab = "Loading 1 SVD", ylab = "Loading 1 robustSVD with outliers")
plot(loadings(resSvd)[,1], loadings(resPPCA)[,1],
xlab = "Loading 1 SVD", ylab = "Loading 1 PPCA with outliers = NA")
plot(loadings(resRobSvd)[,1], loadings(resPPCA)[,1],
xlab = "Loading 1 roubst SVD with outliers", ylab = "Loading 1 svdImpute with outliers = NA")

```