Heter.test {HeterFunctionalData}

Heteroscedastic test for functional data

Description

This function conducts the test of main effect of treatment, interaction effect of treatment and time, main effect of time, and simple effect of treament based on original observations (see Wang and Akritas 2010a) or ranks (see Wang and Akritas 2010b).

Usage

Heter.test(
  data,
  a,
  b,
  mn,
  h = 0.45,
  method = "rank",
  Ca = cbind(as.vector(rep(1, a - 1)), -diag(a - 1))
)

Arguments

data

The data in long format (see the example in function dataformat_wide_to_long( )).

a

The number of treatments.

b

The number of time points or repeated measurements per subject.

mn

The vector of sample sizes in treatments.

h

The h value used in the estimators in Proposition 3.3 of Wang and Akritas (2010a) and Theorem 3.2 of Wang, Higgins, and Blasi (2010). The value of h should be in (0, 0.5) for the test of no main treatment effect or contract effect of treatments. For other tests, h should be in (0, 1). Recommend to use the default value h=0.45 as given in the function. Note: If multiple values are provided to h as a vector, then the calculation will be carried out for each h value, which results in multiple p-values in the returned result for palpha, pbeta, pgamma, pphi.

method

Specifying method='rank' to use rank test. For all other values, the test based on original data will be used.

Ca

Contrast matrix for the contrast effect of the treatments.

Value

a matrix of p-values for the test of no main treatment effect (palpha), test of no main time effect (pbeta), test of no interaction effect between treatment and time (pgamma), and test of no simple effect of treatment (pphi). For each h value, the p-values are in the same column.

References

Haiyan Wang and Michael Akritas (2010a). Inference from heteroscedastic functional data, Journal of Nonparametric Statistics. 22:2, 149-168. DOI: 10.1080/10485250903171621

Haiyan Wang and Michael Akritas (2010b). Rank test for heteroscedastic functional data. Journal of Multivariate Analysis. 101: 1791-1805. https://doi.org/10.1016/j.jmva.2010.03.012

Haiyan Wang, James Higgins, and Dale Blasi (2010). Distribution-Free Tests For No Effect Of Treatment In Heteroscedastic Functional Data Under Both Weak And Long Range Dependence. Statistics and Probability Letters. 80: 390-402. Doi:10.1016/j.spl.2009.11.016

Examples

# Generate a data set that contains data from 3 treatments,
# with 3 subjects in treatment 1, 3 subjects in treatment 2,
# and 4 subjects in treatment 3. Each subject contains m=50
# repeated observations from Poisson distribution. For the 1st treatment,
# the mean vector of the repeated observations from the same subject is
# equal to mu1 plus a random effect vector generated by NorRanGen( ).
# The m is the number of repeated measurements per subject.
 f1<-function(m, mu1, raneff) {
     currentmu=mu1+raneff;
     currentmu[abs(currentmu)<1e-2]=1e-2;
     rpois(m, abs(currentmu))}
 f2<-function(m, mu2, raneff) {
     currentmu=mu2+raneff;
     currentmu[abs(currentmu)<1e-2]=1e-2;
     rpois(m, abs(currentmu))}
f3<- function(m, mu3, raneff){
     currentmu=mu3+raneff;
     currentmu[abs(currentmu)<1e-2]=1e-2;
     rpois(m, abs(currentmu))}
# The a is the number of treatments. The mn stores the number of subjects in treatments.
a=3;  mn=c(3, 3, 4); mu1=3; mu2=3; mu3=3; m=50
# Generate the time effects via random effects with AR(1) structure.
raneff=NorRanGen(m)
# Generate data and store in wide format.
datawide=numeric()
now=0
for (i in 1:a){
   fi=function(x1, x2) f1(m,x1, x2)*(i==1)+f2(m,x1, x2)*(i==2)+f3(m, x1, x2)*(i==3)
   mu=mu1*(i==1)+mu2*(i==2)+mu3*(i==3)
   for (k in 1:mn[i]){
       now=now+1
       datawide<-rbind(datawide, c(k, i, fi(mu, raneff) )   	)
       colnames(datawide)=c("sub", "trt", paste("time", seq(m), sep=""))
       #this is a typical way to store data in practice
   }
} #end of j

# Note:There are different time effects since values in raneff vector are different
dat=dataformat_wide_to_long(datawide)   #dat is in long format
# Define the h value used in Proposition 3.3 of Wang and Akritas (2010a)
h=c(0.45, 0.49)
  #Note: The resulting palpha, pbeta, pgamma, pphi each contains
  #      two p-values, one corresponds to each h value
  #      (see Proposition 3.3 of Wang and Akritas (2010a))
# test based on original data.
(org= Heter.test(dat, a, m, mn, h, method='original') )
#test based on ranks
(rankt= Heter.test(dat, a, m, mn, h, method='rank') )