# Bachleor's Thesis - Functional ANOVA
# Author: Viktor Dolník
# Script for demonstrating the functional one-way ANOVA test with graphical interpretation


# Library by Myllymäki, et al. containing the rank envelope test
library(GET);

#Function generator for demonstration purposes
demo_1 <- function(i,t)       #H_0 is true, Brownian error, unequal variances
{
  result = t/length(t);
  err = rnorm(result,mean=0, sd=i/50);
  err = cumsum(err);
  result = result*(1-result) + err;
  return(result);
}
demo_2 <- function(i,t)       #H_0 is false, Brownian error, unequal variances
{
  result = t/length(t);
  err = rnorm(result,mean=0, sd=i/50);
  err = cumsum(err);
  result = (result**i)*(1-result**(6-i)) + err;
  return(result);
}

graphical_fanova <- function(t, ni, obs_samples, perm_count=2499, is_contrasts=FALSE, var_correction=FALSE, alpha=0.05)
{
  #First test vector - mean functions of K groups
  means <- function(group_id)
  {
    vect = vector(length=K*length(t));
    for (i in (1:K))
    {
      for (l in t)
      {
        vect[(i-1)*length(t)+l] = sum(obs_samples[group_id == i,l])/ni[i];
      }
    }
    return(vect)
  }
  
  #Second test vector - differences of mean functions, (K of 2) combinations
  contrasts <- function(group_id)
  {
    len = 0; pom1 = c(); pom2 = c();
    vect = vector(length=(K*(K-1)/2)*length(t));
    for (i in (1:(K-1)))
    {
      for (j in ((i+1):K))
      {
        for (l in t)
        {
          pom1[l] = sum(obs_samples[group_id==i,l])/ni[i];
          pom2[l] = sum(obs_samples[group_id==j,l])/ni[j];
        }
        pom3 = sqrt(var(pom1)+var(pom2));
        vect[(len*length(t)+1):((len+1)*length(t))]=(pom1-pom2)/pom3;
        len = len + 1;
      }
    }
    return(vect)
  }
  
  K = length(ni);
  N = sum(ni);
  group_id <- factor(rep(1:K,ni));
  
  #Correction of the variance of observations if var_correction == TRUE
  if (var_correction==TRUE)
  {
    mod_obs_samples <- array(dim=c(N,length(t)));
    overall_sample_mean = apply(obs_samples[,],2,mean);
    overall_sample_sd = apply(obs_samples[,],2,sd);
    index = 0;
    for (i in (1:K))
    {
      group_sample_sd = apply(obs_samples[(index+1):(index+ni[i]),],2,sd);
      for (j in (1:ni[i]))
      {
        mod_obs_samples[index+j,] = (obs_samples[index+j,] - overall_sample_mean) *
        (overall_sample_sd/group_sample_sd) + overall_sample_mean;
      }
      index = index + ni[i];
    }
    obs_samples = mod_obs_samples;
  }
  
  #Computes test statistic, allocates permutations matrix
  if(is_contrasts)
  {
    teststat = contrasts(group_id);
    permutations = array(dim=c(K*(K-1)/2*length(t), perm_count));
    t_new = 1:(length(t)*(K*(K-1))/2);
  } else
  {
    teststat = means(group_id);
    permutations = array(dim=c(K*length(t), perm_count));
    t_new = 1:(length(t)*K);
  }
  
  #Calculates permutations
  for (p in (1:perm_count))
  {
    group_id_perm = sample(group_id, size=N, replace=FALSE);
    if(is_contrasts)
    {
      permutations[,p] = contrasts(group_id_perm);
    } else
    {
      permutations[,p] = means(group_id_perm);
    }
  }
  
  #Transforms the data into a data type used in GET package
  cs_obs_samples <- list(t_new, teststat, permutations);
  names(cs_obs_samples) <- c('r','obs','sim_m');
  cs_obs_samples <- create_curve_set(cs_obs_samples);
  
  #Result is computed using rank_envelope function from GET package
  return(rank_envelope(cs_obs_samples));
}

#INPUT
t = 1:100;
ni=c(10,10,10,10,10,10);
demo_obs <- array(dim=c(sum(ni),length(t))); #every row represents one random function

#random sample allocation
set.seed(2018);
index = 0;
for (i in (1:length(ni)))
{
  for (j in (1:ni[i]))
  {
    #Change to demo_2
    demo_obs[index+j,] = demo_1(i,t);
  }
  index = index + ni[i];
}
set.seed(2018);
result = graphical_fanova(t, ni, demo_obs, perm_count = 2499, is_contrasts=FALSE, var_correction=TRUE);
plot(result);



# Script for simulation study, which is located at the end of Chapter 4

# Libraries containing the graphical test and the asymptotic F-test
library(GET);
library(fda.usc);

#Models for simulation study
model_1 <- function(i,t,stdev,is_brownian)  #H_0 true
{
  result = t/length(t);
  err = rnorm(result,mean=0, sd=stdev);
  if (is_brownian) err = cumsum(err);
  result = result*(1-result) + err;
  return(result);
}
model_2 <- function(i,t,stdev,is_brownian)  #H_0 false
{
  result = t/length(t);
  err = rnorm(result,mean=0, sd=stdev);
  if (is_brownian) err = cumsum(err);
  result = (result**i)*(1-result**(6-i)) + err;
  return(result);
}
model_3 <- function(i,t,stdev,is_brownian)  #H_0 false
{
  result = t/length(t);
  err = rnorm(result,mean=0, sd=stdev);
  if (is_brownian) err = cumsum(err);
  result = (result**(i/5))*(1-result**(6-i/5)) + err;
  return(result);
}
model_4 <- function(i,t,stdev,is_brownian)  #H_0 false
{
  result = t/length(t);
  err = rnorm(result,mean=0, sd=stdev);
  if (is_brownian) err = cumsum(err);
  result = 1+i/50 + err;
  return(result);
}

#Function for calculating the rate of rejection in num_sim simulations
simulate <- function(t,ni,sigma, model, is_brownian, num_sim=100, alpha=0.05)
{
  #Error terms fixed for one vector of sigma values
  set.seed(2018);
  p = array(dim=c(5, length(sigma)));
  p[,] = 0;
  for(r in (1:length(sigma)))
  {
    for(q in (1:num_sim))
    {
      obs_samples <- array(dim=c(sum(ni),length(t)));
      K = length(ni);
      group_id = factor(rep(1:K,ni));
      #Allocation of observations
      index = 0;
      for (i in (1:K))
      {
        for (j in (1:ni[i]))
        {
          obs_samples[index+j,] = model(i,t,sigma[r],is_brownian);
        }
        index = index + ni[i];
      }
      #For graphical fANOVA, we use the GET package
      gfanova = graph.fanova(2000,obs_samples,group_id,t,"equal","means");
      p[1,r] = p[1,r] + as.numeric(attributes(gfanova$ranktest)$p_interval[1] < alpha);
      p[2,r] = p[2,r] + as.numeric(attributes(gfanova$ranktest)$p < alpha);
      gfanova = graph.fanova(2000,obs_samples,group_id,t,"equal","contrasts");
      p[3,r] = p[3,r] + as.numeric(attributes(gfanova$ranktest)$p_interval[1] < alpha);
      p[4,r] = p[4,r] + as.numeric(attributes(gfanova$ranktest)$p < alpha);
      #For the asymptotic F-test, we use the fda.usc package
      p[5,r] = p[5,r] + as.numeric(anova.onefactor(fdata(obs_samples),group_id, nboot=300)$pvalue < alpha);
    }
    p[,r] = p[,r]/num_sim;
  }
  #Output
  cat(sprintf(paste("Sigma: ",paste(sigma[],collapse=" ")," \n")));
  cat(sprintf(paste("AsF: ", paste(p[5,],collapse=" "), "\n")));
  cat(sprintf(paste("GFAM-lib: ", paste(p[1,],collapse=" "), "\n")));
  cat(sprintf(paste("GFAM-mid: ", paste(p[2,],collapse=" "), "\n")));
  cat(sprintf(paste("GFAC-lib: ", paste(p[3,],collapse=" "), "\n")));
  cat(sprintf(paste("GFAC-mid: ", paste(p[4,],collapse=" "), "\n")));
}

t = 1:100;
ni = c(10,10,10);
sigma = c(0.05,0.1,0.5,1,2,5);
simulate(t, ni, sigma, mode = model_1, is_brownian = FALSE, num_sim = 100);