run-tvmodels.Rmd

---
title: "Time-variant effects"
author: "Julius J."
date: "`r Sys.Date()`"
params:
  mfrfile: "/mnt/HARVEST/ga_cleaned.csv"
  gtfile: "/mnt/HARVEST/top1-moba30k-dosage.csv.gz"
  mobaresfile: "/mnt/HARVEST/topsnps_moba_summaries.txt"
  i: 1
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(dplyr)
library(ggplot2)
library(survival)
library(mgcv)
library(pammtools)
library(cowplot)
library(knitr)
library(kableExtra)
library(tidyr)
library(broom)

knitr::opts_chunk$set(fig.width=7, fig.height=4)
snpnum = as.numeric(params$i)

# unified palette for plotting
GTpalette =  c("#A6D96A", "#FD8D3C", "#D7191C")
```

# Report on SNP `r snpnum`

## Data

```{r genodata}
# read in filtered phenotypes
mfr_mid = read.table(params$mfrfile, h=T, sep=";")
nrow(mfr_mid)  # 25648

# read in genotypes
gt = data.table::fread(params$gtfile, sep=" ", h=T)
gtinfo = gt[,1:6]
gt = t(gt[,7:ncol(gt)])
gt = data.frame(gt)
gt$SENTRIX_ID = rownames(gt)
rownames(gt) = NULL
gt = filter(gt, !is.na(X1)) # tends to add an empty line in the end

# merge w/ pheno data
merged = inner_join(gt, mfr_mid, by=c("SENTRIX_ID"))
nrow(merged)  # all 25648 for autosomes

# remove irrelevant timeframe for stability
merged$GAc = merged$SVLEN_DG-169

# just some checks
sum(duplicated(merged$SENTRIX_ID))  # 0 
sum(duplicated(merged$PREG_ID_1724))  # 0
range(merged$GAc)  # 1 125

# Prepare covariates:
merged$PARITET_5[merged$PARITET_5==4] = 3  # parity honestly doesn't need 5 levels
merged$PARITET_5 = factor(merged$PARITET_5, levels=c(1,0,2,3))
# NOTE: using height will remove extra 1k missing values, so imputing:
merged$AA87[is.na(merged$AA87)] = mean(merged$AA87, na.rm=T)
# relevel to ensure the oldest and one of the largest batches is ref
merged$BATCH = factor(merged$BATCH, levels=c("M12A","FEB18","JAN15","JUN15",
                                             "M12B","M24","MAY16",
                                             "ROT1","ROT2"))
```

# Analysis

## Select the marker and assign alleles

```{r definesnp, echo=F}
# select the right marker
merged$GT = merged[,paste0("X",snpnum)]
gtinfo$REF.orig = gtinfo$REF
gtinfo$ALT.orig = gtinfo$ALT

# NOTE: flipping alleles to effect=minor alignment,
# because the last analysis is sensitive to that.
if(mean(merged$GT)>1){
  merged$GT = 2-merged$GT
  print("NOTE: allele flipped")
  # flip the gtinfo row accordingly
  gtinfo[snpnum,"REF"]  = gtinfo[snpnum,"ALT.orig"]
  gtinfo[snpnum,"ALT"]  = gtinfo[snpnum,"REF.orig"]
}

merged$GTcat = factor(round(merged$GT))

# print some info
gttable = table(merged$GT, dnn=c("dosage"))
kable(split(gttable, (1:length(gttable))%/%8))
gtinfo[snpnum,] %>% 
  separate(INFO, c(NA, "INFO", "RefPanelAF", "AN", "AC"), sep="[;A-Za-z]*=") %>%
  kable

# merge in some info from the MOBA results table
metares = read.table(params$mobaresfile, h=T)
metares$CHR = as.character(metares$CHR)
gtinfo$CHROM = as.character(gtinfo$CHROM)

select(gtinfo[snpnum,], -one_of(c("REF.orig", "ALT.orig", "INFO", "ALT"))) %>%
  left_join(metares, by=c("CHROM"="CHR", "POS")) %>%
  kable
```

## Some sanity checks

```{r checksnp, echo=F}
ggplot(merged[merged$hadevent,], aes(x=SVLEN_DG, group=GTcat)) +
  geom_vline(xintercept=259, col="grey60") +
  geom_density(aes(col=GTcat, fill=GTcat), alpha=0.05) +
  scale_fill_manual(values=GTpalette) + scale_color_manual(values=GTpalette) +
  theme_bw()
merged %>%
  group_by(GTcat) %>%
  summarize(n(), sum(SVLEN_DG<259), sum(SVLEN_DG<240), mean(SVLEN_DG<259)) %>%
  kable
```

## Basic models

### The most appropriate model:

```{r basicmodels, echo=F}
# to confirm effect size & direction
summary(lm(SVLEN_DG ~ GT + BATCH + poly(MAGE,2) + FAAR + AA87 + KJONN + MISD + PARITET_5,
           data=merged[merged$hadevent,]))
```

### Overview of slightly different models, w/o covariates

```{r basicmodels2, echo=F}
sum1 = bind_rows(
  "lm"=tidy(lm(SVLEN_DG ~ GT, data=merged[merged$hadevent,])),
  "lm.cat"=tidy(lm(SVLEN_DG ~ GTcat, data=merged[merged$hadevent,])),
  "b"=tidy(glm(SVLEN_DG<259 ~ GT, family="binomial", data=merged[merged$hadevent,])),
  "b.cat"=tidy(glm(SVLEN_DG<259 ~ GTcat, family="binomial", data=merged[merged$hadevent,])),
  .id="model")
sum1$batch=" "
sum2 = bind_rows(
  "lm"=tidy(lm(SVLEN_DG ~ GT + BATCH, data=merged[merged$hadevent,])),
  "lm.cat"=tidy(lm(SVLEN_DG ~ GTcat + BATCH, data=merged[merged$hadevent,])),
  "b"=tidy(glm(SVLEN_DG<259 ~ GT + BATCH, family="binomial", data=merged[merged$hadevent,])),
  "b.cat"=tidy(glm(SVLEN_DG<259 ~ GTcat + BATCH, family="binomial", data=merged[merged$hadevent,])),
  .id="model")
sum2$batch="+batch"
bind_rows(sum1, sum2)[,c("model", "batch", "term", "estimate", "std.error", "p.value")] %>%
  filter(term!="(Intercept)") %>%
  mutate(sig=ifelse(p.value<0.05, "*", "")) %>%
  kable
```

## Cox PH model

```{r cox, echo=F}
summary(coxph(Surv(GAc, hadevent) ~ GT, data=merged))
mod.ph = coxph(Surv(GAc, hadevent) ~ GT + BATCH + poly(MAGE,2) + FAAR +
                 AA87 + KJONN + MISD + PARITET_5, data=merged)
summary(mod.ph)

par(mfrow=c(1, 2))
km.ph = survfit(coxph(Surv(GAc, hadevent) ~ GTcat + BATCH + poly(MAGE,2) +
                        FAAR + AA87 + KJONN + MISD + PARITET_5, data=merged),
                newdata=data.frame(GTcat=factor(0:2), BATCH="M12A", MAGE=21, FAAR=2005,
                                   AA87=170, KJONN=1, MISD=F, PARITET_5=factor(1)))
plot(km.ph, conf.int=F, col=4:2, ylab="S(t)")
title(main="K-M curve for coxph(~GTcat + covs)")
legend(5, 0.4, levels(merged$GTcat), col=4:2, lty=1)
plot(km.ph, conf.int=F, col=4:2, cumhaz=T, log=T, ylab="Lambda(t)")
legend(5, 0.4, levels(merged$GTcat), col=4:2, lty=1)
```

## TV by stratified K-M

```{r tv2, echo=F}
# note that we're ignoring covariates here.
# If this is to be present for visuals, need to find a way to incorporate them
par(mfrow=c(1, 2))
km.tv = survfit(Surv(GAc, hadevent) ~ GTcat, data=merged)
plot(km.tv, col=4:2, ylab="S(t)")
title(main="K-M curve, strata by GT, no covs")
legend(5, 0.4, levels(merged$GTcat), col=4:2, lty=1)
plot(km.tv, conf.int=F, col=4:2, cumhaz=T, log=T, ylab="Lambda(t)")
legend(5, 0.4, levels(merged$GTcat), col=4:2, lty=1)
```

## TV by testing Schoenfeld residuals

```{r tvtests}
# Default transform supposedly downweighs outliers.
# Could use the identity one, but that is identical to the interaction w/ t below.
# zp.ph = cox.zph(mod.ph, transform='identity')
zp.ph = cox.zph(mod.ph)
zp.ph
plot(zp.ph, resid=F)
```

## TV by Cox w/ prespecified shapes

```{r coxint}
# via specified forms in survival
mod.tv.1 = coxph(Surv(GAc, hadevent) ~ GT + tt(GT) + BATCH + poly(MAGE,2) +
                   FAAR + AA87 + KJONN + MISD + PARITET_5, data=merged,
                 tt=function(x, t, ...) x*(t/100))
summary(mod.tv.1)

# contrast for PTD
mod.tv.2 = coxph(Surv(GAc, hadevent) ~ GT + tt(GT) + BATCH + poly(MAGE,2) +
                   FAAR + AA87 + KJONN + MISD + PARITET_5, data=merged,
                 tt=function(x, t, ...) x*(t<89))
summary(mod.tv.2)
```


## Piecewise additive models for PH

```{r pammph, fig.height=3}
# NOTE: using cubic regr splines for consistency betw. ti() and s()
# convert to piecewise format - i.e. 1 row per each period per indiv.:
ped = as_ped(merged, Surv(GAc, hadevent)~GT + GTcat + BATCH +
               MAGE + FAAR + AA87 + KJONN + MISD + PARITET_5,
             id = "id", cut=c(0,seq(20, 130, by=7)))
nrow(ped)  # 367824
ped$MISD = factor(ped$MISD) # to prevent gg_slice complaining
ped$KJONN = factor(ped$KJONN)
ped$FAAR = ped$FAAR-2000  # center, for neater plot scales

# fit a piecewise PH model for testing

mod.pc = bam(ped_status ~ s(tend,bs='cr',k=11) + GT + BATCH +
               poly(MAGE,2) + FAAR + AA87 + KJONN + MISD + PARITET_5,
             data=ped, offset=offset, family=poisson())
summary(mod.pc)
plot(mod.pc, pages=2, all.terms=T, scale=0)

# confirm that doesn't differ much from Cox
cbind(coef(mod.pc)[-1], c(coef(mod.ph), rep(NA, 10)))

```

## PAM w/ prespecified shape

```{r pammtvshape, fig.height=3}
# time-varying by prespecified shape
mod.tv.pc.i = bam(ped_status ~ s(tend,bs='cr',k=11) + GT + GT:tend + BATCH +
               poly(MAGE,2) + FAAR + AA87 + KJONN + MISD + PARITET_5,
               data=ped, offset=offset, family=poisson())
summary(mod.tv.pc.i)
plot(mod.tv.pc.i, pages=5, all.terms=T, scale=0)
```

## PAM w/ smooth interactions

```{r pammsmooth, fig.width=7, fig.height=6}
# time-varying smoothly (interaction with spline)
# explanation: ti(tend) is the baseline h(x), and ti(tend, by=GT)
# is a h_1(x)*GT term with potentially different shape.
mod.tv.pc.sm = bam(ped_status ~ ti(tend,bs='cr',k=11) + GT + ti(tend, by=GT,bs='cr') +
                     BATCH + poly(MAGE,2) + FAAR + AA87 + KJONN + MISD + PARITET_5,
                     data=ped, offset=offset, family=poisson())
summary(mod.tv.pc.sm)
plot(mod.tv.pc.sm, pages=3, all.terms=T, scale=0)
# NOTE: bam is used for speed, but be aware that it differs from gam sometimes.

# NOTE: this one is sensitive to allele choice, so should try both ways
# non-linear in time AND non-linear in GT
mod.tv.pc.sm2 = bam(ped_status ~ ti(tend,bs='cr',k=11) + GTcat + ti(tend, by=as.ordered(GTcat),bs='cr') +
                      BATCH + poly(MAGE,2) + FAAR + AA87 + KJONN + MISD + PARITET_5,
                      data=ped, offset=offset, family=poisson())
summary(mod.tv.pc.sm2)
plot(mod.tv.pc.sm2, pages=3, all.terms=T, scale=0)

ped$GTcatRev = factor(round(2-ped$GT))
mod.tv.pc.sm2.r = bam(ped_status ~ ti(tend,bs='cr',k=11) + GTcatRev + ti(tend, by=as.ordered(GTcatRev),bs='cr') +
                        BATCH + poly(MAGE,2) + FAAR + AA87 + KJONN + MISD + PARITET_5,
                        data=ped, offset=offset, family=poisson())
summary(mod.tv.pc.sm2.r)
plot(mod.tv.pc.sm2.r, pages=3, all.terms=T, scale=0)

AIC(mod.pc)
AIC(mod.tv.pc.i); AIC(mod.tv.pc.sm)
AIC(mod.tv.pc.sm2); AIC(mod.tv.pc.sm2.r)
```

## Sensitivity check for batch adjustment

```{r pammbatchtest}
# sensitivity analysis, just to verify that batch adjustment is sufficient:
# limiting batches to 4 large ones that have no differences in offset or such
summary(bam(ped_status ~ ti(tend,bs='cr',k=11) + GTcat + ti(tend, by=as.ordered(GTcat),bs='cr') +
               poly(MAGE,2) + FAAR + AA87 + KJONN + MISD + PARITET_5,
            data=filter(ped, BATCH %in% c("M12A", "M24", "ROT1", "ROT2")),
            offset=offset, family=poisson()))
```
## Sensitivity check for covariate adjustment
```{r pammcovars, fig.height=3}
# best model but w/o covariates
summary(bam(ped_status ~ ti(tend,bs='cr',k=11) + GTcat + ti(tend, by=as.ordered(GTcat),bs='cr') +
              BATCH, data=ped, offset=offset, family=poisson()))
summary(mod.tv.pc.sm2)  # no change on the GT coef or the smooth degrees

# adding smooths on some covariates
mod.s1 = bam(ped_status ~ ti(tend,bs='cr',k=11) + GTcat + ti(tend, by=as.ordered(GTcat),bs='cr') +
               BATCH + poly(MAGE,2) + FAAR + AA87 + MISD + PARITET_5 + 
               KJONN +  ti(tend, by=as.ordered(KJONN),bs='cr'),
               data=ped, offset=offset, family=poisson())
summary(mod.s1)  # no change on the GT coef or the smooth degrees

gg_slice(ped, mod.s1, "KJONN", tend=unique(tend), KJONN=factor(1:2)) + theme_bw()
# notably, these factors have quite some nonlinearities too
```


```{r pammplots}
# ---------------
# Useful for printing the summaries:
# library(itsadug)
# gamtabs(mod.tv.pc.sm2.r, type="HTML")
# report_stats(mod.tv.pc.sm2.r)

# ---------------
# VIZ
# (this will nicely default to the most common level for batch)
par(mfrow=c(1,2))
vis.gam(mod.tv.pc.i, view=c("GT", "tend"), theta=40)
vis.gam(mod.tv.pc.sm, view=c("GT", "tend"), theta=40)
vis.gam(mod.tv.pc.sm2, view=c("GTcat", "tend"), theta=40)
vis.gam(mod.tv.pc.sm2.r, view=c("GTcatRev", "tend"), theta=40)
```

```{r pammplots2, fig.height=3}
# (this should also default to the most common level for batch)
# will plot 2SE intervals by default
p1 = gg_slice(ped, mod.pc, "GT", tend=unique(tend), GT=0:2) + theme_bw() + 
  theme(legend.position = "none")
p2 = gg_slice(ped, mod.tv.pc.i, "GT", tend=unique(tend), GT=0:2) + theme_bw()
plot_grid(p1, p2, rel_widths=c(0.45,0.55))

p1 = gg_slice(ped, mod.tv.pc.sm, "GT", tend=unique(tend), GT=0:2) + theme_bw() +
  theme(legend.position = "none")
p2 = gg_slice(ped, mod.tv.pc.sm, "GT", tend=unique(tend), GT=0:2) + theme_bw() +
  coord_cartesian(ylim=c(-2,2))
plot_grid(p1, p2, rel_widths=c(0.45,0.55))

p1 = gg_slice(ped, mod.tv.pc.sm2, "GTcat", tend=unique(tend), GTcat=factor(0:2)) +
  theme_bw() + theme(legend.position = "none")
p2 = gg_slice(ped, mod.tv.pc.sm2, "GTcat", tend=unique(tend), GTcat=factor(0:2)) +
  theme_bw() + coord_cartesian(ylim=c(-2,2))
plot_grid(p1, p2, rel_widths=c(0.45,0.55))

p1 = gg_slice(ped, mod.tv.pc.sm2.r, "GTcatRev", tend=unique(tend), GTcatRev=factor(0:2)) +
  theme_bw() + theme(legend.position = "none")
p2 = gg_slice(ped, mod.tv.pc.sm2.r, "GTcatRev", tend=unique(tend), GTcatRev=factor(0:2)) +
  theme_bw() + coord_cartesian(ylim=c(-2,2))
plot_grid(p1, p2, rel_widths=c(0.45,0.55))

```