# NOT RUN {
# 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') )
# }
Run the code above in your browser using DataLab