graph.fanova: One-way graphical functional ANOVA

Description

One-way ANOVA tests for functional data with graphical interpretation

Usage

graph.fanova(
  nsim,
  curve_set,
  groups,
  variances = "equal",
  contrasts = FALSE,
  n.aver = 1L,
  mirror = FALSE,
  savefuns = FALSE,
  test.equality = c("mean", "var", "cov"),
  cov.lag = 1,
  ...
)

Arguments

nsim

The number of random permutations.

curve_set

The original data (an array of functions) provided as a curve_set object (see create_curve_set) or a fdata object (see fdata). The curve set should include the argument values for the functions in the component r, and the observed functions in the component obs.

groups

The original groups (a factor vector representing the assignment to groups).

variances

Either "equal" or "unequal". If "unequal", then correction for unequal variances as explained in details will be done.

contrasts

Logical. FALSE and TRUE specify the two test functions as described in description part of this help file.

n.aver

If variances = "unequal", there is a possibility to use variances smoothed by appying moving average to the estimated sample variances. n.aver determines how many values on each side do contribute (incl. value itself).

mirror

The complement of the argument circular of filter.

savefuns

Logical. If TRUE, then the functions from permutations are saved to the attribute simfuns.

test.equality

A character with possible values mean (default), var and cov. If mean, the functional ANOVA is performed to compare the means in the groups. If var, then the equality of variances of the curves in the groups is tested by performing the graphical functional ANOVA test on the functions $$Z_{ij}(r) = T_{ij}(r) - \bar{T}_j(r).$$ If cov, then the equality of lag cov.lag covariance is tested by performing the fANOVA with $$W_{ij}(r) = \sqrt{|V_{ij}(r)|\cdot sign(V_{ij}(r))},$$ where $$V_{ij}(r) = (T_{ij}(r) - \bar{T}_j(r))((T_{ij}(r+s) - \bar{T}_j(r+s))).$$ See Mrkvicka et al. (2018) for more details.

cov.lag

The lag of the covariance for testing the equality of covariances, see test.equality.

...

Additional parameters to be passed to global_envelope_test.

Details

This functions can be used to perform one-way graphical functional ANOVA tests described in Mrkvi<U+010D>ka et al. (2016).

The tests assume that there are $J$ groups which contain $n_1,\dots,n_J$ functions $T_{ij}, i=\dots,J, j=1,\dots,n_j$. The functions should be given in the argument curve_set, and the groups in the argument groups. The tests assume that $T_{ij}, i=1,...,n_j$ is an iid sample from a stochastic process with mean function $\mu_j$ and covariance function $\gamma_j(s,t)$ for s,t in R and j = 1,..., J.

If you want to test the hypothesis $$H_0 : \mu_j(r) \equiv 0, j=1, \dots , J,$$ then you should use the test function $$\mathbf{T} = (\overline{T}_1({\bf r}), \overline{T}_2({\bf r}), \dots , \overline{T}_J({\bf r}))$$ where $\overline{T}_i({\bf r})$ is a vector of mean values of functions in the group j. This test function is used when contrasts = FALSE (default).

An alternative is to test the equivalent hypothesis $$H_0 : \mu_i(r) - \mu_j(r) = 0, i=1,\dots,J-1, j=1,\dots,J.$$ This test corresponds to the post-hoc test done usually after an ANOVA test is significant, but it can be directed tested by mean of the combined rank test (Mrkvi<U+010D>ka et al., 2017), if the test vector is taken to consist of the differences of the group averages of test functions, namely $$\mathbf{T'} = (\overline{T}_1({\bf r})-\overline{T}_2({\bf r}), \overline{T}_1({\bf r})-\overline{T}_3({\bf r}), \dots , \overline{T}_{J-1}({\bf r})-\overline{T}_J({\bf r})).$$ With the option contrasts = TRUE the test will be based on this test vector.

The test as such assumes that the variances are equal across the groups of functions. To deal with unequal variances, the differences are rescaled as the first step as follows $$S_{ij}(r) = \frac{T_{ij}(r) - \overline{T}(r))}{\sqrt{Var(T_j(r))}} \sqrt{Var(T(r))} + \overline{T}(r))$$ where $\overline{T}({\bf r})$ is the overall sample mean and $\sqrt{Var(T(r))}$ is the overall sample standard deviation. This scaling of the test functions can be obtained by giving the argument variances = "unequal".

References

Mrkvi<U+010D>ka, T., Hahn, U. and Myllym<U+00E4>ki, M. A one-way ANOVA test for functional data with graphical interpretation. arXiv:1612.03608 [stat.ME] (http://arxiv.org/abs/1612.03608)

Mrkvi<U+010D>ka, T., Myllym<U+00E4>ki, M., and Hahn, U. (2017). Multiple Monte Carlo testing, with applications in spatial point processes. Statistics and Computing 27 (5): 1239-1255. doi:10.1007/s11222-016-9683-9

Myllym<U+00E4>ki, M and Mrkvi<U+010D>ka, T. (2019). GET: Global envelopes in R. arXiv:1911.06583 [stat.ME]

Examples

Run this code

# NOT RUN {
#-- NOx levels example (see for details Myllymaki and Mrkvicka, 2019)
if(require("fda.usc", quietly=TRUE)) {
  # Prepare data
  data(poblenou)
  Free <- poblenou$df$day.festive == 1 |
    as.integer(poblenou$df$day.week) >= 6
  MonThu <- poblenou$df$day.festive == 0 & poblenou$df$day.week %in% 1:4
  Friday <- poblenou$df$day.festive == 0 & poblenou$df$day.week == 5
  Type <- vector(length=length(Free))
  Type[Free] <- "Free"
  Type[MonThu] <- "MonThu"
  Type[Friday] <- "Fri"
  Type <- factor(Type, levels = c("MonThu", "Fri", "Free"))
# }
# NOT RUN {
  # Plot of data
  if(requireNamespace("ggplot2", quietly=TRUE)) {
    df <- do.call(rbind, lapply(1:24, FUN = function(x) {
      data.frame(Hour = x, NOx = poblenou[['nox']]$data[,x],
                 Type = Type, Date = rownames(poblenou[['nox']]$data))
    }))
    ggplot2::ggplot(df) + ggplot2::geom_line(ggplot2::aes(x = Hour, y = NOx, group = Date)) +
      ggplot2::facet_wrap(ggplot2::vars(Type)) + GET:::ThemePlain()
  }
# }
# NOT RUN {
  # Graphical functional ANOVA
  cset <- create_curve_set(list(r=0:23,
             obs=t(log(poblenou[['nox']][['data']]))))
# }
# NOT RUN {
  res.c <- graph.fanova(nsim = 2999, curve_set = cset,
                        groups = Type, variances = "unequal",
                        contrasts = TRUE)
# }
# NOT RUN {
  plot(res.c, xlab = "Hour", ylab = "Diff.")
}

#-- Centred government expenditure centralization ratios example
# This is an example analysis of the centred GEC in Mrkvicka et al.
data(cgec)

# Number of simulations
# }
# NOT RUN {
nsim <- 2499 # increase to reduce Monte Carlo error

# Test for unequal lag 1 covariances
res.cov1 <- graph.fanova(nsim = nsim, curve_set = cgec,
                         groups = attr(cgec, "group"),
                         test.equality = "cov", cov.lag = 1)
plot(res.cov1, ncol=3,
     labels = paste("Group ", 1:3, sep=""),
     xlab=substitute(paste(i, " (", italic(j), ")", sep=""), list(i="Year", j="r")),
     ylab=expression(italic(bar(W)[i](r))))
# Test for equality of variances among groups
res.var <- graph.fanova(nsim = nsim, curve_set = cgec,
                        groups = attr(cgec, "group"),
                        test.equality = "var")
plot(res.var, ncol=3,
     labels = paste("Group ", 1:3, sep=""),
     xlab=substitute(paste(i, " (", italic(j), ")", sep=""), list(i="Year", j="r")),
     ylab=expression(italic(bar(Z)[i](r))))

# Test for equality of means assuming equality of variances
# a) using 'means'
res <- graph.fanova(nsim = nsim, curve_set = cgec,
                    groups = attr(cgec, "group"),
                    variances = "equal",
                    contrasts = FALSE)
plot(res, ncol=3,
     labels = paste("Group ", 1:3, sep=""),
     xlab=substitute(paste(i, " (", italic(j), ")", sep=""), list(i="Year", j="r")),
     ylab=expression(italic(bar(T)[i](r))))
# b) using 'contrasts'
res2 <- graph.fanova(nsim = nsim, curve_set = cgec,
                    groups = attr(cgec, "group"),
                    variances = "equal",
                    contrasts = TRUE)
# }
# NOT RUN {
plot(res2, ncol=3,
     xlab=substitute(paste(i, " (", italic(j), ")", sep=""), list(i="Year", j="r")),
     ylab=expression(italic(bar(T)[i](r)-bar(T)[j](r))))

# }
# NOT RUN {
#-- Rimov water temperatures example
# This is an example analysis of the water temperature data set
# in Mrkvicka et al. (arXiv:1612.03608v2).
data(rimov)
groups <- factor(c(rep(1, times=12), rep(2, times=12), rep(3, times=12)))
nsim <- 999

# Test for equality of variances in the groups
resV <- graph.fanova(nsim=nsim, curve_set=rimov, groups=groups, contrasts = FALSE,
                     test.equality="var")
plot(resV)
# Test for equality of lag 1 covariances in the groups
resC <- graph.fanova(nsim=nsim, curve_set=rimov, groups=groups, contrasts = FALSE,
                     test.equality="cov", cov.lag=1)
plot(resC)

# Test the equality of means in the groups (fANOVA), assuming equality of variances
res <- graph.fanova(nsim=nsim, curve_set=rimov, groups=groups, contrasts = FALSE)
plot(res)
res2 <- graph.fanova(nsim=nsim, curve_set=rimov, groups=groups, contrasts = TRUE)
plot(res2)
# }