Nothing
R/fanova_tests.R
In fdANOVA: Analysis of Variance for Univariate and Multivariate Functional Data
fanova.tests = function(x = NULL, group.label, test = "ALL", params = NULL,
parallel = FALSE, nslaves = NULL){
if(any(is.na(group.label))){ stop("argument group.label can not contain NA values") }
if(test[1] == "ALL"){
test = c("FP", "CH", "CS", "L2N", "L2B", "L2b", "FN", "FB", "Fb", "GPF", "Fmaxb", "TRP")
}else{
for(i in 1:length(test)){
if(!(test[i] %in% c("FP", "CH", "CS", "L2N", "L2B", "L2b", "FN", "FB", "Fb", "GPF", "Fmaxb", "TRP"))){
stop("argument test must be the subvector of c('FP', 'CH', 'CS', 'L2N', 'L2B', 'L2b', 'FN', 'FB', 'Fb', 'GPF', 'Fmaxb', 'TRP')")
}
}
}
group.label0 = unique(group.label)
l = length(group.label0)
n.i = numeric(l)
for(i in 1:l) n.i[i] = sum(group.label == group.label0[i])
anova.statistic.quick = function(x, group.label){
n = length(c(x))
a = length(unique(group.label))
y = data.frame(cbind(x, group.label))
colnames(y) = c("data", "group.label")
means = (doBy::summaryBy(data ~ group.label, data = y))[,2]
mean.g = mean(x)
SSA = sum(table(group.label) * ((means - mean.g)^2))
SST = sum((x - mean.g)^2)
SSE = SST - SSA
return((n-a)*SSA/SSE/(a-1))
}
ATS.simple = function(x, group.label){
N = length(x); n = table(group.label); ng = length(as.vector(n))
means = matrix(tapply(x, group.label, mean), ncol = 1)
MM = diag(1, ng) - matrix(1, ncol = ng, nrow = ng)/ng
DM = diag(diag(MM))
SN = N * diag(as.vector((tapply(x, group.label, var) * (n - 1)/n^2)))
Delta = diag(as.vector(1/(n - 1)))
test.stat = N * t(means) %*% MM %*% means/sum(diag(MM %*% SN))
f1 = sum(diag(MM %*% SN))^2/sum(diag(MM %*% SN %*% MM %*% SN))
f0 = sum(diag(MM %*% SN))^2/sum(diag(DM %*% DM %*% SN %*% SN %*% Delta))
return(c(test.stat, 1 - pf(test.stat, f1, f0)))
}
WTPSp = function(x, group.label, group.label0, n, n.i, l, perm.WTPS, nrTRP){
means = matrix(0, nrow = l, ncol = 1); variances = numeric(l)
for(i in 1:l){
means[i,] = mean(x[group.label == group.label0[i]])
variances[i] = var(x[group.label == group.label0[i]])
}
V.N = n*diag(variances/n.i)
TT = diag(rep(1, l)) - (1/l)*matrix(1, ncol = l, nrow = l)
Q.N.T = n * t(means) %*% TT %*% MASS::ginv(TT %*% V.N %*% TT) %*% TT %*% means
A = t(rep(1/n.i[1], n.i[1]))
A1 = t(rep(1, n.i[1]))
for(i in 2:l){
A = magic::adiag(A, t(rep(1/n.i[i], n.i[i])))
A1 = magic::adiag(A1, t(rep(1, n.i[i])))
}
xperm = matrix(x[perm.WTPS], ncol = nrTRP)
meansP = A %*% xperm
varsP = (A1 %*% (xperm^2) - n.i * meansP^2)/(n.i * (n.i - 1))
statistics.p = sapply(1:nrTRP, function(arg){
TP = TT %*% MASS::ginv(TT %*% (n * diag(varsP[, arg])) %*% TT) %*% TT
statistics.p = diag(n * t(meansP[, arg]) %*% TP %*% meansP[, arg])
})
return(mean(statistics.p >= c(Q.N.T)))
}
if(any(c("FP", "CH", "CS", "L2b", "Fb", "Fmaxb", "TRP") %in% test)){
if(!parallel){
parallel.method = "parallel.method0"
}else{
if(!("doParallel" %in% rownames(installed.packages()))){
stop("Please install package 'doParallel'")
}
# require(foreach, quietly = TRUE)
requireNamespace("foreach", quietly = TRUE)
nlp = parallel::detectCores()
if(is.null(nslaves)){
if(nlp >= 2){
nslaves = nlp
parallel.method = "parallel.method1"
}else{
parallel.method = "parallel.method0"
}
}else{
if(nlp >= 2){
if(nslaves >= 2){
nslaves = nslaves
}else{
nslaves = nlp
}
parallel.method = "parallel.method1"
}else{
parallel.method = "parallel.method0"
}
}
if(parallel.method == "parallel.method1"){
cl = parallel::makePSOCKcluster(nslaves)
doParallel::registerDoParallel(cl)
}
}
}
if("FP" %in% test){
if(is.null(params$paramFP$basis)){ basis = "Fourier" }else{ basis = params$paramFP$basis }
if(!(basis %in% c("Fourier", "b-spline", "own"))){ stop("argument params$paramFP$basis must be one of the following: 'Fourier', 'b-spline', 'own'") }
if(is.null(params$paramFP$B.FP)){ nrFP = 1000 }else{ nrFP = params$paramFP$B.FP }
if(nrFP < 1){ stop("invalid number of permutations (params$paramFP$B.FP)") }
if(basis == "Fourier"){
if(!("fda" %in% rownames(installed.packages()))){
stop("Please install package 'fda'")
}
x = as.matrix(x); n = ncol(x); p = nrow(x)
if(is.null(params$paramFP$minK)){ minK = 3 }else{ minK = params$paramFP$minK }
if(minK < 1){ stop("invalid argument params$paramFP$minK") }
if(is.null(params$paramFP$maxK)){ if(p %% 2 == 1){ maxK = p - 2 }else{ maxK = p - 1 }
}else{ maxK = params$paramFP$maxK }
if(maxK < 1){ stop("invalid argument params$paramFP$maxK") }
if(p <= maxK){ if(p %% 2 == 1){ maxK = p - 2 }else{ maxK = p - 1 } }
if((maxK %% 2) == 0){ stop("the maximum number of basis functions (params$paramFP$maxK) is even") }
if(n != length(group.label)){
stop("number of observations (number of columns in x) and number of elements
in vector of group labels (group.label) must be the same")
}
if(is.null(params$paramFP$criterion)){ criterion = "BIC" }else{ criterion = params$paramFP$criterion }
if(!(criterion %in% c("BIC", "eBIC", "AIC", "AICc", "NO"))){ stop("argument params$paramFP$criterion must be one of the following: 'BIC', 'eBIC', 'AIC', 'AICc', 'NO'") }
if(criterion == "eBIC"){
if(is.null(params$paramFP$gamma.eBIC)){ gamma.eBIC = 0.5 }else{ gamma.eBIC = params$paramFP$gamma.eBIC }
if((gamma.eBIC < 0)|(gamma.eBIC > 1)){
stop("argument params$paramFP$gamma.eBIC must belong to [0,1]")
}
}
if(criterion != "NO"){
if(is.null(params$paramFP$commonK)){ commonK = "mode" }else{ commonK = params$paramFP$commonK }
if(!(commonK %in% c("mode", "min", "max", "mean"))){ stop("argument params$paramFP$commonK must be one of the following: 'mode', 'min', 'max', 'mean'") }
if((minK %% 2) == 0){ stop("the minimum number of basis functions (params$paramFP$minK) is even") }
if(is.null(params$paramFP$int)){
if(parallel.method == "parallel.method0"){
v = matrix(0, nrow = n, ncol = length(seq(minK, maxK, 2)))
for(i in seq(minK, maxK, 2)){
fbasis = fda::create.fourier.basis(c(0, p), i)
fdata = fda::smooth.basisPar(1:p, x, fbasis)
v[, (i+2-minK)/2] = switch(criterion, "BIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + i*log(p),
"eBIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + i*(log(p) + 2*gamma.eBIC*log(maxK)),
"AIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*i,
"AICc" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*i + (2 * i * (i + 1))/(n - i - 1))
}
}
if(parallel.method == "parallel.method1"){
v = foreach(i = seq(minK, maxK, 2), .combine = cbind) %dopar%
{
fbasis = fda::create.fourier.basis(c(0, p), i)
fdata = fda::smooth.basisPar(1:p, x, fbasis)
switch(criterion, "BIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + i*log(p),
"eBIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + i*(log(p) + 2*gamma.eBIC*log(maxK)),
"AIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*i,
"AICc" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*i + (2 * i * (i + 1))/(n - i - 1))
}
}
KK = numeric(n)
for(j in 1:n) KK[j] = 2*which(v[j,] == min(v[j,])) + minK - 2
if(commonK == "mode"){
temp = table(as.vector(KK))
K = min(as.numeric(names(temp)[temp == max(temp)]))
}
if(commonK == "min"){ K = min(KK) }
if(commonK == "max"){ K = max(KK) }
if(commonK == "mean"){ if(floor(mean(KK)) %% 2 == 1){ K = floor(mean(KK)) }else{ K = floor(mean(KK))-1 } }
}else{
if(length(params$paramFP$int) != 2){ stop("argument params$paramFP$int must be of length two") }
if(parallel.method == "parallel.method0"){
v = matrix(0, nrow = n, ncol = length(seq(minK, maxK, 2)))
for(i in seq(minK, maxK, 2)){
fbasis = fda::create.fourier.basis(rangeval = c(params$paramFP$int[1], params$paramFP$int[2]), nbasis = i)
fdata = fda::smooth.basisPar(seq(params$paramFP$int[1], params$paramFP$int[2], length = p), x, fbasis)
v[, (i+2-minK)/2] = switch(criterion, "BIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + i*log(p),
"eBIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + i*(log(p) + 2*gamma.eBIC*log(maxK)),
"AIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*i,
"AICc" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*i + (2 * i * (i + 1))/(n - i - 1))
}
}
if(parallel.method == "parallel.method1"){
v = foreach(i = seq(minK, maxK, 2), .combine = cbind) %dopar%
{
fbasis = fda::create.fourier.basis(rangeval = c(params$paramFP$int[1], params$paramFP$int[2]), nbasis = i)
fdata = fda::smooth.basisPar(seq(params$paramFP$int[1], params$paramFP$int[2], length = p), x, fbasis)
switch(criterion, "BIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + i*log(p),
"eBIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + i*(log(p) + 2*gamma.eBIC*log(maxK)),
"AIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*i,
"AICc" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*i + (2 * i * (i + 1))/(n - i - 1))
}
}
KK = numeric(n)
for(j in 1:n) KK[j] = 2*which(v[j,] == min(v[j,])) + minK - 2
if(commonK == "mode"){
temp = table(as.vector(KK))
K = min(as.numeric(names(temp)[temp == max(temp)]))
}
if(commonK == "min"){ K = min(KK) }
if(commonK == "max"){ K = max(KK) }
if(commonK == "mean"){ if(floor(mean(KK)) %% 2 == 1){ K = floor(mean(KK)) }else{ K = floor(mean(KK))-1 } }
}
}else{ K = maxK }
if(is.null(params$paramFP$int)){
fbasis = fda::create.fourier.basis(c(0, p), K)
lfdata = vector("list", l)
mfdata = matrix(0, nrow = K, ncol = n)
for(i in 1:l){
lfdata[[i]] = fda::Data2fd(1:p, x[, group.label == group.label0[i]], fbasis)$coefs
mfdata[, group.label == group.label0[i]] = lfdata[[i]]
}
}else{
if(length(params$paramFP$int) != 2){ stop("argument params$paramFP$int must be of length two") }
fbasis = fda::create.fourier.basis(c(params$paramFP$int[1], params$paramFP$int[2]), K)
lfdata = vector("list", l)
mfdata = matrix(0, nrow = K, ncol = n)
for(i in 1:l){
lfdata[[i]] = fda::Data2fd(seq(params$paramFP$int[1], params$paramFP$int[2], length = p), x[, group.label == group.label0[i]], fbasis)$coefs
mfdata[, group.label == group.label0[i]] = lfdata[[i]]
}
}
a = 0; b = 0; c = 0
for(i in 1:l){
a = a + sum(t(lfdata[[i]]) %*% lfdata[[i]])/n.i[i]
for(j in 1:l) b = b + sum(t(lfdata[[i]]) %*% (lfdata[[j]]))
c = c + sum((lfdata[[i]])^2)
}
b = b/n; statFP0 = (a-b)/(c-a); statFP = statFP0*(n-l)/(l-1)
if(parallel.method == "parallel.method0"){
FP.p = numeric(nrFP)
for(ii in 1:nrFP){
a = numeric(l)
Perm = sample(1:n)
for(j in 1:l){
if(j == 1){ l1 = Perm[1:n.i[1]] }else{ l1 = Perm[(sum(n.i[1:(j-1)])+1):sum(n.i[1:j])] }
a[j] = sum(t(mfdata[, l1]) %*% mfdata[, l1])/n.i[j]
}
FP.p[ii] = (sum(a)-b)/(c-sum(a))
}
}
if(parallel.method == "parallel.method1"){
FP.p = foreach(ii = 1:nrFP, .combine = "c") %dopar%
{
a = numeric(l)
Perm = sample(1:n)
for(j in 1:l){
if(j == 1){ l1 = Perm[1:n.i[1]] }else{ l1 = Perm[(sum(n.i[1:(j-1)])+1):sum(n.i[1:j])] }
a[j] = sum(t(mfdata[, l1]) %*% mfdata[, l1])/n.i[j]
}
(sum(a)-b)/(c-sum(a))
}
}
pvalueFP = mean(FP.p >= statFP0)
}
if(basis == "b-spline"){
if(!("fda" %in% rownames(installed.packages()))){
stop("Please install package 'fda'")
}
x = as.matrix(x); n = ncol(x); p = nrow(x)
if(is.null(params$paramFP$norder)){ norder = 4 }else{ norder = params$paramFP$norder }
if(norder < 1){ stop("invalid argument params$paramFP$norder") }
if(is.null(params$paramFP$minK)){ minK = norder }else{ minK = params$paramFP$minK }
if(minK < 1){ stop("invalid argument params$paramFP$minK") }
if(minK < norder){ minK = norder }
if(is.null(params$paramFP$maxK)){ maxK = p }else{ maxK = params$paramFP$maxK }
if(maxK < 1){ stop("invalid argument params$paramFP$maxK") }
if(p <= maxK){ maxK = p }
if(n != length(group.label)){
stop("number of observations (number of columns in x) and number of elements
in vector of group labels (group.label) must be the same")
}
if(is.null(params$paramFP$criterion)){ criterion = "BIC" }else{ criterion = params$paramFP$criterion }
if(!(criterion %in% c("BIC", "eBIC", "AIC", "AICc", "NO"))){ stop("argument params$paramFP$criterion must be one of the following: 'BIC', 'eBIC', 'AIC', 'AICc', 'NO'") }
if(criterion == "eBIC"){
if(is.null(params$paramFP$gamma.eBIC)){ gamma.eBIC = 0.5 }else{ gamma.eBIC = params$paramFP$gamma.eBIC }
if((gamma.eBIC < 0)|(gamma.eBIC > 1)){
stop("argument params$paramFP$gamma.eBIC must belong to [0,1]")
}
}
if(criterion != "NO"){
if(is.null(params$paramFP$commonK)){ commonK = "mode" }else{ commonK = params$paramFP$commonK }
if(!(commonK %in% c("mode", "min", "max", "mean"))){ stop("argument params$paramFP$commonK must be one of the following: 'mode', 'min', 'max', 'mean'") }
mmm = minK:maxK
if(is.null(params$paramFP$int)){
if(parallel.method == "parallel.method0"){
v = matrix(0, nrow = n, ncol = length(mmm))
for(i in 1:length(mmm)){
fbasis = fda::create.bspline.basis(c(0, p), mmm[i], norder = norder)
fdata = fda::smooth.basisPar(1:p, x, fbasis)
v[, i] = switch(criterion, "BIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + mmm[i]*log(p),
"eBIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + mmm[i]*(log(p) + 2*gamma.eBIC*log(maxK)),
"AIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*mmm[i],
"AICc" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*mmm[i] + (2 * mmm[i] * (mmm[i] + 1))/(n - mmm[i] - 1))
}
}
if(parallel.method == "parallel.method1"){
v = foreach(i = 1:length(mmm), .combine = cbind) %dopar%
{
fbasis = fda::create.bspline.basis(c(0, p), mmm[i], norder = norder)
fdata = fda::smooth.basisPar(1:p, x, fbasis)
switch(criterion, "BIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + mmm[i]*log(p),
"eBIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + mmm[i]*(log(p) + 2*gamma.eBIC*log(maxK)),
"AIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*mmm[i],
"AICc" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*mmm[i] + (2 * mmm[i] * (mmm[i] + 1))/(n - mmm[i] - 1))
}
}
KK = numeric(n)
for(j in 1:n) KK[j] = mmm[which(v[j,] == min(v[j,]))]
if(commonK == "mode"){
temp = table(as.vector(KK))
K = min(as.numeric(names(temp)[temp == max(temp)]))
}
if(commonK == "min"){ K = min(KK) }
if(commonK == "max"){ K = max(KK) }
if(commonK == "mean"){ K = floor(mean(KK)) }
}else{
if(length(params$paramFP$int) != 2){ stop("argument params$paramFP$int must be of length two") }
if(parallel.method == "parallel.method0"){
v = matrix(0, nrow = n, ncol = length(mmm))
for(i in 1:length(mmm)){
fbasis = fda::create.bspline.basis(rangeval = c(params$paramFP$int[1], params$paramFP$int[2]), nbasis = mmm[i], norder = norder)
fdata = fda::smooth.basisPar(seq(params$paramFP$int[1], params$paramFP$int[2], length = p), x, fbasis)
v[, i] = switch(criterion, "BIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + mmm[i]*log(p),
"eBIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + mmm[i]*(log(p) + 2*gamma.eBIC*log(maxK)),
"AIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*mmm[i],
"AICc" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*mmm[i] + (2 * mmm[i] * (mmm[i] + 1))/(n - mmm[i] - 1))
}
}
if(parallel.method == "parallel.method1"){
v = foreach(i = 1:length(mmm), .combine = cbind) %dopar%
{
fbasis = fda::create.bspline.basis(rangeval = c(params$paramFP$int[1], params$paramFP$int[2]), nbasis = mmm[i], norder = norder)
fdata = fda::smooth.basisPar(seq(params$paramFP$int[1], params$paramFP$int[2], length = p), x, fbasis)
switch(criterion, "BIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + mmm[i]*log(p),
"eBIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + mmm[i]*(log(p) + 2*gamma.eBIC*log(maxK)),
"AIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*mmm[i],
"AICc" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*mmm[i] + (2 * mmm[i] * (mmm[i] + 1))/(n - mmm[i] - 1))
}
}
KK = numeric(n)
for(j in 1:n) KK[j] = mmm[which(v[j,] == min(v[j,]))]
if(commonK == "mode"){
temp = table(as.vector(KK))
K = min(as.numeric(names(temp)[temp == max(temp)]))
}
if(commonK == "min"){ K = min(KK) }
if(commonK == "max"){ K = max(KK) }
if(commonK == "mean"){ K = floor(mean(KK)) }
}
}else{ K = maxK }
if(is.null(params$paramFP$int)){
fbasis = fda::create.bspline.basis(c(0, p), K, norder = norder)
lfdata = vector("list", l)
mfdata = matrix(0, nrow = K, ncol = n)
for(i in 1:l){
lfdata[[i]] = fda::Data2fd(1:p, x[, group.label == group.label0[i]], fbasis)$coefs
mfdata[, group.label == group.label0[i]] = lfdata[[i]]
}
}else{
if(length(params$paramFP$int) != 2){ stop("argument params$paramFP$int must be of length two") }
fbasis = fda::create.bspline.basis(c(params$paramFP$int[1], params$paramFP$int[2]), K, norder = norder)
lfdata = vector("list", l)
mfdata = matrix(0, nrow = K, ncol = n)
for(i in 1:l){
lfdata[[i]] = fda::Data2fd(seq(params$paramFP$int[1], params$paramFP$int[2], length = p), x[, group.label == group.label0[i]], fbasis)$coefs
mfdata[, group.label == group.label0[i]] = lfdata[[i]]
}
}
bspline.cross.prod.matrix = fda::inprod(fbasis, fbasis)
a = 0; b = 0; c = 0
for(i in 1:l){
a = a + sum(t(lfdata[[i]]) %*% bspline.cross.prod.matrix %*% lfdata[[i]])/n.i[i]
for(j in 1:l) b = b + sum(t(lfdata[[i]]) %*% bspline.cross.prod.matrix %*% (lfdata[[j]]))
c = c + sum(diag(t(lfdata[[i]]) %*% bspline.cross.prod.matrix %*% lfdata[[i]]))
}
b = b/n; statFP0 = (a-b)/(c-a); statFP = statFP0*(n-l)/(l-1)
if(parallel.method == "parallel.method0"){
FP.p = numeric(nrFP)
for(ii in 1:nrFP){
a = numeric(l)
Perm = sample(1:n)
for(j in 1:l){
if(j == 1){ l1 = Perm[1:n.i[1]] }else{ l1 = Perm[(sum(n.i[1:(j-1)])+1):sum(n.i[1:j])] }
a[j] = sum(t(mfdata[, l1]) %*% bspline.cross.prod.matrix %*% mfdata[, l1])/n.i[j]
}
FP.p[ii] = (sum(a)-b)/(c-sum(a))
}
}
if(parallel.method == "parallel.method1"){
FP.p = foreach(ii = 1:nrFP, .combine = "c") %dopar%
{
a = numeric(l)
Perm = sample(1:n)
for(j in 1:l){
if(j == 1){ l1 = Perm[1:n.i[1]] }else{ l1 = Perm[(sum(n.i[1:(j-1)])+1):sum(n.i[1:j])] }
a[j] = sum(t(mfdata[, l1]) %*% bspline.cross.prod.matrix %*% mfdata[, l1])/n.i[j]
}
(sum(a)-b)/(c-sum(a))
}
}
pvalueFP = mean(FP.p >= statFP0)
}
if(basis == "own"){
own.basis = params$paramFP$own.basis
mfdata = own.basis
n = ncol(mfdata); K = nrow(mfdata)
if(n != length(group.label)){ stop("the number of observations is not equal to the number of elements in vector of labels") }
lfdata = vector("list", l)
for(i in 1:l){ lfdata[[i]] = mfdata[, group.label == group.label0[i]] }
own.cross.prod.mat = params$paramFP$own.cross.prod.mat
if(K != nrow(own.cross.prod.mat)){ stop("invalid argument params$paramFP$own.cross.prod.mat") }
if(K != ncol(own.cross.prod.mat)){ stop("invalid argument params$paramFP$own.cross.prod.mat") }
a = 0; b = 0; c = 0
for(i in 1:l){
a = a + sum(t(lfdata[[i]]) %*% own.cross.prod.mat %*% lfdata[[i]])/n.i[i]
for(j in 1:l) b = b + sum(t(lfdata[[i]]) %*% own.cross.prod.mat %*% (lfdata[[j]]))
c = c + sum(diag(t(lfdata[[i]]) %*% own.cross.prod.mat %*% lfdata[[i]]))
}
b = b/n; statFP0 = (a-b)/(c-a); statFP = statFP0*(n-l)/(l-1)
if(parallel.method == "parallel.method0"){
FP.p = numeric(nrFP)
for(ii in 1:nrFP){
a = numeric(l)
Perm = sample(1:n)
for(j in 1:l){
if(j == 1){ l1 = Perm[1:n.i[1]] }else{ l1 = Perm[(sum(n.i[1:(j-1)])+1):sum(n.i[1:j])] }
a[j] = sum(t(mfdata[, l1]) %*% own.cross.prod.mat %*% mfdata[, l1])/n.i[j]
}
FP.p[ii] = (sum(a)-b)/(c-sum(a))
}
}
if(parallel.method == "parallel.method1"){
FP.p = foreach(ii = 1:nrFP, .combine = "c") %dopar%
{
a = numeric(l)
Perm = sample(1:n)
for(j in 1:l){
if(j == 1){ l1 = Perm[1:n.i[1]] }else{ l1 = Perm[(sum(n.i[1:(j-1)])+1):sum(n.i[1:j])] }
a[j] = sum(t(mfdata[, l1]) %*% own.cross.prod.mat %*% mfdata[, l1])/n.i[j]
}
(sum(a)-b)/(c-sum(a))
}
}
pvalueFP = mean(FP.p >= statFP0)
}
if(basis == "own"){
resultFP = list(statFP = statFP, pvalueFP = pvalueFP, B.FP = nrFP, basis = basis, own.basis = own.basis, own.cross.prod.mat = own.cross.prod.mat)
}else{
if(criterion == "NO"){
if(is.null(params$paramFP$int)){
if(basis == "Fourier"){
resultFP = list(statFP = statFP, pvalueFP = pvalueFP, B.FP = nrFP, basis = basis, criterion = criterion, K = K, maxK = maxK)
}
if(basis == "b-spline"){
resultFP = list(statFP = statFP, pvalueFP = pvalueFP, B.FP = nrFP, basis = basis, criterion = criterion, K = K, maxK = maxK, norder = norder)
}
}else{
if(basis == "Fourier"){
resultFP = list(statFP = statFP, pvalueFP = pvalueFP, int = params$paramFP$int, B.FP = nrFP, basis = basis, criterion = criterion, K = K, maxK = maxK)
}
if(basis == "b-spline"){
resultFP = list(statFP = statFP, pvalueFP = pvalueFP, int = params$paramFP$int, B.FP = nrFP, basis = basis, criterion = criterion, K = K, maxK = maxK, norder = norder)
}
}
}else{
if(is.null(params$paramFP$int)){
if(criterion == "eBIC"){
if(basis == "Fourier"){
resultFP = list(statFP = statFP, pvalueFP = pvalueFP, B.FP = nrFP, basis = basis, criterion = criterion, commonK = commonK, K = K, minK = minK, maxK = maxK, gamma.eBIC = gamma.eBIC)
}
if(basis == "b-spline"){
resultFP = list(statFP = statFP, pvalueFP = pvalueFP, B.FP = nrFP, basis = basis, criterion = criterion, commonK = commonK, K = K, minK = minK, maxK = maxK, norder = norder, gamma.eBIC = gamma.eBIC)
}
}else{
if(basis == "Fourier"){
resultFP = list(statFP = statFP, pvalueFP = pvalueFP, B.FP = nrFP, basis = basis, criterion = criterion, commonK = commonK, K = K, minK = minK, maxK = maxK)
}
if(basis == "b-spline"){
resultFP = list(statFP = statFP, pvalueFP = pvalueFP, B.FP = nrFP, basis = basis, criterion = criterion, commonK = commonK, K = K, minK = minK, maxK = maxK, norder = norder)
}
}
}else{
if(criterion == "eBIC"){
if(basis == "Fourier"){
resultFP = list(statFP = statFP, pvalueFP = pvalueFP, int = params$paramFP$int, B.FP = nrFP, basis = basis, criterion = criterion, commonK = commonK, K = K, minK = minK, maxK = maxK, gamma.eBIC = gamma.eBIC)
}
if(basis == "b-spline"){
resultFP = list(statFP = statFP, pvalueFP = pvalueFP, int = params$paramFP$int, B.FP = nrFP, basis = basis, criterion = criterion, commonK = commonK, K = K, minK = minK, maxK = maxK, norder = norder, gamma.eBIC = gamma.eBIC)
}
}else{
if(basis == "Fourier"){
resultFP = list(statFP = statFP, pvalueFP = pvalueFP, int = params$paramFP$int, B.FP = nrFP, basis = basis, criterion = criterion, commonK = commonK, K = K, minK = minK, maxK = maxK)
}
if(basis == "b-spline"){
resultFP = list(statFP = statFP, pvalueFP = pvalueFP, int = params$paramFP$int, B.FP = nrFP, basis = basis, criterion = criterion, commonK = commonK, K = K, minK = minK, maxK = maxK, norder = norder)
}
}
}
}
}
}
if(any(c("CH", "CS") %in% test)){
if(!("MASS" %in% rownames(installed.packages()))){
stop("Please install package 'MASS'")
}
x = as.matrix(x); n = ncol(x)
if(n != length(group.label)){
stop("number of observations (number of columns in x) and number of elements
in vector of group labels (group.label) must be the same")
}
statCHCS = 0
for(i in 1:(l-1)){
for(j in (i+1):l){
statCHCS = statCHCS + n.i[i]*sum((rowMeans(x[, group.label == group.label0[i]]) -
rowMeans(x[, group.label == group.label0[j]]))^2)
}
}
if("CH" %in% test){
if(is.null(params$paramCH)){ nrCH = 10000 }else{ nrCH = params$paramCH }
if(nrCH < 1){ stop("invalid number of discretized artificial trajectories for each process Z_i(t) (params$paramCH)") }
cov.fun = ((n-1)/(n-l))*var(t(x)); Z = vector("list", l)
for(i in 1:l) Z[[i]] = MASS::mvrnorm(nrCH, numeric(nrow(cov.fun)), cov.fun)
if(parallel.method == "parallel.method0"){
V = numeric(nrCH)
for(m in 1:nrCH){
for(i in 1:(l-1)){
for(j in (i+1):l){
V[m] = V[m] + sum((Z[[i]][m,]-sqrt(n.i[i]/n.i[j])*Z[[j]][m,])^2)
}
}
}
}
if(parallel.method == "parallel.method1"){
V = foreach(m = 1:nrCH, .combine = "c") %dopar%
{
V = 0
for(i in 1:(l-1)){
for(j in (i+1):l){
V = V + sum((Z[[i]][m,]-sqrt(n.i[i]/n.i[j])*Z[[j]][m,])^2)
}
}
V
}
}
resultCH = list(statCH = statCHCS, pvalueCH = mean(V > statCHCS), paramCH = nrCH)
}
if("CS" %in% test){
if(is.null(params$paramCS)){ nrCS = 10000 }else{ nrCS = params$paramCS }
if(nrCS < 1){ stop("invalid number of discretized artificial trajectories for each process Z_i(t) (params$paramCS)") }
Z = vector("list", l)
for(i in 1:l){
cov.fun = var(t(x[, group.label == group.label0[i]]))
Z[[i]] = MASS::mvrnorm(nrCS, numeric(nrow(cov.fun)), cov.fun)
}
if(parallel.method == "parallel.method0"){
V = numeric(nrCS)
for(m in 1:nrCS){
for(i in 1:(l-1)){
for(j in (i+1):l){
V[m] = V[m] + sum((Z[[i]][m,]-sqrt(n.i[i]/n.i[j])*Z[[j]][m,])^2)
}
}
}
}
if(parallel.method == "parallel.method1"){
V = foreach(m = 1:nrCS, .combine = "c") %dopar%
{
V = 0
for(i in 1:(l-1)){
for(j in (i+1):l){
V = V + sum((Z[[i]][m,]-sqrt(n.i[i]/n.i[j])*Z[[j]][m,])^2)
}
}
V
}
}
resultCS = list(statCS = statCHCS, pvalueCS = mean(V > statCHCS), paramCS = nrCS)
}
}
if(any(c("L2N", "L2B", "L2b", "FN", "FB", "Fb", "GPF", "Fmaxb") %in% test)){
x = as.matrix(x); n = ncol(x); p = nrow(x)
if(n != length(group.label)){
stop("number of observations (number of columns in x) and number of elements
in vector of group labels (group.label) must be the same")
}
mu0 = rowMeans(x)
vmu = matrix(0, nrow = l, ncol = p)
z = matrix(0, nrow = n, ncol = p)
SSR = 0; SSE = 0
for(i in 1:l){
xi = x[, group.label == group.label0[i]]
mui = rowMeans(xi); vmu[i,] = mui
zi = t(xi) - as.matrix(rep(1, n.i[i])) %*% mui
if(i==1){ z[1:n.i[i],] = zi }else{ z[(cumsum(n.i)[i-1]+1):cumsum(n.i)[i],] = zi }
SSR = SSR + n.i[i]*(mui - mu0)^2
SSE = SSE + colSums(zi^2)
}
if(n > p){ gamma.z = t(z) %*% z/(n-l) }else{ gamma.z = z %*% t(z)/(n-l) }
A = sum(diag(gamma.z)); B = sum(diag(gamma.z %*% gamma.z))
A2N = A^2; B2N = B
A2B = (n-l)*(n-l+1)/(n-l-1)/(n-l+2)*(A^2-2*B/(n-l+1))
B2B = (n-l)^2/(n-l-1)/(n-l+2)*(B-A^2/(n-l))
if(any(c("L2N", "L2B", "L2b") %in% test)){
statL2 = sum(SSR)
if("L2N" %in% test){
betaL2N = B2N/A; kappaL2N = A2N/B2N
pvalueL2N = 1-pchisq(statL2/betaL2N, (l-1)*kappaL2N)
resultL2N = list(statL2 = statL2, pvalueL2N = pvalueL2N, betaL2N = betaL2N, dL2N = (l-1)*kappaL2N)
}
if("L2B" %in% test){
betaL2B = B2B/A; kappaL2B = A2B/B2B
pvalueL2B = 1-pchisq(statL2/betaL2B, (l-1)*kappaL2B)
resultL2B = list(statL2 = statL2, pvalueL2B = pvalueL2B, betaL2B = betaL2B, dL2B = (l-1)*kappaL2B)
}
if("L2b" %in% test){
if(is.null(params$paramL2b)){ nrL2b = 10000 }else{ nrL2b = params$paramL2b }
if(nrL2b < 1){ stop("invalid number of bootstrap replicates (params$paramL2b)") }
if(parallel.method == "parallel.method0"){
statL2boot = numeric(nrL2b)
for(ii in 1:nrL2b){
vmuboot = matrix(0, nrow = l, ncol = p)
for(i in 1:l){
xi = x[, group.label == group.label0[i]]
xiboot = xi[, floor(runif(n.i[i])*(n.i[i]-1)) + 1]
vmuboot[i,] = rowMeans(xiboot)-vmu[i,]
}
mu0boot = n.i %*% vmuboot/n
SSRboot = 0
for(i in 1:l) SSRboot = SSRboot + n.i[i]*(vmuboot[i,]-mu0boot)^2
statL2boot[ii] = sum(SSRboot)
}
}
if(parallel.method == "parallel.method1"){
statL2boot = foreach(ii = 1:nrL2b, .combine = "c") %dopar%
{
vmuboot = matrix(0, nrow = l, ncol = p)
for(i in 1:l){
xi = x[, group.label == group.label0[i]]
xiboot = xi[, floor(runif(n.i[i])*(n.i[i]-1)) + 1]
vmuboot[i,] = rowMeans(xiboot)-vmu[i,]
}
mu0boot = n.i %*% vmuboot/n
SSRboot = 0
for(i in 1:l) SSRboot = SSRboot + n.i[i]*(vmuboot[i,]-mu0boot)^2
sum(SSRboot)
}
}
pvalueL2b = mean(statL2boot >= statL2)
resultL2b = list(statL2 = statL2, pvalueL2b = pvalueL2b, paramL2b = nrL2b)
}
}
if(any(c("FN", "FB", "Fb") %in% test)){
statF = sum(SSR)/sum(SSE)*(n-l)/(l-1)
if("FN" %in% test){
kappaFN = A2N/B2N
pvalueFN = 1-pf(statF, (l-1)*kappaFN, (n-l)*kappaFN)
resultFN = list(statF = statF, pvalueFN = pvalueFN, d1FN = (l-1)*kappaFN, d2FN = (n-l)*kappaFN)
}
if("FB" %in% test){
kappaFB = A2B/B2B
pvalueFB = 1-pf(statF, (l-1)*kappaFB, (n-l)*kappaFB)
resultFB = list(statF = statF, pvalueFB = pvalueFB, d1FB = (l-1)*kappaFB, d2FB = (n-l)*kappaFB)
}
if("Fb" %in% test){
if(is.null(params$paramFb)){ nrFb = 10000 }else{ nrFb = params$paramFb }
if(nrFb < 1){ stop("invalid number of bootstrap replicates (params$paramFb)") }
if(parallel.method == "parallel.method0"){
statFboot = numeric(nrFb)
for(ii in 1:nrFb){
vmuboot = matrix(0, nrow = l, ncol = p)
zboot = matrix(0, nrow = n, ncol = p)
for(i in 1:l){
xi = x[, group.label == group.label0[i]]
xiboot = xi[, floor(runif(n.i[i])*(n.i[i]-1)) + 1]
muiboot = rowMeans(xiboot)
ziboot = t(xiboot) - as.matrix(rep(1, n.i[i])) %*% muiboot
if(i==1){ zboot[1:n.i[i],] = ziboot }else{ zboot[(cumsum(n.i)[i-1]+1):cumsum(n.i)[i],] = ziboot }
vmuboot[i,] = rowMeans(xiboot)-vmu[i,]
}
mu0boot = n.i %*% vmuboot/n
SSRboot = 0
for(i in 1:l) SSRboot = SSRboot + n.i[i]*(vmuboot[i,]-mu0boot)^2
if(n > p){ gammaboot = t(zboot) %*% zboot/(n-l) }else{ gammaboot = zboot %*% t(zboot)/(n-l) }
statFboot[ii] = sum(SSRboot)/sum(diag(gammaboot))/(l-1)
}
}
if(parallel.method == "parallel.method1"){
statFboot = foreach(ii = 1:nrFb, .combine = "c") %dopar%
{
vmuboot = matrix(0, nrow = l, ncol = p)
zboot = matrix(0, nrow = n, ncol = p)
for(i in 1:l){
xi = x[, group.label == group.label0[i]]
xiboot = xi[, floor(runif(n.i[i])*(n.i[i]-1)) + 1]
muiboot = rowMeans(xiboot)
ziboot = t(xiboot) - as.matrix(rep(1, n.i[i])) %*% muiboot
if(i==1){ zboot[1:n.i[i],] = ziboot }else{ zboot[(cumsum(n.i)[i-1]+1):cumsum(n.i)[i],] = ziboot }
vmuboot[i,] = rowMeans(xiboot)-vmu[i,]
}
mu0boot = n.i %*% vmuboot/n
SSRboot = 0
for(i in 1:l) SSRboot = SSRboot + n.i[i]*(vmuboot[i,]-mu0boot)^2
if(n > p){ gammaboot = t(zboot) %*% zboot/(n-l) }else{ gammaboot = zboot %*% t(zboot)/(n-l) }
sum(SSRboot)/sum(diag(gammaboot))/(l-1)
}
}
pvalueFb = mean(statFboot >= statF)
resultFb = list(statF = statF, pvalueFb = pvalueFb, paramFb = nrFb)
}
}
if("GPF" %in% test){
statGPF = mean(SSR/SSE*(n-l)/(l-1))
z = z/(as.matrix(rep(1, n)) %*% sqrt(colSums(z^2)))
if(n >= p){ z = t(z) %*% z }else{ z = z %*% t(z) }
betaGPF = (sum(diag(z %*% z))/p^2/(l-1))/((n-l)/(n-l-2))
dGPF = ((n-l)/(n-l-2))^2/(sum(diag(z %*% z))/p^2/(l-1))
pvalueGPF = 1-pchisq(statGPF/betaGPF,dGPF)
resultGPF = list(statGPF = statGPF, pvalueGPF = pvalueGPF, betaGPF = betaGPF, dGPF = dGPF)
}
if("Fmaxb" %in% test){
if(is.null(params$paramFmaxb)){ nrFmaxb = 10000 }else{ nrFmaxb = params$paramFmaxb }
if(nrFmaxb < 1){ stop("invalid number of bootstrap replicates (params$paramFmaxb)") }
statFmax = max(SSR/SSE*(n-l)/(l-1))
xs = matrix(0, nrow = n, ncol = p)
for(i in 1:l){ xs[group.label == group.label0[i],] = t(x[, group.label == group.label0[i]]) - as.matrix(rep(1, n.i[i])) %*% vmu[i,] }
if(parallel.method == "parallel.method0"){
statFmaxboot = numeric(nrFmaxb)
for(ii in 1:nrFmaxb){
vmuboot = matrix(0, nrow = l, ncol = p)
zboot = matrix(0, nrow = n, ncol = p)
xboot = xs[sample(1:n, replace = TRUE),]
mu0boot = colMeans(xboot)
SSRboot = 0
for(i in 1:l){
xiboot = xboot[group.label == group.label0[i],]
muiboot = colMeans(xiboot)
vmuboot[i,] = muiboot
ziboot = xiboot - as.matrix(rep(1, n.i[i])) %*% muiboot
if(i==1){ zboot[1:n.i[i],] = ziboot }else{ zboot[(cumsum(n.i)[i-1]+1):cumsum(n.i)[i],] = ziboot }
SSRboot = SSRboot + n.i[i]*(muiboot - mu0boot)^2
}
SSEboot = diag(t(zboot) %*% zboot)
statFmaxboot[ii] = max(SSRboot/SSEboot*(n-l)/(l-1))
}
}
if(parallel.method == "parallel.method1"){
statFmaxboot = foreach(ii = 1:nrFmaxb, .combine = "c") %dopar%
{
vmuboot = matrix(0, nrow = l, ncol = p)
zboot = matrix(0, nrow = n, ncol = p)
xboot = xs[sample(1:n, replace = TRUE),]
mu0boot = colMeans(xboot)
SSRboot = 0
for(i in 1:l){
xiboot = xboot[group.label == group.label0[i],]
muiboot = colMeans(xiboot)
vmuboot[i,] = muiboot
ziboot = xiboot - as.matrix(rep(1, n.i[i])) %*% muiboot
if(i==1){ zboot[1:n.i[i],] = ziboot }else{ zboot[(cumsum(n.i)[i-1]+1):cumsum(n.i)[i],] = ziboot }
SSRboot = SSRboot + n.i[i]*(muiboot - mu0boot)^2
}
SSEboot = diag(t(zboot) %*% zboot)
max(SSRboot/SSEboot*(n-l)/(l-1))
}
}
pvalueFmaxb = mean(statFmaxboot >= statFmax)
resultFmaxb = list(statFmax = statFmax, pvalueFmaxb = pvalueFmaxb, paramFmaxb = nrFmaxb)
}
}
if ("TRP" %in% test) {
if (!("doBy" %in% rownames(installed.packages()))) {
stop("Please install package 'doBY'")
}
if (!("MASS" %in% rownames(installed.packages()))) {
stop("Please install package 'MASS'")
}
if (!("magic" %in% rownames(installed.packages()))) {
stop("Please install package 'magic'")
}
x <- as.matrix(x); n <- ncol(x); p <- nrow(x)
if (n != length(group.label)) {
stop("number of observations (number of columns in x) and number of elements
in vector of group labels (group.label) must be the same")
}
if (is.null(params$paramTRP$projection)) {
projection <- "GAUSS"
}else{
projection <- params$paramTRP$projection
}
if (!(projection %in% c("GAUSS", "BM"))) {
stop("argument params$paramTRP$projection must be one of the following: 'GAUSS', 'BM'")
}
if (is.null(params$paramTRP$permutation)) {
permutationTRP <- FALSE
}else{
permutationTRP <- params$paramTRP$permutation
}
if (! is.logical(permutationTRP)) {
stop("argument permutation is not logical (params$paramTRP$permutation)")
}
if (is.null(params$paramTRP$independent.projection.tests)) {
independent.projection.tests <- TRUE
}else{
independent.projection.tests <- params$paramTRP$independent.projection.tests
}
if (! is.logical(independent.projection.tests)) {
stop("argument independent.projection.tests is not logical
(params$paramTRP$independent.projection.tests)")
}
if (is.null(params$paramTRP$B.TRP)) {
nrTRP <- 10000
}else{
nrTRP <- params$paramTRP$B.TRP
}
if (nrTRP < 1) {
stop("invalid number of permutations (params$paramTRP$B.TRP)")
}
if (is.null(params$paramTRP$k)) {
k <- 30
}else{
k <- params$paramTRP$k
}
if (any(k < 1)) {
stop("invalid number of projections (k)")
}
pvalues.anova <- numeric(length(k))
pvalues.ATS.simple <- numeric(length(k))
pvalues.WTPS <- numeric(length(k))
modulo <- function(z){ sqrt(sum(z^2)) }
if (independent.projection.tests) {
all.data.proj <- list()
iik <- 0
for(ik in k){
ik.data.proj <- matrix(0, nrow = n, ncol = ik)
if(projection == "GAUSS"){
z <- matrix(rnorm(p * ik), nrow = ik, ncol = p)
modu <- apply(z, 1, modulo)
z <- z/modu
if(parallel.method == "parallel.method0"){
anova.p <- numeric(ik)
ATS.simple.p <- numeric(ik)
WTPS.p <- numeric(ik)
for(j in 1:ik){
ik.data.proj[, j] <- t(x) %*% z[j,]
data.proj <- ik.data.proj[, j]
perm.WTPS <- matrix(0, nrow = n, ncol = nrTRP)
for(i.perm in 1:nrTRP){
perm.WTPS[, i.perm] <- sample(1:n)
}
WTPS.p[j] <- WTPSp(data.proj, group.label, group.label0, n, n.i, l, perm.WTPS = perm.WTPS, nrTRP)
if(permutationTRP == FALSE){
anova.p[j] <- 1-pf(anova.statistic.quick(data.proj, group.label), l-1, n-l)
ATS.simple.p[j] <- ATS.simple(data.proj, group.label)[2]
}else{
anova.s <- anova.statistic.quick(data.proj, group.label)
ATS.simple.s <- ATS.simple(data.proj, group.label)[1]
anova.perm <- numeric(nrTRP)
ATS.simple.perm <- numeric(nrTRP)
for(i.perm in 1:nrTRP){
anova.perm[i.perm] <- anova.statistic.quick(data.proj, sample(group.label))
ATS.simple.perm[i.perm] <- ATS.simple(data.proj, sample(group.label))[1]
}
anova.p[j] <- mean(anova.perm >= anova.s)
ATS.simple.p[j] <- mean(ATS.simple.perm >= ATS.simple.s)
}
}
iik <- iik + 1
all.data.proj[[iik]] <- ik.data.proj
pvalues.anova[iik] <- min(ik*anova.p[order(anova.p)]/1:ik)
pvalues.ATS.simple[iik] <- min(ik*ATS.simple.p[order(ATS.simple.p)]/1:ik)
pvalues.WTPS[iik] <- min(ik*WTPS.p[order(WTPS.p)]/1:ik)
}
if(parallel.method == "parallel.method1"){
rs <- foreach(j = 1:ik, .combine = rbind) %dopar%
{
ik.data.proj[, j] <- t(x) %*% z[j,]
data.proj <- ik.data.proj[, j]
perm.WTPS <- matrix(0, nrow = n, ncol = nrTRP)
for(i.perm in 1:nrTRP){
perm.WTPS[, i.perm] <- sample(1:n)
}
WTPS.p <- WTPSp(data.proj, group.label, group.label0, n, n.i, l, perm.WTPS = perm.WTPS, nrTRP)
if(permutationTRP == FALSE){
c(1-pf(anova.statistic.quick(data.proj, group.label), l-1, n-l),
ATS.simple(data.proj, group.label)[2], WTPS.p, data.proj)
}else{
anova.s <- anova.statistic.quick(data.proj, group.label)
ATS.simple.s <- ATS.simple(data.proj, group.label)[1]
anova.perm <- numeric(nrTRP)
ATS.simple.perm <- numeric(nrTRP)
for(i.perm in 1:nrTRP){
anova.perm[i.perm] <- anova.statistic.quick(data.proj, sample(group.label))
ATS.simple.perm[i.perm] <- ATS.simple(data.proj, sample(group.label))[1]
}
c(mean(anova.perm >= anova.s),
mean(ATS.simple.perm >= ATS.simple.s), WTPS.p, data.proj)
}
}
if(ik == 1) rs <- matrix(rs, nrow = 1, ncol = n+3)
iik <- iik + 1
if(ik == 1){
all.data.proj[[iik]] <- matrix(rs[, 4:(n+3)], ncol = 1)
}else{
all.data.proj[[iik]] <- t(rs[, 4:(n+3)])
}
pvalues.anova[iik] <- min(ik*(rs[, 1])[order(rs[, 1])]/1:ik)
pvalues.ATS.simple[iik] <- min(ik*(rs[, 2])[order(rs[, 2])]/1:ik)
pvalues.WTPS[iik] <- min(ik*(rs[, 3])[order(rs[, 3])]/1:ik)
}
}else{
if(parallel.method == "parallel.method0"){
anova.p <- numeric(ik)
ATS.simple.p <- numeric(ik)
WTPS.p <- numeric(ik)
for(j in 1:ik){
bm.p <- cumsum(rnorm(p, mean = 0, sd = 1))/sqrt(p)
bm.p <- bm.p/modulo(bm.p)
ik.data.proj[, j] <- t(x) %*% as.matrix(bm.p)
data.proj <- ik.data.proj[, j]
perm.WTPS <- matrix(0, nrow = n, ncol = nrTRP)
for(i.perm in 1:nrTRP){
perm.WTPS[, i.perm] <- sample(1:n)
}
WTPS.p[j] <- WTPSp(data.proj, group.label, group.label0, n, n.i, l, perm.WTPS = perm.WTPS, nrTRP)
if(permutationTRP == FALSE){
anova.p[j] <- 1-pf(anova.statistic.quick(data.proj, group.label), l-1, n-l)
ATS.simple.p[j] <- ATS.simple(data.proj, group.label)[2]
}else{
anova.s <- anova.statistic.quick(data.proj, group.label)
ATS.simple.s <- ATS.simple(data.proj, group.label)[1]
anova.perm <- numeric(nrTRP)
ATS.simple.perm <- numeric(nrTRP)
for(i.perm in 1:nrTRP){
anova.perm[i.perm] <- anova.statistic.quick(data.proj, sample(group.label))
ATS.simple.perm[i.perm] <- ATS.simple(data.proj, sample(group.label))[1]
}
anova.p[j] <- mean(anova.perm >= anova.s)
ATS.simple.p[j] <- mean(ATS.simple.perm >= ATS.simple.s)
}
}
iik <- iik + 1
all.data.proj[[iik]] <- ik.data.proj
pvalues.anova[iik] <- min(ik*anova.p[order(anova.p)]/1:ik)
pvalues.ATS.simple[iik] <- min(ik*ATS.simple.p[order(ATS.simple.p)]/1:ik)
pvalues.WTPS[iik] <- min(ik*WTPS.p[order(WTPS.p)]/1:ik)
}
if(parallel.method == "parallel.method1"){
rs <- foreach(j = 1:ik, .combine = rbind) %dopar%
{
bm.p <- cumsum(rnorm(p, mean = 0, sd = 1))/sqrt(p)
bm.p <- bm.p/modulo(bm.p)
ik.data.proj[, j] <- t(x) %*% as.matrix(bm.p)
data.proj <- ik.data.proj[, j]
perm.WTPS <- matrix(0, nrow = n, ncol = nrTRP)
for(i.perm in 1:nrTRP){
perm.WTPS[, i.perm] <- sample(1:n)
}
WTPS.p <- WTPSp(data.proj, group.label, group.label0, n, n.i, l, perm.WTPS = perm.WTPS, nrTRP)
if(permutationTRP == FALSE){
c(1-pf(anova.statistic.quick(data.proj, group.label), l-1, n-l),
ATS.simple(data.proj, group.label)[2], WTPS.p, data.proj)
}else{
anova.s <- anova.statistic.quick(data.proj, group.label)
ATS.simple.s <- ATS.simple(data.proj, group.label)[1]
anova.perm <- numeric(nrTRP)
ATS.simple.perm <- numeric(nrTRP)
for(i.perm in 1:nrTRP){
anova.perm[i.perm] <- anova.statistic.quick(data.proj, sample(group.label))
ATS.simple.perm[i.perm] <- ATS.simple(data.proj, sample(group.label))[1]
}
c(mean(anova.perm >= anova.s),
mean(ATS.simple.perm >= ATS.simple.s), WTPS.p, data.proj)
}
}
if(ik == 1) rs <- matrix(rs, nrow = 1, ncol = n+3)
iik <- iik + 1
if(ik == 1){
all.data.proj[[iik]] <- matrix(rs[, 4:(n+3)], ncol = 1)
}else{
all.data.proj[[iik]] <- t(rs[, 4:(n+3)])
}
pvalues.anova[iik] <- min(ik*(rs[, 1])[order(rs[, 1])]/1:ik)
pvalues.ATS.simple[iik] <- min(ik*(rs[, 2])[order(rs[, 2])]/1:ik)
pvalues.WTPS[iik] <- min(ik*(rs[, 3])[order(rs[, 3])]/1:ik)
}
}
}
}else{
ik <- max(k)
ik.data.proj <- matrix(0, nrow = n, ncol = ik)
if (projection == "GAUSS") {
z <- matrix(rnorm(p * ik), nrow = ik, ncol = p)
modu <- apply(z, 1, modulo)
z <- z/modu
if(parallel.method == "parallel.method0"){
anova.p <- numeric(ik)
ATS.simple.p <- numeric(ik)
WTPS.p <- numeric(ik)
for(j in 1:ik){
ik.data.proj[, j] <- t(x) %*% z[j,]
data.proj <- ik.data.proj[, j]
perm.WTPS <- matrix(0, nrow = n, ncol = nrTRP)
for(i.perm in 1:nrTRP){
perm.WTPS[, i.perm] <- sample(1:n)
}
WTPS.p[j] <- WTPSp(data.proj, group.label, group.label0, n, n.i, l, perm.WTPS = perm.WTPS, nrTRP)
if(permutationTRP == FALSE){
anova.p[j] <- 1-pf(anova.statistic.quick(data.proj, group.label), l-1, n-l)
ATS.simple.p[j] <- ATS.simple(data.proj, group.label)[2]
}else{
anova.s <- anova.statistic.quick(data.proj, group.label)
ATS.simple.s <- ATS.simple(data.proj, group.label)[1]
anova.perm <- numeric(nrTRP)
ATS.simple.perm <- numeric(nrTRP)
for(i.perm in 1:nrTRP){
anova.perm[i.perm] <- anova.statistic.quick(data.proj, sample(group.label))
ATS.simple.perm[i.perm] <- ATS.simple(data.proj, sample(group.label))[1]
}
anova.p[j] <- mean(anova.perm >= anova.s)
ATS.simple.p[j] <- mean(ATS.simple.perm >= ATS.simple.s)
}
}
all.data.proj <- ik.data.proj
iik <- 0
for (ik in k) {
iik <- iik + 1
pvalues.anova[iik] <- min(ik * anova.p[1:ik][order(anova.p[1:ik])] / 1:ik)
pvalues.ATS.simple[iik] <- min(ik * ATS.simple.p[1:ik][order(ATS.simple.p[1:ik])] / 1:ik)
pvalues.WTPS[iik] <- min(ik * WTPS.p[1:ik][order(WTPS.p[1:ik])] / 1:ik)
}
}
if(parallel.method == "parallel.method1"){
rs <- foreach(j = 1:ik, .combine = rbind) %dopar%
{
ik.data.proj[, j] <- t(x) %*% z[j,]
data.proj <- ik.data.proj[, j]
perm.WTPS <- matrix(0, nrow = n, ncol = nrTRP)
for(i.perm in 1:nrTRP){
perm.WTPS[, i.perm] <- sample(1:n)
}
WTPS.p <- WTPSp(data.proj, group.label, group.label0, n, n.i, l, perm.WTPS = perm.WTPS, nrTRP)
if(permutationTRP == FALSE){
c(1-pf(anova.statistic.quick(data.proj, group.label), l-1, n-l),
ATS.simple(data.proj, group.label)[2], WTPS.p, data.proj)
}else{
anova.s <- anova.statistic.quick(data.proj, group.label)
ATS.simple.s <- ATS.simple(data.proj, group.label)[1]
anova.perm <- numeric(nrTRP)
ATS.simple.perm <- numeric(nrTRP)
for(i.perm in 1:nrTRP){
anova.perm[i.perm] <- anova.statistic.quick(data.proj, sample(group.label))
ATS.simple.perm[i.perm] <- ATS.simple(data.proj, sample(group.label))[1]
}
c(mean(anova.perm >= anova.s),
mean(ATS.simple.perm >= ATS.simple.s), WTPS.p, data.proj)
}
}
if (ik == 1) rs <- matrix(rs, nrow = 1, ncol = n+3)
if (ik == 1) {
all.data.proj <- matrix(rs[, 4:(n+3)], ncol = 1)
} else {
all.data.proj <- t(rs[, 4:(n+3)])
}
iik <- 0
for (ik in k) {
iik <- iik + 1
pvalues.anova[iik] <- min(ik * (rs[1:ik, 1])[order(rs[1:ik, 1])] / 1:ik)
pvalues.ATS.simple[iik] <- min(ik * (rs[1:ik, 2])[order(rs[1:ik, 2])] / 1:ik)
pvalues.WTPS[iik] <- min(ik * (rs[1:ik, 3])[order(rs[1:ik, 3])] / 1:ik)
}
}
} else {
if(parallel.method == "parallel.method0"){
anova.p <- numeric(ik)
ATS.simple.p <- numeric(ik)
WTPS.p <- numeric(ik)
for(j in 1:ik){
bm.p <- cumsum(rnorm(p, mean = 0, sd = 1))/sqrt(p)
bm.p <- bm.p/modulo(bm.p)
ik.data.proj[, j] <- t(x) %*% as.matrix(bm.p)
data.proj <- ik.data.proj[, j]
perm.WTPS <- matrix(0, nrow = n, ncol = nrTRP)
for(i.perm in 1:nrTRP){
perm.WTPS[, i.perm] <- sample(1:n)
}
WTPS.p[j] <- WTPSp(data.proj, group.label, group.label0, n, n.i, l, perm.WTPS = perm.WTPS, nrTRP)
if(permutationTRP == FALSE){
anova.p[j] <- 1-pf(anova.statistic.quick(data.proj, group.label), l-1, n-l)
ATS.simple.p[j] <- ATS.simple(data.proj, group.label)[2]
}else{
anova.s <- anova.statistic.quick(data.proj, group.label)
ATS.simple.s <- ATS.simple(data.proj, group.label)[1]
anova.perm <- numeric(nrTRP)
ATS.simple.perm <- numeric(nrTRP)
for(i.perm in 1:nrTRP){
anova.perm[i.perm] <- anova.statistic.quick(data.proj, sample(group.label))
ATS.simple.perm[i.perm] <- ATS.simple(data.proj, sample(group.label))[1]
}
anova.p[j] <- mean(anova.perm >= anova.s)
ATS.simple.p[j] <- mean(ATS.simple.perm >= ATS.simple.s)
}
}
all.data.proj <- ik.data.proj
iik <- 0
for (ik in k) {
iik <- iik + 1
pvalues.anova[iik] <- min(ik * anova.p[1:ik][order(anova.p[1:ik])] / 1:ik)
pvalues.ATS.simple[iik] <- min(ik * ATS.simple.p[1:ik][order(ATS.simple.p[1:ik])] / 1:ik)
pvalues.WTPS[iik] <- min(ik * WTPS.p[1:ik][order(WTPS.p[1:ik])] / 1:ik)
}
}
if(parallel.method == "parallel.method1"){
rs <- foreach(j = 1:ik, .combine = rbind) %dopar%
{
bm.p <- cumsum(rnorm(p, mean = 0, sd = 1))/sqrt(p)
bm.p <- bm.p/modulo(bm.p)
ik.data.proj[, j] <- t(x) %*% as.matrix(bm.p)
data.proj <- ik.data.proj[, j]
perm.WTPS <- matrix(0, nrow = n, ncol = nrTRP)
for(i.perm in 1:nrTRP){
perm.WTPS[, i.perm] <- sample(1:n)
}
WTPS.p <- WTPSp(data.proj, group.label, group.label0, n, n.i, l, perm.WTPS = perm.WTPS, nrTRP)
if(permutationTRP == FALSE){
c(1-pf(anova.statistic.quick(data.proj, group.label), l-1, n-l),
ATS.simple(data.proj, group.label)[2], WTPS.p, data.proj)
}else{
anova.s <- anova.statistic.quick(data.proj, group.label)
ATS.simple.s <- ATS.simple(data.proj, group.label)[1]
anova.perm <- numeric(nrTRP)
ATS.simple.perm <- numeric(nrTRP)
for(i.perm in 1:nrTRP){
anova.perm[i.perm] <- anova.statistic.quick(data.proj, sample(group.label))
ATS.simple.perm[i.perm] <- ATS.simple(data.proj, sample(group.label))[1]
}
c(mean(anova.perm >= anova.s),
mean(ATS.simple.perm >= ATS.simple.s), WTPS.p, data.proj)
}
}
if (ik == 1) rs <- matrix(rs, nrow = 1, ncol = n+3)
if (ik == 1) {
all.data.proj <- matrix(rs[, 4:(n+3)], ncol = 1)
} else {
all.data.proj <- t(rs[, 4:(n+3)])
}
iik <- 0
for (ik in k) {
iik <- iik + 1
pvalues.anova[iik] <- min(ik * (rs[1:ik, 1])[order(rs[1:ik, 1])] / 1:ik)
pvalues.ATS.simple[iik] <- min(ik * (rs[1:ik, 2])[order(rs[1:ik, 2])] / 1:ik)
pvalues.WTPS[iik] <- min(ik * (rs[1:ik, 3])[order(rs[1:ik, 3])] / 1:ik)
}
}
}
}
resultTRP <- list(pvalues.anova = pvalues.anova, pvalues.ATS = pvalues.ATS.simple,
pvalues.WTPS = pvalues.WTPS, data.projections = all.data.proj,
k = k, projection = projection,
permutation = permutationTRP, B.TRP = nrTRP,
independent.projection.tests = independent.projection.tests)
}
if(any(c("FP", "CH", "CS", "L2b", "Fb", "Fmaxb", "TRP") %in% test)){
if(parallel.method == "parallel.method1"){
parallel::stopCluster(cl)
}
}
list.results = list()
for(ii in test){
if(ii == "FP"){
list.results[[which(test == ii)]] = resultFP
names(list.results)[which(test == ii)] = "FP"
}
if(ii == "CH"){
list.results[[which(test == ii)]] = resultCH
names(list.results)[which(test == ii)] = "CH"
}
if(ii == "CS"){
list.results[[which(test == ii)]] = resultCS
names(list.results)[which(test == ii)] = "CS"
}
if(ii == "L2N"){
list.results[[which(test == ii)]] = resultL2N
names(list.results)[which(test == ii)] = "L2N"
}
if(ii == "L2B"){
list.results[[which(test == ii)]] = resultL2B
names(list.results)[which(test == ii)] = "L2B"
}
if(ii == "L2b"){
list.results[[which(test == ii)]] = resultL2b
names(list.results)[which(test == ii)] = "L2b"
}
if(ii == "FN"){
list.results[[which(test == ii)]] = resultFN
names(list.results)[which(test == ii)] = "FN"
}
if(ii == "FB"){
list.results[[which(test == ii)]] = resultFB
names(list.results)[which(test == ii)] = "FB"
}
if(ii == "Fb"){
list.results[[which(test == ii)]] = resultFb
names(list.results)[which(test == ii)] = "Fb"
}
if(ii == "GPF"){
list.results[[which(test == ii)]] = resultGPF
names(list.results)[which(test == ii)] = "GPF"
}
if(ii == "Fmaxb"){
list.results[[which(test == ii)]] = resultFmaxb
names(list.results)[which(test == ii)] = "Fmaxb"
}
if(ii == "TRP"){
list.results[[which(test == ii)]] = resultTRP
names(list.results)[which(test == ii)] = "TRP"
}
}
list.results$data = x
list.results$group.label = group.label
list.results$parallel = parallel
list.results$nslaves = nslaves
class(list.results) = "fanovatests"
return(list.results)
}
Try the fdANOVA package in your browser
Any scripts or data that you put into this service are public.
fdANOVA documentation built on May 2, 2019, 2:43 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
fanova.tests = function(x = NULL, group.label, test = "ALL", params = NULL,
parallel = FALSE, nslaves = NULL){
if(any(is.na(group.label))){ stop("argument group.label can not contain NA values") }
if(test[1] == "ALL"){
test = c("FP", "CH", "CS", "L2N", "L2B", "L2b", "FN", "FB", "Fb", "GPF", "Fmaxb", "TRP")
}else{
for(i in 1:length(test)){
if(!(test[i] %in% c("FP", "CH", "CS", "L2N", "L2B", "L2b", "FN", "FB", "Fb", "GPF", "Fmaxb", "TRP"))){
stop("argument test must be the subvector of c('FP', 'CH', 'CS', 'L2N', 'L2B', 'L2b', 'FN', 'FB', 'Fb', 'GPF', 'Fmaxb', 'TRP')")
}
}
}
group.label0 = unique(group.label)
l = length(group.label0)
n.i = numeric(l)
for(i in 1:l) n.i[i] = sum(group.label == group.label0[i])
anova.statistic.quick = function(x, group.label){
n = length(c(x))
a = length(unique(group.label))
y = data.frame(cbind(x, group.label))
colnames(y) = c("data", "group.label")
means = (doBy::summaryBy(data ~ group.label, data = y))[,2]
mean.g = mean(x)
SSA = sum(table(group.label) * ((means - mean.g)^2))
SST = sum((x - mean.g)^2)
SSE = SST - SSA
return((n-a)*SSA/SSE/(a-1))
}
ATS.simple = function(x, group.label){
N = length(x); n = table(group.label); ng = length(as.vector(n))
means = matrix(tapply(x, group.label, mean), ncol = 1)
MM = diag(1, ng) - matrix(1, ncol = ng, nrow = ng)/ng
DM = diag(diag(MM))
SN = N * diag(as.vector((tapply(x, group.label, var) * (n - 1)/n^2)))
Delta = diag(as.vector(1/(n - 1)))
test.stat = N * t(means) %*% MM %*% means/sum(diag(MM %*% SN))
f1 = sum(diag(MM %*% SN))^2/sum(diag(MM %*% SN %*% MM %*% SN))
f0 = sum(diag(MM %*% SN))^2/sum(diag(DM %*% DM %*% SN %*% SN %*% Delta))
return(c(test.stat, 1 - pf(test.stat, f1, f0)))
}
WTPSp = function(x, group.label, group.label0, n, n.i, l, perm.WTPS, nrTRP){
means = matrix(0, nrow = l, ncol = 1); variances = numeric(l)
for(i in 1:l){
means[i,] = mean(x[group.label == group.label0[i]])
variances[i] = var(x[group.label == group.label0[i]])
}
V.N = n*diag(variances/n.i)
TT = diag(rep(1, l)) - (1/l)*matrix(1, ncol = l, nrow = l)
Q.N.T = n * t(means) %*% TT %*% MASS::ginv(TT %*% V.N %*% TT) %*% TT %*% means
A = t(rep(1/n.i[1], n.i[1]))
A1 = t(rep(1, n.i[1]))
for(i in 2:l){
A = magic::adiag(A, t(rep(1/n.i[i], n.i[i])))
A1 = magic::adiag(A1, t(rep(1, n.i[i])))
}
xperm = matrix(x[perm.WTPS], ncol = nrTRP)
meansP = A %*% xperm
varsP = (A1 %*% (xperm^2) - n.i * meansP^2)/(n.i * (n.i - 1))
statistics.p = sapply(1:nrTRP, function(arg){
TP = TT %*% MASS::ginv(TT %*% (n * diag(varsP[, arg])) %*% TT) %*% TT
statistics.p = diag(n * t(meansP[, arg]) %*% TP %*% meansP[, arg])
})
return(mean(statistics.p >= c(Q.N.T)))
}
if(any(c("FP", "CH", "CS", "L2b", "Fb", "Fmaxb", "TRP") %in% test)){
if(!parallel){
parallel.method = "parallel.method0"
}else{
if(!("doParallel" %in% rownames(installed.packages()))){
stop("Please install package 'doParallel'")
}
# require(foreach, quietly = TRUE)
requireNamespace("foreach", quietly = TRUE)
nlp = parallel::detectCores()
if(is.null(nslaves)){
if(nlp >= 2){
nslaves = nlp
parallel.method = "parallel.method1"
}else{
parallel.method = "parallel.method0"
}
}else{
if(nlp >= 2){
if(nslaves >= 2){
nslaves = nslaves
}else{
nslaves = nlp
}
parallel.method = "parallel.method1"
}else{
parallel.method = "parallel.method0"
}
}
if(parallel.method == "parallel.method1"){
cl = parallel::makePSOCKcluster(nslaves)
doParallel::registerDoParallel(cl)
}
}
}
if("FP" %in% test){
if(is.null(params$paramFP$basis)){ basis = "Fourier" }else{ basis = params$paramFP$basis }
if(!(basis %in% c("Fourier", "b-spline", "own"))){ stop("argument params$paramFP$basis must be one of the following: 'Fourier', 'b-spline', 'own'") }
if(is.null(params$paramFP$B.FP)){ nrFP = 1000 }else{ nrFP = params$paramFP$B.FP }
if(nrFP < 1){ stop("invalid number of permutations (params$paramFP$B.FP)") }
if(basis == "Fourier"){
if(!("fda" %in% rownames(installed.packages()))){
stop("Please install package 'fda'")
}
x = as.matrix(x); n = ncol(x); p = nrow(x)
if(is.null(params$paramFP$minK)){ minK = 3 }else{ minK = params$paramFP$minK }
if(minK < 1){ stop("invalid argument params$paramFP$minK") }
if(is.null(params$paramFP$maxK)){ if(p %% 2 == 1){ maxK = p - 2 }else{ maxK = p - 1 }
}else{ maxK = params$paramFP$maxK }
if(maxK < 1){ stop("invalid argument params$paramFP$maxK") }
if(p <= maxK){ if(p %% 2 == 1){ maxK = p - 2 }else{ maxK = p - 1 } }
if((maxK %% 2) == 0){ stop("the maximum number of basis functions (params$paramFP$maxK) is even") }
if(n != length(group.label)){
stop("number of observations (number of columns in x) and number of elements
in vector of group labels (group.label) must be the same")
}
if(is.null(params$paramFP$criterion)){ criterion = "BIC" }else{ criterion = params$paramFP$criterion }
if(!(criterion %in% c("BIC", "eBIC", "AIC", "AICc", "NO"))){ stop("argument params$paramFP$criterion must be one of the following: 'BIC', 'eBIC', 'AIC', 'AICc', 'NO'") }
if(criterion == "eBIC"){
if(is.null(params$paramFP$gamma.eBIC)){ gamma.eBIC = 0.5 }else{ gamma.eBIC = params$paramFP$gamma.eBIC }
if((gamma.eBIC < 0)|(gamma.eBIC > 1)){
stop("argument params$paramFP$gamma.eBIC must belong to [0,1]")
}
}
if(criterion != "NO"){
if(is.null(params$paramFP$commonK)){ commonK = "mode" }else{ commonK = params$paramFP$commonK }
if(!(commonK %in% c("mode", "min", "max", "mean"))){ stop("argument params$paramFP$commonK must be one of the following: 'mode', 'min', 'max', 'mean'") }
if((minK %% 2) == 0){ stop("the minimum number of basis functions (params$paramFP$minK) is even") }
if(is.null(params$paramFP$int)){
if(parallel.method == "parallel.method0"){
v = matrix(0, nrow = n, ncol = length(seq(minK, maxK, 2)))
for(i in seq(minK, maxK, 2)){
fbasis = fda::create.fourier.basis(c(0, p), i)
fdata = fda::smooth.basisPar(1:p, x, fbasis)
v[, (i+2-minK)/2] = switch(criterion, "BIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + i*log(p),
"eBIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + i*(log(p) + 2*gamma.eBIC*log(maxK)),
"AIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*i,
"AICc" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*i + (2 * i * (i + 1))/(n - i - 1))
}
}
if(parallel.method == "parallel.method1"){
v = foreach(i = seq(minK, maxK, 2), .combine = cbind) %dopar%
{
fbasis = fda::create.fourier.basis(c(0, p), i)
fdata = fda::smooth.basisPar(1:p, x, fbasis)
switch(criterion, "BIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + i*log(p),
"eBIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + i*(log(p) + 2*gamma.eBIC*log(maxK)),
"AIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*i,
"AICc" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*i + (2 * i * (i + 1))/(n - i - 1))
}
}
KK = numeric(n)
for(j in 1:n) KK[j] = 2*which(v[j,] == min(v[j,])) + minK - 2
if(commonK == "mode"){
temp = table(as.vector(KK))
K = min(as.numeric(names(temp)[temp == max(temp)]))
}
if(commonK == "min"){ K = min(KK) }
if(commonK == "max"){ K = max(KK) }
if(commonK == "mean"){ if(floor(mean(KK)) %% 2 == 1){ K = floor(mean(KK)) }else{ K = floor(mean(KK))-1 } }
}else{
if(length(params$paramFP$int) != 2){ stop("argument params$paramFP$int must be of length two") }
if(parallel.method == "parallel.method0"){
v = matrix(0, nrow = n, ncol = length(seq(minK, maxK, 2)))
for(i in seq(minK, maxK, 2)){
fbasis = fda::create.fourier.basis(rangeval = c(params$paramFP$int[1], params$paramFP$int[2]), nbasis = i)
fdata = fda::smooth.basisPar(seq(params$paramFP$int[1], params$paramFP$int[2], length = p), x, fbasis)
v[, (i+2-minK)/2] = switch(criterion, "BIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + i*log(p),
"eBIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + i*(log(p) + 2*gamma.eBIC*log(maxK)),
"AIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*i,
"AICc" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*i + (2 * i * (i + 1))/(n - i - 1))
}
}
if(parallel.method == "parallel.method1"){
v = foreach(i = seq(minK, maxK, 2), .combine = cbind) %dopar%
{
fbasis = fda::create.fourier.basis(rangeval = c(params$paramFP$int[1], params$paramFP$int[2]), nbasis = i)
fdata = fda::smooth.basisPar(seq(params$paramFP$int[1], params$paramFP$int[2], length = p), x, fbasis)
switch(criterion, "BIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + i*log(p),
"eBIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + i*(log(p) + 2*gamma.eBIC*log(maxK)),
"AIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*i,
"AICc" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*i + (2 * i * (i + 1))/(n - i - 1))
}
}
KK = numeric(n)
for(j in 1:n) KK[j] = 2*which(v[j,] == min(v[j,])) + minK - 2
if(commonK == "mode"){
temp = table(as.vector(KK))
K = min(as.numeric(names(temp)[temp == max(temp)]))
}
if(commonK == "min"){ K = min(KK) }
if(commonK == "max"){ K = max(KK) }
if(commonK == "mean"){ if(floor(mean(KK)) %% 2 == 1){ K = floor(mean(KK)) }else{ K = floor(mean(KK))-1 } }
}
}else{ K = maxK }
if(is.null(params$paramFP$int)){
fbasis = fda::create.fourier.basis(c(0, p), K)
lfdata = vector("list", l)
mfdata = matrix(0, nrow = K, ncol = n)
for(i in 1:l){
lfdata[[i]] = fda::Data2fd(1:p, x[, group.label == group.label0[i]], fbasis)$coefs
mfdata[, group.label == group.label0[i]] = lfdata[[i]]
}
}else{
if(length(params$paramFP$int) != 2){ stop("argument params$paramFP$int must be of length two") }
fbasis = fda::create.fourier.basis(c(params$paramFP$int[1], params$paramFP$int[2]), K)
lfdata = vector("list", l)
mfdata = matrix(0, nrow = K, ncol = n)
for(i in 1:l){
lfdata[[i]] = fda::Data2fd(seq(params$paramFP$int[1], params$paramFP$int[2], length = p), x[, group.label == group.label0[i]], fbasis)$coefs
mfdata[, group.label == group.label0[i]] = lfdata[[i]]
}
}
a = 0; b = 0; c = 0
for(i in 1:l){
a = a + sum(t(lfdata[[i]]) %*% lfdata[[i]])/n.i[i]
for(j in 1:l) b = b + sum(t(lfdata[[i]]) %*% (lfdata[[j]]))
c = c + sum((lfdata[[i]])^2)
}
b = b/n; statFP0 = (a-b)/(c-a); statFP = statFP0*(n-l)/(l-1)
if(parallel.method == "parallel.method0"){
FP.p = numeric(nrFP)
for(ii in 1:nrFP){
a = numeric(l)
Perm = sample(1:n)
for(j in 1:l){
if(j == 1){ l1 = Perm[1:n.i[1]] }else{ l1 = Perm[(sum(n.i[1:(j-1)])+1):sum(n.i[1:j])] }
a[j] = sum(t(mfdata[, l1]) %*% mfdata[, l1])/n.i[j]
}
FP.p[ii] = (sum(a)-b)/(c-sum(a))
}
}
if(parallel.method == "parallel.method1"){
FP.p = foreach(ii = 1:nrFP, .combine = "c") %dopar%
{
a = numeric(l)
Perm = sample(1:n)
for(j in 1:l){
if(j == 1){ l1 = Perm[1:n.i[1]] }else{ l1 = Perm[(sum(n.i[1:(j-1)])+1):sum(n.i[1:j])] }
a[j] = sum(t(mfdata[, l1]) %*% mfdata[, l1])/n.i[j]
}
(sum(a)-b)/(c-sum(a))
}
}
pvalueFP = mean(FP.p >= statFP0)
}
if(basis == "b-spline"){
if(!("fda" %in% rownames(installed.packages()))){
stop("Please install package 'fda'")
}
x = as.matrix(x); n = ncol(x); p = nrow(x)
if(is.null(params$paramFP$norder)){ norder = 4 }else{ norder = params$paramFP$norder }
if(norder < 1){ stop("invalid argument params$paramFP$norder") }
if(is.null(params$paramFP$minK)){ minK = norder }else{ minK = params$paramFP$minK }
if(minK < 1){ stop("invalid argument params$paramFP$minK") }
if(minK < norder){ minK = norder }
if(is.null(params$paramFP$maxK)){ maxK = p }else{ maxK = params$paramFP$maxK }
if(maxK < 1){ stop("invalid argument params$paramFP$maxK") }
if(p <= maxK){ maxK = p }
if(n != length(group.label)){
stop("number of observations (number of columns in x) and number of elements
in vector of group labels (group.label) must be the same")
}
if(is.null(params$paramFP$criterion)){ criterion = "BIC" }else{ criterion = params$paramFP$criterion }
if(!(criterion %in% c("BIC", "eBIC", "AIC", "AICc", "NO"))){ stop("argument params$paramFP$criterion must be one of the following: 'BIC', 'eBIC', 'AIC', 'AICc', 'NO'") }
if(criterion == "eBIC"){
if(is.null(params$paramFP$gamma.eBIC)){ gamma.eBIC = 0.5 }else{ gamma.eBIC = params$paramFP$gamma.eBIC }
if((gamma.eBIC < 0)|(gamma.eBIC > 1)){
stop("argument params$paramFP$gamma.eBIC must belong to [0,1]")
}
}
if(criterion != "NO"){
if(is.null(params$paramFP$commonK)){ commonK = "mode" }else{ commonK = params$paramFP$commonK }
if(!(commonK %in% c("mode", "min", "max", "mean"))){ stop("argument params$paramFP$commonK must be one of the following: 'mode', 'min', 'max', 'mean'") }
mmm = minK:maxK
if(is.null(params$paramFP$int)){
if(parallel.method == "parallel.method0"){
v = matrix(0, nrow = n, ncol = length(mmm))
for(i in 1:length(mmm)){
fbasis = fda::create.bspline.basis(c(0, p), mmm[i], norder = norder)
fdata = fda::smooth.basisPar(1:p, x, fbasis)
v[, i] = switch(criterion, "BIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + mmm[i]*log(p),
"eBIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + mmm[i]*(log(p) + 2*gamma.eBIC*log(maxK)),
"AIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*mmm[i],
"AICc" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*mmm[i] + (2 * mmm[i] * (mmm[i] + 1))/(n - mmm[i] - 1))
}
}
if(parallel.method == "parallel.method1"){
v = foreach(i = 1:length(mmm), .combine = cbind) %dopar%
{
fbasis = fda::create.bspline.basis(c(0, p), mmm[i], norder = norder)
fdata = fda::smooth.basisPar(1:p, x, fbasis)
switch(criterion, "BIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + mmm[i]*log(p),
"eBIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + mmm[i]*(log(p) + 2*gamma.eBIC*log(maxK)),
"AIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*mmm[i],
"AICc" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*mmm[i] + (2 * mmm[i] * (mmm[i] + 1))/(n - mmm[i] - 1))
}
}
KK = numeric(n)
for(j in 1:n) KK[j] = mmm[which(v[j,] == min(v[j,]))]
if(commonK == "mode"){
temp = table(as.vector(KK))
K = min(as.numeric(names(temp)[temp == max(temp)]))
}
if(commonK == "min"){ K = min(KK) }
if(commonK == "max"){ K = max(KK) }
if(commonK == "mean"){ K = floor(mean(KK)) }
}else{
if(length(params$paramFP$int) != 2){ stop("argument params$paramFP$int must be of length two") }
if(parallel.method == "parallel.method0"){
v = matrix(0, nrow = n, ncol = length(mmm))
for(i in 1:length(mmm)){
fbasis = fda::create.bspline.basis(rangeval = c(params$paramFP$int[1], params$paramFP$int[2]), nbasis = mmm[i], norder = norder)
fdata = fda::smooth.basisPar(seq(params$paramFP$int[1], params$paramFP$int[2], length = p), x, fbasis)
v[, i] = switch(criterion, "BIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + mmm[i]*log(p),
"eBIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + mmm[i]*(log(p) + 2*gamma.eBIC*log(maxK)),
"AIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*mmm[i],
"AICc" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*mmm[i] + (2 * mmm[i] * (mmm[i] + 1))/(n - mmm[i] - 1))
}
}
if(parallel.method == "parallel.method1"){
v = foreach(i = 1:length(mmm), .combine = cbind) %dopar%
{
fbasis = fda::create.bspline.basis(rangeval = c(params$paramFP$int[1], params$paramFP$int[2]), nbasis = mmm[i], norder = norder)
fdata = fda::smooth.basisPar(seq(params$paramFP$int[1], params$paramFP$int[2], length = p), x, fbasis)
switch(criterion, "BIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + mmm[i]*log(p),
"eBIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + mmm[i]*(log(p) + 2*gamma.eBIC*log(maxK)),
"AIC" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*mmm[i],
"AICc" = p*log(p*(1-fdata$df/p)^2*fdata$gcv/p) + 2*mmm[i] + (2 * mmm[i] * (mmm[i] + 1))/(n - mmm[i] - 1))
}
}
KK = numeric(n)
for(j in 1:n) KK[j] = mmm[which(v[j,] == min(v[j,]))]
if(commonK == "mode"){
temp = table(as.vector(KK))
K = min(as.numeric(names(temp)[temp == max(temp)]))
}
if(commonK == "min"){ K = min(KK) }
if(commonK == "max"){ K = max(KK) }
if(commonK == "mean"){ K = floor(mean(KK)) }
}
}else{ K = maxK }
if(is.null(params$paramFP$int)){
fbasis = fda::create.bspline.basis(c(0, p), K, norder = norder)
lfdata = vector("list", l)
mfdata = matrix(0, nrow = K, ncol = n)
for(i in 1:l){
lfdata[[i]] = fda::Data2fd(1:p, x[, group.label == group.label0[i]], fbasis)$coefs
mfdata[, group.label == group.label0[i]] = lfdata[[i]]
}
}else{
if(length(params$paramFP$int) != 2){ stop("argument params$paramFP$int must be of length two") }
fbasis = fda::create.bspline.basis(c(params$paramFP$int[1], params$paramFP$int[2]), K, norder = norder)
lfdata = vector("list", l)
mfdata = matrix(0, nrow = K, ncol = n)
for(i in 1:l){
lfdata[[i]] = fda::Data2fd(seq(params$paramFP$int[1], params$paramFP$int[2], length = p), x[, group.label == group.label0[i]], fbasis)$coefs
mfdata[, group.label == group.label0[i]] = lfdata[[i]]
}
}
bspline.cross.prod.matrix = fda::inprod(fbasis, fbasis)
a = 0; b = 0; c = 0
for(i in 1:l){
a = a + sum(t(lfdata[[i]]) %*% bspline.cross.prod.matrix %*% lfdata[[i]])/n.i[i]
for(j in 1:l) b = b + sum(t(lfdata[[i]]) %*% bspline.cross.prod.matrix %*% (lfdata[[j]]))
c = c + sum(diag(t(lfdata[[i]]) %*% bspline.cross.prod.matrix %*% lfdata[[i]]))
}
b = b/n; statFP0 = (a-b)/(c-a); statFP = statFP0*(n-l)/(l-1)
if(parallel.method == "parallel.method0"){
FP.p = numeric(nrFP)
for(ii in 1:nrFP){
a = numeric(l)
Perm = sample(1:n)
for(j in 1:l){
if(j == 1){ l1 = Perm[1:n.i[1]] }else{ l1 = Perm[(sum(n.i[1:(j-1)])+1):sum(n.i[1:j])] }
a[j] = sum(t(mfdata[, l1]) %*% bspline.cross.prod.matrix %*% mfdata[, l1])/n.i[j]
}
FP.p[ii] = (sum(a)-b)/(c-sum(a))
}
}
if(parallel.method == "parallel.method1"){
FP.p = foreach(ii = 1:nrFP, .combine = "c") %dopar%
{
a = numeric(l)
Perm = sample(1:n)
for(j in 1:l){
if(j == 1){ l1 = Perm[1:n.i[1]] }else{ l1 = Perm[(sum(n.i[1:(j-1)])+1):sum(n.i[1:j])] }
a[j] = sum(t(mfdata[, l1]) %*% bspline.cross.prod.matrix %*% mfdata[, l1])/n.i[j]
}
(sum(a)-b)/(c-sum(a))
}
}
pvalueFP = mean(FP.p >= statFP0)
}
if(basis == "own"){
own.basis = params$paramFP$own.basis
mfdata = own.basis
n = ncol(mfdata); K = nrow(mfdata)
if(n != length(group.label)){ stop("the number of observations is not equal to the number of elements in vector of labels") }
lfdata = vector("list", l)
for(i in 1:l){ lfdata[[i]] = mfdata[, group.label == group.label0[i]] }
own.cross.prod.mat = params$paramFP$own.cross.prod.mat
if(K != nrow(own.cross.prod.mat)){ stop("invalid argument params$paramFP$own.cross.prod.mat") }
if(K != ncol(own.cross.prod.mat)){ stop("invalid argument params$paramFP$own.cross.prod.mat") }
a = 0; b = 0; c = 0
for(i in 1:l){
a = a + sum(t(lfdata[[i]]) %*% own.cross.prod.mat %*% lfdata[[i]])/n.i[i]
for(j in 1:l) b = b + sum(t(lfdata[[i]]) %*% own.cross.prod.mat %*% (lfdata[[j]]))
c = c + sum(diag(t(lfdata[[i]]) %*% own.cross.prod.mat %*% lfdata[[i]]))
}
b = b/n; statFP0 = (a-b)/(c-a); statFP = statFP0*(n-l)/(l-1)
if(parallel.method == "parallel.method0"){
FP.p = numeric(nrFP)
for(ii in 1:nrFP){
a = numeric(l)
Perm = sample(1:n)
for(j in 1:l){
if(j == 1){ l1 = Perm[1:n.i[1]] }else{ l1 = Perm[(sum(n.i[1:(j-1)])+1):sum(n.i[1:j])] }
a[j] = sum(t(mfdata[, l1]) %*% own.cross.prod.mat %*% mfdata[, l1])/n.i[j]
}
FP.p[ii] = (sum(a)-b)/(c-sum(a))
}
}
if(parallel.method == "parallel.method1"){
FP.p = foreach(ii = 1:nrFP, .combine = "c") %dopar%
{
a = numeric(l)
Perm = sample(1:n)
for(j in 1:l){
if(j == 1){ l1 = Perm[1:n.i[1]] }else{ l1 = Perm[(sum(n.i[1:(j-1)])+1):sum(n.i[1:j])] }
a[j] = sum(t(mfdata[, l1]) %*% own.cross.prod.mat %*% mfdata[, l1])/n.i[j]
}
(sum(a)-b)/(c-sum(a))
}
}
pvalueFP = mean(FP.p >= statFP0)
}
if(basis == "own"){
resultFP = list(statFP = statFP, pvalueFP = pvalueFP, B.FP = nrFP, basis = basis, own.basis = own.basis, own.cross.prod.mat = own.cross.prod.mat)
}else{
if(criterion == "NO"){
if(is.null(params$paramFP$int)){
if(basis == "Fourier"){
resultFP = list(statFP = statFP, pvalueFP = pvalueFP, B.FP = nrFP, basis = basis, criterion = criterion, K = K, maxK = maxK)
}
if(basis == "b-spline"){
resultFP = list(statFP = statFP, pvalueFP = pvalueFP, B.FP = nrFP, basis = basis, criterion = criterion, K = K, maxK = maxK, norder = norder)
}
}else{
if(basis == "Fourier"){
resultFP = list(statFP = statFP, pvalueFP = pvalueFP, int = params$paramFP$int, B.FP = nrFP, basis = basis, criterion = criterion, K = K, maxK = maxK)
}
if(basis == "b-spline"){
resultFP = list(statFP = statFP, pvalueFP = pvalueFP, int = params$paramFP$int, B.FP = nrFP, basis = basis, criterion = criterion, K = K, maxK = maxK, norder = norder)
}
}
}else{
if(is.null(params$paramFP$int)){
if(criterion == "eBIC"){
if(basis == "Fourier"){
resultFP = list(statFP = statFP, pvalueFP = pvalueFP, B.FP = nrFP, basis = basis, criterion = criterion, commonK = commonK, K = K, minK = minK, maxK = maxK, gamma.eBIC = gamma.eBIC)
}
if(basis == "b-spline"){
resultFP = list(statFP = statFP, pvalueFP = pvalueFP, B.FP = nrFP, basis = basis, criterion = criterion, commonK = commonK, K = K, minK = minK, maxK = maxK, norder = norder, gamma.eBIC = gamma.eBIC)
}
}else{
if(basis == "Fourier"){
resultFP = list(statFP = statFP, pvalueFP = pvalueFP, B.FP = nrFP, basis = basis, criterion = criterion, commonK = commonK, K = K, minK = minK, maxK = maxK)
}
if(basis == "b-spline"){
resultFP = list(statFP = statFP, pvalueFP = pvalueFP, B.FP = nrFP, basis = basis, criterion = criterion, commonK = commonK, K = K, minK = minK, maxK = maxK, norder = norder)
}
}
}else{
if(criterion == "eBIC"){
if(basis == "Fourier"){
resultFP = list(statFP = statFP, pvalueFP = pvalueFP, int = params$paramFP$int, B.FP = nrFP, basis = basis, criterion = criterion, commonK = commonK, K = K, minK = minK, maxK = maxK, gamma.eBIC = gamma.eBIC)
}
if(basis == "b-spline"){
resultFP = list(statFP = statFP, pvalueFP = pvalueFP, int = params$paramFP$int, B.FP = nrFP, basis = basis, criterion = criterion, commonK = commonK, K = K, minK = minK, maxK = maxK, norder = norder, gamma.eBIC = gamma.eBIC)
}
}else{
if(basis == "Fourier"){
resultFP = list(statFP = statFP, pvalueFP = pvalueFP, int = params$paramFP$int, B.FP = nrFP, basis = basis, criterion = criterion, commonK = commonK, K = K, minK = minK, maxK = maxK)
}
if(basis == "b-spline"){
resultFP = list(statFP = statFP, pvalueFP = pvalueFP, int = params$paramFP$int, B.FP = nrFP, basis = basis, criterion = criterion, commonK = commonK, K = K, minK = minK, maxK = maxK, norder = norder)
}
}
}
}
}
}
if(any(c("CH", "CS") %in% test)){
if(!("MASS" %in% rownames(installed.packages()))){
stop("Please install package 'MASS'")
}
x = as.matrix(x); n = ncol(x)
if(n != length(group.label)){
stop("number of observations (number of columns in x) and number of elements
in vector of group labels (group.label) must be the same")
}
statCHCS = 0
for(i in 1:(l-1)){
for(j in (i+1):l){
statCHCS = statCHCS + n.i[i]*sum((rowMeans(x[, group.label == group.label0[i]]) -
rowMeans(x[, group.label == group.label0[j]]))^2)
}
}
if("CH" %in% test){
if(is.null(params$paramCH)){ nrCH = 10000 }else{ nrCH = params$paramCH }
if(nrCH < 1){ stop("invalid number of discretized artificial trajectories for each process Z_i(t) (params$paramCH)") }
cov.fun = ((n-1)/(n-l))*var(t(x)); Z = vector("list", l)
for(i in 1:l) Z[[i]] = MASS::mvrnorm(nrCH, numeric(nrow(cov.fun)), cov.fun)
if(parallel.method == "parallel.method0"){
V = numeric(nrCH)
for(m in 1:nrCH){
for(i in 1:(l-1)){
for(j in (i+1):l){
V[m] = V[m] + sum((Z[[i]][m,]-sqrt(n.i[i]/n.i[j])*Z[[j]][m,])^2)
}
}
}
}
if(parallel.method == "parallel.method1"){
V = foreach(m = 1:nrCH, .combine = "c") %dopar%
{
V = 0
for(i in 1:(l-1)){
for(j in (i+1):l){
V = V + sum((Z[[i]][m,]-sqrt(n.i[i]/n.i[j])*Z[[j]][m,])^2)
}
}
V
}
}
resultCH = list(statCH = statCHCS, pvalueCH = mean(V > statCHCS), paramCH = nrCH)
}
if("CS" %in% test){
if(is.null(params$paramCS)){ nrCS = 10000 }else{ nrCS = params$paramCS }
if(nrCS < 1){ stop("invalid number of discretized artificial trajectories for each process Z_i(t) (params$paramCS)") }
Z = vector("list", l)
for(i in 1:l){
cov.fun = var(t(x[, group.label == group.label0[i]]))
Z[[i]] = MASS::mvrnorm(nrCS, numeric(nrow(cov.fun)), cov.fun)
}
if(parallel.method == "parallel.method0"){
V = numeric(nrCS)
for(m in 1:nrCS){
for(i in 1:(l-1)){
for(j in (i+1):l){
V[m] = V[m] + sum((Z[[i]][m,]-sqrt(n.i[i]/n.i[j])*Z[[j]][m,])^2)
}
}
}
}
if(parallel.method == "parallel.method1"){
V = foreach(m = 1:nrCS, .combine = "c") %dopar%
{
V = 0
for(i in 1:(l-1)){
for(j in (i+1):l){
V = V + sum((Z[[i]][m,]-sqrt(n.i[i]/n.i[j])*Z[[j]][m,])^2)
}
}
V
}
}
resultCS = list(statCS = statCHCS, pvalueCS = mean(V > statCHCS), paramCS = nrCS)
}
}
if(any(c("L2N", "L2B", "L2b", "FN", "FB", "Fb", "GPF", "Fmaxb") %in% test)){
x = as.matrix(x); n = ncol(x); p = nrow(x)
if(n != length(group.label)){
stop("number of observations (number of columns in x) and number of elements
in vector of group labels (group.label) must be the same")
}
mu0 = rowMeans(x)
vmu = matrix(0, nrow = l, ncol = p)
z = matrix(0, nrow = n, ncol = p)
SSR = 0; SSE = 0
for(i in 1:l){
xi = x[, group.label == group.label0[i]]
mui = rowMeans(xi); vmu[i,] = mui
zi = t(xi) - as.matrix(rep(1, n.i[i])) %*% mui
if(i==1){ z[1:n.i[i],] = zi }else{ z[(cumsum(n.i)[i-1]+1):cumsum(n.i)[i],] = zi }
SSR = SSR + n.i[i]*(mui - mu0)^2
SSE = SSE + colSums(zi^2)
}
if(n > p){ gamma.z = t(z) %*% z/(n-l) }else{ gamma.z = z %*% t(z)/(n-l) }
A = sum(diag(gamma.z)); B = sum(diag(gamma.z %*% gamma.z))
A2N = A^2; B2N = B
A2B = (n-l)*(n-l+1)/(n-l-1)/(n-l+2)*(A^2-2*B/(n-l+1))
B2B = (n-l)^2/(n-l-1)/(n-l+2)*(B-A^2/(n-l))
if(any(c("L2N", "L2B", "L2b") %in% test)){
statL2 = sum(SSR)
if("L2N" %in% test){
betaL2N = B2N/A; kappaL2N = A2N/B2N
pvalueL2N = 1-pchisq(statL2/betaL2N, (l-1)*kappaL2N)
resultL2N = list(statL2 = statL2, pvalueL2N = pvalueL2N, betaL2N = betaL2N, dL2N = (l-1)*kappaL2N)
}
if("L2B" %in% test){
betaL2B = B2B/A; kappaL2B = A2B/B2B
pvalueL2B = 1-pchisq(statL2/betaL2B, (l-1)*kappaL2B)
resultL2B = list(statL2 = statL2, pvalueL2B = pvalueL2B, betaL2B = betaL2B, dL2B = (l-1)*kappaL2B)
}
if("L2b" %in% test){
if(is.null(params$paramL2b)){ nrL2b = 10000 }else{ nrL2b = params$paramL2b }
if(nrL2b < 1){ stop("invalid number of bootstrap replicates (params$paramL2b)") }
if(parallel.method == "parallel.method0"){
statL2boot = numeric(nrL2b)
for(ii in 1:nrL2b){
vmuboot = matrix(0, nrow = l, ncol = p)
for(i in 1:l){
xi = x[, group.label == group.label0[i]]
xiboot = xi[, floor(runif(n.i[i])*(n.i[i]-1)) + 1]
vmuboot[i,] = rowMeans(xiboot)-vmu[i,]
}
mu0boot = n.i %*% vmuboot/n
SSRboot = 0
for(i in 1:l) SSRboot = SSRboot + n.i[i]*(vmuboot[i,]-mu0boot)^2
statL2boot[ii] = sum(SSRboot)
}
}
if(parallel.method == "parallel.method1"){
statL2boot = foreach(ii = 1:nrL2b, .combine = "c") %dopar%
{
vmuboot = matrix(0, nrow = l, ncol = p)
for(i in 1:l){
xi = x[, group.label == group.label0[i]]
xiboot = xi[, floor(runif(n.i[i])*(n.i[i]-1)) + 1]
vmuboot[i,] = rowMeans(xiboot)-vmu[i,]
}
mu0boot = n.i %*% vmuboot/n
SSRboot = 0
for(i in 1:l) SSRboot = SSRboot + n.i[i]*(vmuboot[i,]-mu0boot)^2
sum(SSRboot)
}
}
pvalueL2b = mean(statL2boot >= statL2)
resultL2b = list(statL2 = statL2, pvalueL2b = pvalueL2b, paramL2b = nrL2b)
}
}
if(any(c("FN", "FB", "Fb") %in% test)){
statF = sum(SSR)/sum(SSE)*(n-l)/(l-1)
if("FN" %in% test){
kappaFN = A2N/B2N
pvalueFN = 1-pf(statF, (l-1)*kappaFN, (n-l)*kappaFN)
resultFN = list(statF = statF, pvalueFN = pvalueFN, d1FN = (l-1)*kappaFN, d2FN = (n-l)*kappaFN)
}
if("FB" %in% test){
kappaFB = A2B/B2B
pvalueFB = 1-pf(statF, (l-1)*kappaFB, (n-l)*kappaFB)
resultFB = list(statF = statF, pvalueFB = pvalueFB, d1FB = (l-1)*kappaFB, d2FB = (n-l)*kappaFB)
}
if("Fb" %in% test){
if(is.null(params$paramFb)){ nrFb = 10000 }else{ nrFb = params$paramFb }
if(nrFb < 1){ stop("invalid number of bootstrap replicates (params$paramFb)") }
if(parallel.method == "parallel.method0"){
statFboot = numeric(nrFb)
for(ii in 1:nrFb){
vmuboot = matrix(0, nrow = l, ncol = p)
zboot = matrix(0, nrow = n, ncol = p)
for(i in 1:l){
xi = x[, group.label == group.label0[i]]
xiboot = xi[, floor(runif(n.i[i])*(n.i[i]-1)) + 1]
muiboot = rowMeans(xiboot)
ziboot = t(xiboot) - as.matrix(rep(1, n.i[i])) %*% muiboot
if(i==1){ zboot[1:n.i[i],] = ziboot }else{ zboot[(cumsum(n.i)[i-1]+1):cumsum(n.i)[i],] = ziboot }
vmuboot[i,] = rowMeans(xiboot)-vmu[i,]
}
mu0boot = n.i %*% vmuboot/n
SSRboot = 0
for(i in 1:l) SSRboot = SSRboot + n.i[i]*(vmuboot[i,]-mu0boot)^2
if(n > p){ gammaboot = t(zboot) %*% zboot/(n-l) }else{ gammaboot = zboot %*% t(zboot)/(n-l) }
statFboot[ii] = sum(SSRboot)/sum(diag(gammaboot))/(l-1)
}
}
if(parallel.method == "parallel.method1"){
statFboot = foreach(ii = 1:nrFb, .combine = "c") %dopar%
{
vmuboot = matrix(0, nrow = l, ncol = p)
zboot = matrix(0, nrow = n, ncol = p)
for(i in 1:l){
xi = x[, group.label == group.label0[i]]
xiboot = xi[, floor(runif(n.i[i])*(n.i[i]-1)) + 1]
muiboot = rowMeans(xiboot)
ziboot = t(xiboot) - as.matrix(rep(1, n.i[i])) %*% muiboot
if(i==1){ zboot[1:n.i[i],] = ziboot }else{ zboot[(cumsum(n.i)[i-1]+1):cumsum(n.i)[i],] = ziboot }
vmuboot[i,] = rowMeans(xiboot)-vmu[i,]
}
mu0boot = n.i %*% vmuboot/n
SSRboot = 0
for(i in 1:l) SSRboot = SSRboot + n.i[i]*(vmuboot[i,]-mu0boot)^2
if(n > p){ gammaboot = t(zboot) %*% zboot/(n-l) }else{ gammaboot = zboot %*% t(zboot)/(n-l) }
sum(SSRboot)/sum(diag(gammaboot))/(l-1)
}
}
pvalueFb = mean(statFboot >= statF)
resultFb = list(statF = statF, pvalueFb = pvalueFb, paramFb = nrFb)
}
}
if("GPF" %in% test){
statGPF = mean(SSR/SSE*(n-l)/(l-1))
z = z/(as.matrix(rep(1, n)) %*% sqrt(colSums(z^2)))
if(n >= p){ z = t(z) %*% z }else{ z = z %*% t(z) }
betaGPF = (sum(diag(z %*% z))/p^2/(l-1))/((n-l)/(n-l-2))
dGPF = ((n-l)/(n-l-2))^2/(sum(diag(z %*% z))/p^2/(l-1))
pvalueGPF = 1-pchisq(statGPF/betaGPF,dGPF)
resultGPF = list(statGPF = statGPF, pvalueGPF = pvalueGPF, betaGPF = betaGPF, dGPF = dGPF)
}
if("Fmaxb" %in% test){
if(is.null(params$paramFmaxb)){ nrFmaxb = 10000 }else{ nrFmaxb = params$paramFmaxb }
if(nrFmaxb < 1){ stop("invalid number of bootstrap replicates (params$paramFmaxb)") }
statFmax = max(SSR/SSE*(n-l)/(l-1))
xs = matrix(0, nrow = n, ncol = p)
for(i in 1:l){ xs[group.label == group.label0[i],] = t(x[, group.label == group.label0[i]]) - as.matrix(rep(1, n.i[i])) %*% vmu[i,] }
if(parallel.method == "parallel.method0"){
statFmaxboot = numeric(nrFmaxb)
for(ii in 1:nrFmaxb){
vmuboot = matrix(0, nrow = l, ncol = p)
zboot = matrix(0, nrow = n, ncol = p)
xboot = xs[sample(1:n, replace = TRUE),]
mu0boot = colMeans(xboot)
SSRboot = 0
for(i in 1:l){
xiboot = xboot[group.label == group.label0[i],]
muiboot = colMeans(xiboot)
vmuboot[i,] = muiboot
ziboot = xiboot - as.matrix(rep(1, n.i[i])) %*% muiboot
if(i==1){ zboot[1:n.i[i],] = ziboot }else{ zboot[(cumsum(n.i)[i-1]+1):cumsum(n.i)[i],] = ziboot }
SSRboot = SSRboot + n.i[i]*(muiboot - mu0boot)^2
}
SSEboot = diag(t(zboot) %*% zboot)
statFmaxboot[ii] = max(SSRboot/SSEboot*(n-l)/(l-1))
}
}
if(parallel.method == "parallel.method1"){
statFmaxboot = foreach(ii = 1:nrFmaxb, .combine = "c") %dopar%
{
vmuboot = matrix(0, nrow = l, ncol = p)
zboot = matrix(0, nrow = n, ncol = p)
xboot = xs[sample(1:n, replace = TRUE),]
mu0boot = colMeans(xboot)
SSRboot = 0
for(i in 1:l){
xiboot = xboot[group.label == group.label0[i],]
muiboot = colMeans(xiboot)
vmuboot[i,] = muiboot
ziboot = xiboot - as.matrix(rep(1, n.i[i])) %*% muiboot
if(i==1){ zboot[1:n.i[i],] = ziboot }else{ zboot[(cumsum(n.i)[i-1]+1):cumsum(n.i)[i],] = ziboot }
SSRboot = SSRboot + n.i[i]*(muiboot - mu0boot)^2
}
SSEboot = diag(t(zboot) %*% zboot)
max(SSRboot/SSEboot*(n-l)/(l-1))
}
}
pvalueFmaxb = mean(statFmaxboot >= statFmax)
resultFmaxb = list(statFmax = statFmax, pvalueFmaxb = pvalueFmaxb, paramFmaxb = nrFmaxb)
}
}
if ("TRP" %in% test) {
if (!("doBy" %in% rownames(installed.packages()))) {
stop("Please install package 'doBY'")
}
if (!("MASS" %in% rownames(installed.packages()))) {
stop("Please install package 'MASS'")
}
if (!("magic" %in% rownames(installed.packages()))) {
stop("Please install package 'magic'")
}
x <- as.matrix(x); n <- ncol(x); p <- nrow(x)
if (n != length(group.label)) {
stop("number of observations (number of columns in x) and number of elements
in vector of group labels (group.label) must be the same")
}
if (is.null(params$paramTRP$projection)) {
projection <- "GAUSS"
}else{
projection <- params$paramTRP$projection
}
if (!(projection %in% c("GAUSS", "BM"))) {
stop("argument params$paramTRP$projection must be one of the following: 'GAUSS', 'BM'")
}
if (is.null(params$paramTRP$permutation)) {
permutationTRP <- FALSE
}else{
permutationTRP <- params$paramTRP$permutation
}
if (! is.logical(permutationTRP)) {
stop("argument permutation is not logical (params$paramTRP$permutation)")
}
if (is.null(params$paramTRP$independent.projection.tests)) {
independent.projection.tests <- TRUE
}else{
independent.projection.tests <- params$paramTRP$independent.projection.tests
}
if (! is.logical(independent.projection.tests)) {
stop("argument independent.projection.tests is not logical
(params$paramTRP$independent.projection.tests)")
}
if (is.null(params$paramTRP$B.TRP)) {
nrTRP <- 10000
}else{
nrTRP <- params$paramTRP$B.TRP
}
if (nrTRP < 1) {
stop("invalid number of permutations (params$paramTRP$B.TRP)")
}
if (is.null(params$paramTRP$k)) {
k <- 30
}else{
k <- params$paramTRP$k
}
if (any(k < 1)) {
stop("invalid number of projections (k)")
}
pvalues.anova <- numeric(length(k))
pvalues.ATS.simple <- numeric(length(k))
pvalues.WTPS <- numeric(length(k))
modulo <- function(z){ sqrt(sum(z^2)) }
if (independent.projection.tests) {
all.data.proj <- list()
iik <- 0
for(ik in k){
ik.data.proj <- matrix(0, nrow = n, ncol = ik)
if(projection == "GAUSS"){
z <- matrix(rnorm(p * ik), nrow = ik, ncol = p)
modu <- apply(z, 1, modulo)
z <- z/modu
if(parallel.method == "parallel.method0"){
anova.p <- numeric(ik)
ATS.simple.p <- numeric(ik)
WTPS.p <- numeric(ik)
for(j in 1:ik){
ik.data.proj[, j] <- t(x) %*% z[j,]
data.proj <- ik.data.proj[, j]
perm.WTPS <- matrix(0, nrow = n, ncol = nrTRP)
for(i.perm in 1:nrTRP){
perm.WTPS[, i.perm] <- sample(1:n)
}
WTPS.p[j] <- WTPSp(data.proj, group.label, group.label0, n, n.i, l, perm.WTPS = perm.WTPS, nrTRP)
if(permutationTRP == FALSE){
anova.p[j] <- 1-pf(anova.statistic.quick(data.proj, group.label), l-1, n-l)
ATS.simple.p[j] <- ATS.simple(data.proj, group.label)[2]
}else{
anova.s <- anova.statistic.quick(data.proj, group.label)
ATS.simple.s <- ATS.simple(data.proj, group.label)[1]
anova.perm <- numeric(nrTRP)
ATS.simple.perm <- numeric(nrTRP)
for(i.perm in 1:nrTRP){
anova.perm[i.perm] <- anova.statistic.quick(data.proj, sample(group.label))
ATS.simple.perm[i.perm] <- ATS.simple(data.proj, sample(group.label))[1]
}
anova.p[j] <- mean(anova.perm >= anova.s)
ATS.simple.p[j] <- mean(ATS.simple.perm >= ATS.simple.s)
}
}
iik <- iik + 1
all.data.proj[[iik]] <- ik.data.proj
pvalues.anova[iik] <- min(ik*anova.p[order(anova.p)]/1:ik)
pvalues.ATS.simple[iik] <- min(ik*ATS.simple.p[order(ATS.simple.p)]/1:ik)
pvalues.WTPS[iik] <- min(ik*WTPS.p[order(WTPS.p)]/1:ik)
}
if(parallel.method == "parallel.method1"){
rs <- foreach(j = 1:ik, .combine = rbind) %dopar%
{
ik.data.proj[, j] <- t(x) %*% z[j,]
data.proj <- ik.data.proj[, j]
perm.WTPS <- matrix(0, nrow = n, ncol = nrTRP)
for(i.perm in 1:nrTRP){
perm.WTPS[, i.perm] <- sample(1:n)
}
WTPS.p <- WTPSp(data.proj, group.label, group.label0, n, n.i, l, perm.WTPS = perm.WTPS, nrTRP)
if(permutationTRP == FALSE){
c(1-pf(anova.statistic.quick(data.proj, group.label), l-1, n-l),
ATS.simple(data.proj, group.label)[2], WTPS.p, data.proj)
}else{
anova.s <- anova.statistic.quick(data.proj, group.label)
ATS.simple.s <- ATS.simple(data.proj, group.label)[1]
anova.perm <- numeric(nrTRP)
ATS.simple.perm <- numeric(nrTRP)
for(i.perm in 1:nrTRP){
anova.perm[i.perm] <- anova.statistic.quick(data.proj, sample(group.label))
ATS.simple.perm[i.perm] <- ATS.simple(data.proj, sample(group.label))[1]
}
c(mean(anova.perm >= anova.s),
mean(ATS.simple.perm >= ATS.simple.s), WTPS.p, data.proj)
}
}
if(ik == 1) rs <- matrix(rs, nrow = 1, ncol = n+3)
iik <- iik + 1
if(ik == 1){
all.data.proj[[iik]] <- matrix(rs[, 4:(n+3)], ncol = 1)
}else{
all.data.proj[[iik]] <- t(rs[, 4:(n+3)])
}
pvalues.anova[iik] <- min(ik*(rs[, 1])[order(rs[, 1])]/1:ik)
pvalues.ATS.simple[iik] <- min(ik*(rs[, 2])[order(rs[, 2])]/1:ik)
pvalues.WTPS[iik] <- min(ik*(rs[, 3])[order(rs[, 3])]/1:ik)
}
}else{
if(parallel.method == "parallel.method0"){
anova.p <- numeric(ik)
ATS.simple.p <- numeric(ik)
WTPS.p <- numeric(ik)
for(j in 1:ik){
bm.p <- cumsum(rnorm(p, mean = 0, sd = 1))/sqrt(p)
bm.p <- bm.p/modulo(bm.p)
ik.data.proj[, j] <- t(x) %*% as.matrix(bm.p)
data.proj <- ik.data.proj[, j]
perm.WTPS <- matrix(0, nrow = n, ncol = nrTRP)
for(i.perm in 1:nrTRP){
perm.WTPS[, i.perm] <- sample(1:n)
}
WTPS.p[j] <- WTPSp(data.proj, group.label, group.label0, n, n.i, l, perm.WTPS = perm.WTPS, nrTRP)
if(permutationTRP == FALSE){
anova.p[j] <- 1-pf(anova.statistic.quick(data.proj, group.label), l-1, n-l)
ATS.simple.p[j] <- ATS.simple(data.proj, group.label)[2]
}else{
anova.s <- anova.statistic.quick(data.proj, group.label)
ATS.simple.s <- ATS.simple(data.proj, group.label)[1]
anova.perm <- numeric(nrTRP)
ATS.simple.perm <- numeric(nrTRP)
for(i.perm in 1:nrTRP){
anova.perm[i.perm] <- anova.statistic.quick(data.proj, sample(group.label))
ATS.simple.perm[i.perm] <- ATS.simple(data.proj, sample(group.label))[1]
}
anova.p[j] <- mean(anova.perm >= anova.s)
ATS.simple.p[j] <- mean(ATS.simple.perm >= ATS.simple.s)
}
}
iik <- iik + 1
all.data.proj[[iik]] <- ik.data.proj
pvalues.anova[iik] <- min(ik*anova.p[order(anova.p)]/1:ik)
pvalues.ATS.simple[iik] <- min(ik*ATS.simple.p[order(ATS.simple.p)]/1:ik)
pvalues.WTPS[iik] <- min(ik*WTPS.p[order(WTPS.p)]/1:ik)
}
if(parallel.method == "parallel.method1"){
rs <- foreach(j = 1:ik, .combine = rbind) %dopar%
{
bm.p <- cumsum(rnorm(p, mean = 0, sd = 1))/sqrt(p)
bm.p <- bm.p/modulo(bm.p)
ik.data.proj[, j] <- t(x) %*% as.matrix(bm.p)
data.proj <- ik.data.proj[, j]
perm.WTPS <- matrix(0, nrow = n, ncol = nrTRP)
for(i.perm in 1:nrTRP){
perm.WTPS[, i.perm] <- sample(1:n)
}
WTPS.p <- WTPSp(data.proj, group.label, group.label0, n, n.i, l, perm.WTPS = perm.WTPS, nrTRP)
if(permutationTRP == FALSE){
c(1-pf(anova.statistic.quick(data.proj, group.label), l-1, n-l),
ATS.simple(data.proj, group.label)[2], WTPS.p, data.proj)
}else{
anova.s <- anova.statistic.quick(data.proj, group.label)
ATS.simple.s <- ATS.simple(data.proj, group.label)[1]
anova.perm <- numeric(nrTRP)
ATS.simple.perm <- numeric(nrTRP)
for(i.perm in 1:nrTRP){
anova.perm[i.perm] <- anova.statistic.quick(data.proj, sample(group.label))
ATS.simple.perm[i.perm] <- ATS.simple(data.proj, sample(group.label))[1]
}
c(mean(anova.perm >= anova.s),
mean(ATS.simple.perm >= ATS.simple.s), WTPS.p, data.proj)
}
}
if(ik == 1) rs <- matrix(rs, nrow = 1, ncol = n+3)
iik <- iik + 1
if(ik == 1){
all.data.proj[[iik]] <- matrix(rs[, 4:(n+3)], ncol = 1)
}else{
all.data.proj[[iik]] <- t(rs[, 4:(n+3)])
}
pvalues.anova[iik] <- min(ik*(rs[, 1])[order(rs[, 1])]/1:ik)
pvalues.ATS.simple[iik] <- min(ik*(rs[, 2])[order(rs[, 2])]/1:ik)
pvalues.WTPS[iik] <- min(ik*(rs[, 3])[order(rs[, 3])]/1:ik)
}
}
}
}else{
ik <- max(k)
ik.data.proj <- matrix(0, nrow = n, ncol = ik)
if (projection == "GAUSS") {
z <- matrix(rnorm(p * ik), nrow = ik, ncol = p)
modu <- apply(z, 1, modulo)
z <- z/modu
if(parallel.method == "parallel.method0"){
anova.p <- numeric(ik)
ATS.simple.p <- numeric(ik)
WTPS.p <- numeric(ik)
for(j in 1:ik){
ik.data.proj[, j] <- t(x) %*% z[j,]
data.proj <- ik.data.proj[, j]
perm.WTPS <- matrix(0, nrow = n, ncol = nrTRP)
for(i.perm in 1:nrTRP){
perm.WTPS[, i.perm] <- sample(1:n)
}
WTPS.p[j] <- WTPSp(data.proj, group.label, group.label0, n, n.i, l, perm.WTPS = perm.WTPS, nrTRP)
if(permutationTRP == FALSE){
anova.p[j] <- 1-pf(anova.statistic.quick(data.proj, group.label), l-1, n-l)
ATS.simple.p[j] <- ATS.simple(data.proj, group.label)[2]
}else{
anova.s <- anova.statistic.quick(data.proj, group.label)
ATS.simple.s <- ATS.simple(data.proj, group.label)[1]
anova.perm <- numeric(nrTRP)
ATS.simple.perm <- numeric(nrTRP)
for(i.perm in 1:nrTRP){
anova.perm[i.perm] <- anova.statistic.quick(data.proj, sample(group.label))
ATS.simple.perm[i.perm] <- ATS.simple(data.proj, sample(group.label))[1]
}
anova.p[j] <- mean(anova.perm >= anova.s)
ATS.simple.p[j] <- mean(ATS.simple.perm >= ATS.simple.s)
}
}
all.data.proj <- ik.data.proj
iik <- 0
for (ik in k) {
iik <- iik + 1
pvalues.anova[iik] <- min(ik * anova.p[1:ik][order(anova.p[1:ik])] / 1:ik)
pvalues.ATS.simple[iik] <- min(ik * ATS.simple.p[1:ik][order(ATS.simple.p[1:ik])] / 1:ik)
pvalues.WTPS[iik] <- min(ik * WTPS.p[1:ik][order(WTPS.p[1:ik])] / 1:ik)
}
}
if(parallel.method == "parallel.method1"){
rs <- foreach(j = 1:ik, .combine = rbind) %dopar%
{
ik.data.proj[, j] <- t(x) %*% z[j,]
data.proj <- ik.data.proj[, j]
perm.WTPS <- matrix(0, nrow = n, ncol = nrTRP)
for(i.perm in 1:nrTRP){
perm.WTPS[, i.perm] <- sample(1:n)
}
WTPS.p <- WTPSp(data.proj, group.label, group.label0, n, n.i, l, perm.WTPS = perm.WTPS, nrTRP)
if(permutationTRP == FALSE){
c(1-pf(anova.statistic.quick(data.proj, group.label), l-1, n-l),
ATS.simple(data.proj, group.label)[2], WTPS.p, data.proj)
}else{
anova.s <- anova.statistic.quick(data.proj, group.label)
ATS.simple.s <- ATS.simple(data.proj, group.label)[1]
anova.perm <- numeric(nrTRP)
ATS.simple.perm <- numeric(nrTRP)
for(i.perm in 1:nrTRP){
anova.perm[i.perm] <- anova.statistic.quick(data.proj, sample(group.label))
ATS.simple.perm[i.perm] <- ATS.simple(data.proj, sample(group.label))[1]
}
c(mean(anova.perm >= anova.s),
mean(ATS.simple.perm >= ATS.simple.s), WTPS.p, data.proj)
}
}
if (ik == 1) rs <- matrix(rs, nrow = 1, ncol = n+3)
if (ik == 1) {
all.data.proj <- matrix(rs[, 4:(n+3)], ncol = 1)
} else {
all.data.proj <- t(rs[, 4:(n+3)])
}
iik <- 0
for (ik in k) {
iik <- iik + 1
pvalues.anova[iik] <- min(ik * (rs[1:ik, 1])[order(rs[1:ik, 1])] / 1:ik)
pvalues.ATS.simple[iik] <- min(ik * (rs[1:ik, 2])[order(rs[1:ik, 2])] / 1:ik)
pvalues.WTPS[iik] <- min(ik * (rs[1:ik, 3])[order(rs[1:ik, 3])] / 1:ik)
}
}
} else {
if(parallel.method == "parallel.method0"){
anova.p <- numeric(ik)
ATS.simple.p <- numeric(ik)
WTPS.p <- numeric(ik)
for(j in 1:ik){
bm.p <- cumsum(rnorm(p, mean = 0, sd = 1))/sqrt(p)
bm.p <- bm.p/modulo(bm.p)
ik.data.proj[, j] <- t(x) %*% as.matrix(bm.p)
data.proj <- ik.data.proj[, j]
perm.WTPS <- matrix(0, nrow = n, ncol = nrTRP)
for(i.perm in 1:nrTRP){
perm.WTPS[, i.perm] <- sample(1:n)
}
WTPS.p[j] <- WTPSp(data.proj, group.label, group.label0, n, n.i, l, perm.WTPS = perm.WTPS, nrTRP)
if(permutationTRP == FALSE){
anova.p[j] <- 1-pf(anova.statistic.quick(data.proj, group.label), l-1, n-l)
ATS.simple.p[j] <- ATS.simple(data.proj, group.label)[2]
}else{
anova.s <- anova.statistic.quick(data.proj, group.label)
ATS.simple.s <- ATS.simple(data.proj, group.label)[1]
anova.perm <- numeric(nrTRP)
ATS.simple.perm <- numeric(nrTRP)
for(i.perm in 1:nrTRP){
anova.perm[i.perm] <- anova.statistic.quick(data.proj, sample(group.label))
ATS.simple.perm[i.perm] <- ATS.simple(data.proj, sample(group.label))[1]
}
anova.p[j] <- mean(anova.perm >= anova.s)
ATS.simple.p[j] <- mean(ATS.simple.perm >= ATS.simple.s)
}
}
all.data.proj <- ik.data.proj
iik <- 0
for (ik in k) {
iik <- iik + 1
pvalues.anova[iik] <- min(ik * anova.p[1:ik][order(anova.p[1:ik])] / 1:ik)
pvalues.ATS.simple[iik] <- min(ik * ATS.simple.p[1:ik][order(ATS.simple.p[1:ik])] / 1:ik)
pvalues.WTPS[iik] <- min(ik * WTPS.p[1:ik][order(WTPS.p[1:ik])] / 1:ik)
}
}
if(parallel.method == "parallel.method1"){
rs <- foreach(j = 1:ik, .combine = rbind) %dopar%
{
bm.p <- cumsum(rnorm(p, mean = 0, sd = 1))/sqrt(p)
bm.p <- bm.p/modulo(bm.p)
ik.data.proj[, j] <- t(x) %*% as.matrix(bm.p)
data.proj <- ik.data.proj[, j]
perm.WTPS <- matrix(0, nrow = n, ncol = nrTRP)
for(i.perm in 1:nrTRP){
perm.WTPS[, i.perm] <- sample(1:n)
}
WTPS.p <- WTPSp(data.proj, group.label, group.label0, n, n.i, l, perm.WTPS = perm.WTPS, nrTRP)
if(permutationTRP == FALSE){
c(1-pf(anova.statistic.quick(data.proj, group.label), l-1, n-l),
ATS.simple(data.proj, group.label)[2], WTPS.p, data.proj)
}else{
anova.s <- anova.statistic.quick(data.proj, group.label)
ATS.simple.s <- ATS.simple(data.proj, group.label)[1]
anova.perm <- numeric(nrTRP)
ATS.simple.perm <- numeric(nrTRP)
for(i.perm in 1:nrTRP){
anova.perm[i.perm] <- anova.statistic.quick(data.proj, sample(group.label))
ATS.simple.perm[i.perm] <- ATS.simple(data.proj, sample(group.label))[1]
}
c(mean(anova.perm >= anova.s),
mean(ATS.simple.perm >= ATS.simple.s), WTPS.p, data.proj)
}
}
if (ik == 1) rs <- matrix(rs, nrow = 1, ncol = n+3)
if (ik == 1) {
all.data.proj <- matrix(rs[, 4:(n+3)], ncol = 1)
} else {
all.data.proj <- t(rs[, 4:(n+3)])
}
iik <- 0
for (ik in k) {
iik <- iik + 1
pvalues.anova[iik] <- min(ik * (rs[1:ik, 1])[order(rs[1:ik, 1])] / 1:ik)
pvalues.ATS.simple[iik] <- min(ik * (rs[1:ik, 2])[order(rs[1:ik, 2])] / 1:ik)
pvalues.WTPS[iik] <- min(ik * (rs[1:ik, 3])[order(rs[1:ik, 3])] / 1:ik)
}
}
}
}
resultTRP <- list(pvalues.anova = pvalues.anova, pvalues.ATS = pvalues.ATS.simple,
pvalues.WTPS = pvalues.WTPS, data.projections = all.data.proj,
k = k, projection = projection,
permutation = permutationTRP, B.TRP = nrTRP,
independent.projection.tests = independent.projection.tests)
}
if(any(c("FP", "CH", "CS", "L2b", "Fb", "Fmaxb", "TRP") %in% test)){
if(parallel.method == "parallel.method1"){
parallel::stopCluster(cl)
}
}
list.results = list()
for(ii in test){
if(ii == "FP"){
list.results[[which(test == ii)]] = resultFP
names(list.results)[which(test == ii)] = "FP"
}
if(ii == "CH"){
list.results[[which(test == ii)]] = resultCH
names(list.results)[which(test == ii)] = "CH"
}
if(ii == "CS"){
list.results[[which(test == ii)]] = resultCS
names(list.results)[which(test == ii)] = "CS"
}
if(ii == "L2N"){
list.results[[which(test == ii)]] = resultL2N
names(list.results)[which(test == ii)] = "L2N"
}
if(ii == "L2B"){
list.results[[which(test == ii)]] = resultL2B
names(list.results)[which(test == ii)] = "L2B"
}
if(ii == "L2b"){
list.results[[which(test == ii)]] = resultL2b
names(list.results)[which(test == ii)] = "L2b"
}
if(ii == "FN"){
list.results[[which(test == ii)]] = resultFN
names(list.results)[which(test == ii)] = "FN"
}
if(ii == "FB"){
list.results[[which(test == ii)]] = resultFB
names(list.results)[which(test == ii)] = "FB"
}
if(ii == "Fb"){
list.results[[which(test == ii)]] = resultFb
names(list.results)[which(test == ii)] = "Fb"
}
if(ii == "GPF"){
list.results[[which(test == ii)]] = resultGPF
names(list.results)[which(test == ii)] = "GPF"
}
if(ii == "Fmaxb"){
list.results[[which(test == ii)]] = resultFmaxb
names(list.results)[which(test == ii)] = "Fmaxb"
}
if(ii == "TRP"){
list.results[[which(test == ii)]] = resultTRP
names(list.results)[which(test == ii)] = "TRP"
}
}
list.results$data = x
list.results$group.label = group.label
list.results$parallel = parallel
list.results$nslaves = nslaves
class(list.results) = "fanovatests"
return(list.results)
}
Try the fdANOVA package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.