---
title: "MD3F_Simulation_Code"
author: "Alexander V. Alekseyenko"
date: "`r date()`"
output: 
  word_document: default
  html_document: default
---

```{r}
library(ade4)
```

# Simulation
This function simulates a $d$-dimensional sample of $n$ observations. Assuming that the starting position is the origin. $s$ is the speed of drift away from the starting position. Without loss of generality the drift is in equal magnitude along every dimension of the space. $t$ is the time interval that elapses for each Post observation. $sigman$ is the parameter that controls per observation diffusion rate. In most of the following examples $sigman$ is constant for all observations.
```{r}
simulatePrePostS1 = function(n, s=0, d=1, t = rep(1,n), sigmad = rep(1,d)){
  m = s*t/sqrt(d)
  post = sapply(sigmad, function(s2) rnorm(n, mean=m, sd=sqrt(s2)) )
  data.frame(t=c(0, cumsum(t)), 
             X=rbind(rep(0, d), apply(post, 2, cumsum)) )
}

sim2dXdT.interpoint = function(sdf){
  dT = diff(sdf$t)
  dX = apply(apply(as.matrix(sdf[,-1]), 2, diff), 1, function(r) sqrt(sum(r^2)))
  list(dT=dT, dX=dX)
}

```

MD3F and alternative methods:
```{r}
MD3F_linear_estimator = function(sdf){
  sdXdT = sim2dXdT.interpoint(sdf)
  dX = sdXdT$dX
  dT = sdXdT$dT
  est = lm(dX~dT)
  summary(est)$coef[2,c(1,4)]
}

MD3F_linear_estimator2 = #squared distances and time. 
  function(sdf){
  sdXdT = sim2dXdT.interpoint(sdf)
  dX = sdXdT$dX^2
  dT = sdXdT$dT^2
  est = lm(dX~dT)
  summary(est)$coef[2,c(1,4)]
}

univariate_Raw_estimator = function(sdf){
  apply(sdf[,-1, drop=F], 2, function(x) summary(lm(x~sdf$t))$coef[2, c(1,4)])
}

univariate_FDR_estimator = function(uRaw){
  q = p.adjust(uRaw[2,], method="fdr")
  est = sqrt(sum(uRaw[1,]^2))
  c(est, min(q)) 
}

univariate_Holm_estimator = function(uRaw){
  q = p.adjust(uRaw[2,], method="holm")
  est = sqrt(sum(uRaw[1,]^2))
  c(est, min(q))  
}


pcoa_Raw_estimator = function(sdf){
  sdf.pco = data.frame(
    t = sdf$t,
    X = dudi.pco(dist(sdf[,-1,drop=F]), full=T, scannf=F)$li)
  univariate_Raw_estimator(sdf.pco)
}

pcoa_N_estimator = function(pcoRaw, n = 1){
  if(n > ncol(pcoRaw))
    n = ncol(pcoRaw)
  q = p.adjust(pcoRaw[2, 1:n], method="bonf")
  est = sqrt(sum(pcoRaw[1, 1:n]^2))
  c(est, min(q))      
}

set.seed(57.55)
n=#number of simulation replicates; 10000 was used in reported results. 
dt = rexp(n)

sdf = simulatePrePostS1(n, d=4, t=dt)
pcoRaw = pcoa_Raw_estimator(sdf)
pcoa_N_estimator(pcoRaw, 1)
pcoa_N_estimator(pcoRaw, 2)
pcoa_N_estimator(pcoRaw, ncol(pcoRaw))

uRaw = univariate_Raw_estimator(sdf)
univariate_FDR_estimator(uRaw)
univariate_Holm_estimator(uRaw)

MD3F_linear_estimator(sdf)
MD3F_scaled_linear_estimator(sdf)
```

```{r}
simRes = as.data.frame(simRes)
colnames(simRes) = c("s", "d", "n", "repl", 
                     "PCo1.est", "PCo1.p", 
                     "PCo2.est", "PCo2.p",
                     "PCoA.est", "PCoA.p",
                     "uFDR.est", "uFDR.p",
                     "uHolm.est", "uHolm.p",
                     "MD3Flinear.est", "MD3Flinear.p",
                     "MD3Flinear2.est", "MD3Flinear2.p")
```





