tau.Rmd

---
title: "tau optimization"
author: "Lluís Revilla Sancho"
output:
  html_document:
    code_folding: show
    toc: true
    toc_float: true
date: "`r date()`"
---

```{r setup, include=FALSE, cache=FALSE}
knitr::opts_chunk$set(echo = TRUE, cache = TRUE, warning = FALSE, fig.width = 120, fig.height = 50, collapse = TRUE)
library("RGCCA")
library("gliomaData")
library("ggplot2")
library("patchwork")
theme_set(theme_bw())
```

# Introduction

Prepare the data

```{r prepare}
# Load data ####
data(ge_cgh_locIGR, package = "gliomaData")
A <- ge_cgh_locIGR$multiblocks
Loc <- factor(ge_cgh_locIGR$y)
levels(Loc) <- colnames(ge_cgh_locIGR$multiblocks$y)
ncomp <- rep(1, length(A))
# rgcca algorithm using the dual formulation for X1 and X2 
# and the dual formulation for X3
A[[3]] = A[[3]][, -3]
C <-  matrix(c(0, 0, 1, 0, 0, 1, 1, 1, 0), 3, 3)
dimnames(C) <- list(names(A), names(A))
C
```

# Preparing the taus

Testing the tau effect on the AVE

```{r taus}
(min_shrinkage <- sapply(A, function(x) {
  1 / sqrt(ncol(x))
}))
taus <- lapply(min_shrinkage, seq, to = 1, length.out = 10)
# ts <- lapply(A, tau.estimate) # On the computer
ts <- c(GE = 0.503293099215374, CGH = 0.636243186088126, y = 0.0398631441261977)
taus[["GE"]] <- c(taus[["GE"]], ts["GE"])
taus[["CGH"]] <- c(taus[["CGH"]], ts["CGH"])
taus.combn <- expand.grid(taus)
keep <- apply(taus.combn, 1, function(x){sum(x != 0)})
taus.combn <- taus.combn[keep >= 2, ]
# Add the example values at the bottom
taus.combn <- rbind(taus.combn, c(.071,.2, 1))
head(taus.combn)
```

# Calculate the AVE values

We iterate over the list of taus
```{r}
sgcca_eval <- function(x, scheme){
  result.sgcca <- RGCCA::sgcca(A, C, c1 = x, ncomp = ncomp, 
                        scheme = scheme, verbose = FALSE)
  
  ind <- apply(which(upper.tri(result.sgcca$C), arr.ind = TRUE), 1, 
               paste0, collapse = "")
  
  # Correlation between Y
  Y <- simplify2array(result.sgcca$Y, higher = FALSE)
  cY <- cor(Y)
  d <- cY * x
  var <- cY[upper.tri(cY)]
  names(var) <- paste0("vs", ind)
  # Canonical correlation
  cc <- sum(d[upper.tri(d)])
  
  c(unlist(result.sgcca$AVE[c("AVE_inner", "AVE_outer")]), "cc1" = cc)
}
```

```{r centroid}
out <- apply(taus.combn, 1, sgcca_eval, scheme = "centroid")
centroid <- cbind.data.frame(t(out), taus.combn)
saveRDS(centroid, "centroid_tau.RDS")
library("dplyr")
library("broom")
centroid %>%
  filter(AVE_inner == max(AVE_inner))
```

```{r}
# centroid <- cbind(centroid, scheme = "centroid")
# lmC0 <- lm(AVE_inner ~ 0+GE+CGH+y, data = centroid)
# lmC <- lm(AVE_inner ~ 0+GE*CGH*y, data = centroid)
# anova(lmC, lmC0)
# glance(lmC)
# tidy(lmC) %>% arrange(desc(abs(estimate)))
```

```{r factorial}
out <- apply(taus.combn, 1, sgcca_eval, scheme = "factorial")
factorial <- cbind.data.frame(t(out), taus.combn)
saveRDS(factorial, "factorial_tau.RDS")
factorial %>% 
  filter(AVE_inner == max(AVE_inner))
```


```{r}
factorial <- cbind(factorial, scheme = "factorial")
lmF0 <- lm(AVE_inner ~ 0+ GE+CGH+y, data = factorial)
lmF <- lm(AVE_inner ~ 0+ GE*CGH*y, data = factorial)
anova(lmF, lmF0)
glance(lmF)
tidy(lmF) %>% arrange(desc(abs(estimate)))
```

```{r horst}
out <- apply(taus.combn, 1, sgcca_eval, scheme = "horst")
horst <- cbind.data.frame(t(out), taus.combn)
saveRDS(horst, "horst_tau.RDS")
horst %>% 
  filter(AVE_inner == max(AVE_inner))
```

```{r}
horst <- cbind(horst, scheme = "horst")
lmH0 <- lm(AVE_inner ~ 0 +GE+CGH+y, data = horst)
lmH <- lm(AVE_inner ~ 0 +GE*CGH*y, data = horst)
anova(lmH, lmH0)
glance(lmH)
tidy(lmH) %>% arrange(desc(abs(estimate)))
```


```{r}
library("broom")
library("dplyr")
all_data <- rbind(horst, factorial, centroid)
lmALL0 <- lm(AVE_inner ~ 0+GE+CGH+y, data = all_data)
lmALL <- lm(AVE_inner ~ 0+GE*CGH*y, data = all_data)
glance(lmALL)
tidy(lmALL) %>% arrange(desc(abs(estimate)))
lmALL20 <- lm(AVE_inner ~ 0+scheme+GE+CGH+y, data = all_data)
lmALL2 <- lm(AVE_inner ~ 0+scheme*GE*CGH*y, data = all_data)
o <- anova(lmALL0, lmALL, lmALL2, lmALL20)
tidy(lmALL2) %>% arrange(desc(abs(estimate)))
glance(lmALL2)
tidy(lmALL2)
```

# Plot the result

```{r plot.centroid}
p <- ggplot(centroid)
p1 <- p + geom_point(aes(GE, AVE_inner, col = "inner")) + ylab("AVE") +
  geom_smooth(aes(GE, AVE_inner, col = "inner")) +
  geom_smooth(aes(GE, AVE_outer, col = "outer")) +
  geom_point(aes(GE, AVE_outer, col = "outer")) + guides(col = FALSE)
p2 <- p + geom_point(aes(CGH, AVE_inner, col = "inner")) + ylab("") +
  geom_smooth(aes(CGH, AVE_inner, col = "inner")) +
  geom_smooth(aes(CGH, AVE_outer, col = "outer")) +
  geom_point(aes(CGH, AVE_outer, col = "outer")) + guides(col = FALSE)
p3 <- p + geom_point(aes(y, AVE_inner, col = "inner")) + ylab("") +
  geom_smooth(aes(y, AVE_inner, col = "inner")) +
  geom_smooth(aes(y, AVE_outer, col = "outer")) +
  geom_point(aes(y, AVE_outer, col = "outer")) + labs(col = 'AVE type') 

q <- p + geom_count(aes(AVE_outer, AVE_inner)) + 
  xlab("AVE outer") + ylab("AVE inner")

(p1 + p2 + p3) / (q)
```

```{r plot.factorial}
p <- ggplot(factorial)
p1 <- p + geom_point(aes(GE, AVE_inner, col = "inner")) + ylab("AVE") +
  geom_smooth(aes(GE, AVE_inner, col = "inner")) +
  geom_smooth(aes(GE, AVE_outer, col = "outer")) +
  geom_point(aes(GE, AVE_outer, col = "outer")) + guides(col = FALSE)
p2 <- p + geom_point(aes(CGH, AVE_inner, col = "inner")) + ylab("") +
  geom_smooth(aes(CGH, AVE_inner, col = "inner")) +
  geom_smooth(aes(CGH, AVE_outer, col = "outer")) +
  geom_point(aes(CGH, AVE_outer, col = "outer")) + guides(col = FALSE)
p3 <- p + geom_point(aes(y, AVE_inner, col = "inner")) + ylab("") +
  geom_smooth(aes(y, AVE_inner, col = "inner")) +
  geom_smooth(aes(y, AVE_outer, col = "outer")) +
  geom_point(aes(y, AVE_outer, col = "outer")) + labs(col = 'AVE type') 

q <- p + geom_count(aes(AVE_outer, AVE_inner)) + 
  xlab("AVE outer") + ylab("AVE inner")

(p1 + p2 + p3) / (q)
```

```{r plot.horst}
p <- ggplot(horst)
p1 <- p + geom_point(aes(GE, AVE_inner, col = "inner")) + ylab("AVE") +
  geom_smooth(aes(GE, AVE_inner, col = "inner")) +
  geom_smooth(aes(GE, AVE_outer, col = "outer")) +
  geom_point(aes(GE, AVE_outer, col = "outer")) + guides(col = FALSE)
p2 <- p + geom_point(aes(CGH, AVE_inner, col = "inner")) + ylab("") +
  geom_smooth(aes(CGH, AVE_inner, col = "inner")) +
  geom_smooth(aes(CGH, AVE_outer, col = "outer")) +
  geom_point(aes(CGH, AVE_outer, col = "outer")) + guides(col = FALSE)
p3 <- p + geom_point(aes(y, AVE_inner, col = "inner")) + ylab("") +
  geom_smooth(aes(y, AVE_inner, col = "inner")) +
  geom_smooth(aes(y, AVE_outer, col = "outer")) +
  geom_point(aes(y, AVE_outer, col = "outer")) + labs(col = 'AVE type') 

q <- p + geom_count(aes(AVE_outer, AVE_inner)) + 
  xlab("AVE outer") + ylab("AVE inner")

(p1 + p2 + p3) / (q)
```


```{r}
out <- cbind(taus.combn, horst)
colnames(out) <- c("tau.GE", "tau.CGH", "tau.y", colnames(out)[4:9])
ggplot(out) + 
  geom_point(aes(tau.y, AVE_inner, col = tau.CGH)) + 
  scale_color_viridis_c() + 
  ylab("inner AVE") + 
  ggtitle("Colored by tau of CGH") + 
ggplot(out) + 
  geom_point(aes(tau.y, AVE_inner, col = tau.GE)) + 
  scale_color_viridis_c() + 
  ylab("inner AVE") + 
  ggtitle("Colored by tau of GE")

ggplot(out) + 
  geom_point(aes(tau.y, AVE_outer, col = tau.CGH)) + 
  scale_color_viridis_c() + 
  ylab("outer AVE") + 
  ggtitle("Colored by tau of CGH") + 
ggplot(out) + 
  geom_point(aes(tau.y, AVE_outer, col = tau.GE)) + 
  scale_color_viridis_c() + 
  ylab("outer AVE") + 
  ggtitle("Colored by tau of GE")

```

## Session Info {-}

```{r session, echo=FALSE}
devtools::session_info()
```