R/NPtest.R

In eRm: Extended Rasch Modeling

NPtest<-function(obj, n=NULL, method="T1", ...){

   dots<-as.list(substitute(list(...)))[-1]
   nn<-names(dots)
   for (i in seq(along=dots)) assign(nn[i],dots[[i]])

   if(!exists("burn_in", inherits = FALSE)) burn_in <- 256
   if(!("step" %in% nn)) step<-32
   if(!exists("seed", inherits = FALSE)) seed<-0
   if(is.null(n)) n <- 500

   if(is.matrix(obj) || is.data.frame(obj)){ # input is datamatrix -  RaschSampler object is generated
      if (!all(obj %in% 0:1)) stop("Data matrix must be binary, NAs not allowed")
      itscor<-colSums(obj) # rh 2011-03-03
      itcol<-(itscor==0|itscor==nrow(obj))
      if (any(itcol)){
        cat("The following columns in the data show complete 0/full responses: \n")
        cat((1:ncol(obj))[itcol],sep=", ")
        cat("\n")
        stop("NPtest using these items is meaningless. Delete them first!")
      }
      obj<-rsampler(obj,rsctrl(burn_in=burn_in, n_eff=n, step=step, seed=seed))
#browser()

   } else if(!("RSmpl" %in% class(obj))){
        stop("Input object must be data matrix/data frame or output from RaschSampler")
   }

   if(exists("RSinfo", inherits = FALSE)) if(get("RSinfo")) summary(obj)

  switch(method,
    "T1"    = T1(obj, ...),
    "T1l"   = T1l(obj, ...),
    "T1m"   = T1m(obj, ...),
    "Tmd"   = Tmd(obj, ...),
    "T2"    = T2(obj, ...),
    "T2m"   = T2m(obj, ...),
    "T4"    = T4(obj, ...),
#    "T7"    = T7(obj, ...),
#    "T7a"   = T7a(obj, ...),
    "T10"   = T10(obj, ...),
    "T11"   = T11(obj, ...),
    "Tpbis" = Tpbis(obj, ...),
    "MLoef" = MLoef.x(obj, ...),
	"Q3h" = Q3h(obj, ...),
	"Q3l" = Q3l(obj, ...)
  )
}

MLoef.x<-function(rsobj, splitcr=NULL, ...){
     # user function
     MLexact<-function(X,splitcr){
       rmod<-RM(X)
       LR<-MLoef(rmod,splitcr)$LR
       LR
     }
     #if(!exists("splitcr", inherits = FALSE)) splitcr="median"
     if(is.null(splitcr)) splitcr="median"
     res <- rstats(rsextrobj(rsobj, 2), MLexact, splitcr)

     rmod<-RM(rsextrmat(rsobj,1))                     # MLoef for original data
     MLres<-MLoef(rmod,splitcr)
     class(MLres)<-c(class(MLres),"MLx")              # for printing without blank line
     res1<-MLres$LR

     n_eff<-rsobj$n_eff                         # number of simulated matrices
     res<-unlist(res)
     prop<-sum((res[1:n_eff]>=res1)/n_eff)

     result<-list(MLres=MLres, n_eff=n_eff, prop=prop, MLoefvec=res) # MLobj
     class(result)<-"MLobj"
     result
}

Tpbis <- function(rsobj, idxt=NULL, idxs=NULL, ...){ # fixed 2013-08-09
  Tpbis.stat <- function(x){
    rb <- rowSums(x[, idxs, drop = FALSE])     # all raw scores
    t  <- x[, idxt]                            # dichotomous item
    r  <- tapply(rb, t, sum, simplify = FALSE) # raw scores by item; simplify = FALSE to be on the safe side
    n1 <- sum(t)                               # n_1 = sum of raw scores with t == 1
    n0 <- sum(1 - t)                           # n_0 = sum of raw scores with t == 0
    return(n0 * r[[2L]][1L] - n1*r[[1L]][1L])  # n_0 * sum(r_1) - n_1 * sum(r_0)
  }

  if(is.null(idxs)) stop("No item(s) for subscale  specified (use idxs!)")
  if(is.null(idxt)) stop("No test item for testing against subscale specified (use idx!)")
  li1 <- length(idxt)
  li2 <- length(idxs)
  k   <- rsobj$k
  if(li1 > 1L ||li2 >= k || (li1 + li2) > k || any(idxt %in% idxs) || any(c(idxt,idxs) > k)){
    stop("Subscale and/or test item incorrectly specified.")
  }

  n_eff <- rsobj$n_eff                   # number of simulated matrices
  n_tot <- rsobj$n_tot                   # number of simulated matrices

  res     <- rstats(rsobj, Tpbis.stat)              # calculates statistic for each matrix
  corrvec <- do.call(cbind, lapply(res, as.vector)) # converts result list to matrix

  prop <- sum(corrvec[2L:(n_tot)] <= corrvec[1L]) / n_eff   # T(A_s) >= T(A_0)

  # Tpbisobj
  result <- list("n_eff"    = n_eff,
                 "prop"     = prop,
                 "idxt"     = idxt,
                 "idxs"     = idxs,
                 "Tpbisvec" = corrvec)
  class(result)<-"Tpbisobj"
  return(result)
}

Tmd<-function(rsobj, idx1=NULL, idx2=NULL, ...){
     Tmd.stat<-function(x){
        r1<-rowSums(x[,idx1, drop=FALSE])
        r2<-rowSums(x[,idx2, drop=FALSE])
        corr<-cor(r1,r2)
        corr
     }

     if(is.null(idx1))
         stop("No item(s) for subscale 1 specified (use idx1!)")
     if(is.null(idx2))
         stop("No item(s) for subscale 2 specified (use idx2!)")
     li1<-length(idx1)
     li2<-length(idx2)
     k<-rsobj$k
     if(li1>=k ||li2>=k || li1+li2>k || any(idx1 %in% idx2))
         stop("Subscale(s) incorrectly specified.")

     n_eff<-rsobj$n_eff                         # number of simulated matrices
     n_tot<-rsobj$n_tot                         # number of simulated matrices

     res<-rstats(rsobj,Tmd.stat)               # calculates statistic for each matrix
     corrvec<-do.call(cbind, lapply(res,as.vector)) # converts result list to matrix

     prop<-sum(corrvec[2:(n_tot)]<=corrvec[1])/n_eff

     result<-list(n_eff=n_eff, prop=prop, idx1=idx1, idx2=idx2, Tmdvec=corrvec)   # Tmdobj
     class(result)<-"Tmdobj"
     result
}


T1m<-function(rsobj, ...){
     T1mstat<-function(x){      # calculates statistic T1m
        unlist(lapply(1:(k-1),function(i) lapply((i+1):k, function(j) sum(x[,i]==x[,j]))))
     }
     n_eff<-rsobj$n_eff                         # number of simulated matrices
     n_tot<-rsobj$n_tot                         # number of simulated matrices
     k<-rsobj$k                                 # number of columns of matrices

     res<-rstats(rsobj,T1mstat)                  # calculates statistic for each matrix

     res<-do.call(cbind, lapply(res,as.vector)) # converts result list to matrix
     T1mvec<-apply(res, 1, function(x) sum(x[2:(n_tot)]<=x[1])/n_eff)
     T1mmat<-matrix(,k,k)
     T1mmat[lower.tri(T1mmat)] <- T1mvec           # lower triangular matrix of p-values
     result<-list(n_eff=n_eff, prop=T1mvec, T1mmat=T1mmat) # T1mobj
     class(result)<-"T1mobj"
     result
}

T1<-function(rsobj, ...){
     T1stat<-function(x){      # calculates statistic T1
        unlist(lapply(1:(k-1),function(i) lapply((i+1):k, function(j) sum(x[,i]==x[,j]))))
     }
     n_eff<-rsobj$n_eff                         # number of simulated matrices
     n_tot<-rsobj$n_tot                         # number of simulated matrices
     k<-rsobj$k                                 # number of columns of matrices

     res<-rstats(rsobj,T1stat)                  # calculates statistic for each matrix

     res<-do.call(cbind, lapply(res,as.vector)) # converts result list to matrix
     T1vec<-apply(res, 1, function(x) sum(x[2:(n_tot)]>=x[1])/n_eff)
     T1mat<-matrix(,k,k)
     T1mat[lower.tri(T1mat)] <- T1vec           # lower triangular matrix of p-values
     result<-list(n_eff=n_eff, prop=T1vec, T1mat=T1mat) # T1obj
     class(result)<-"T1obj"
     result
}

T1l<-function(rsobj, ...){
     T1lstat<-function(x){      # calculates statistic T1
        unlist(lapply(1:(k-1),function(i) lapply((i+1):k, function(j) sum(x[,i] & x[,j]))))
     }
     n_eff<-rsobj$n_eff                         # number of simulated matrices
     n_tot<-rsobj$n_tot                         # number of simulated matrices
     k<-rsobj$k                                 # number of columns of matrices

     res<-rstats(rsobj,T1lstat)                  # calculates statistic for each matrix

     res<-do.call(cbind, lapply(res,as.vector)) # converts result list to matrix
     T1lvec<-apply(res, 1, function(x) sum(x[2:(n_tot)]>=x[1])/n_eff)
     T1lmat<-matrix(,k,k)
     T1lmat[lower.tri(T1lmat)] <- T1lvec           # lower triangular matrix of p-values
     result<-list(n_eff=n_eff, prop=T1lvec, T1lmat=T1lmat) # T1obj
     class(result)<-"T1lobj"
     result
}
T2<-function(rsobj,idx=NULL,stat="var", ...){

     T2.Var.stat<-function(x){       # calculates statistic T2
        var(rowSums(x[,idx, drop=FALSE]))
     }
     T2.MAD1.stat<-function(x){       # calculates statistic T2
        y<-rowSums(x[,idx, drop=FALSE])           # mean absolute deviation
        mean(abs(y-mean(y)))
     }
     T2.MAD2.stat<-function(x){       # calculates statistic T2
        mad(rowSums(x[,idx, drop=FALSE]),constant=1) # unscaled median absolute deviation
     }
     T2.Range.stat<-function(x){     # calculates statistic T2
        diff(range(rowSums(x[,idx, drop=FALSE])))
     }
     n<-rsobj$n
     n_eff<-rsobj$n_eff
     k<-rsobj$k                      # number of columns of matrices
     if(is.null(idx))
         stop("No item(s) for subscale specified (use idx!)")
     res<-switch(stat,
          "var"=rstats(rsobj,T2.Var.stat),
          "mad1"=rstats(rsobj,T2.MAD1.stat),
          "mad2"=rstats(rsobj,T2.MAD2.stat),
          "range"=rstats(rsobj,T2.Range.stat),
          stop("stat must be one of \"var\", \"mad1\", \"mad2\", \"range\"")
     )
     res<-unlist(res)
     prop<-sum(res[2:(n_eff+1)]>=res[1])/n_eff
     result<-list(n_eff=n_eff, prop=prop, idx=idx, stat=stat, T2vec=res) # T2obj
     class(result)<-"T2obj"
     result
}

T2m<-function(rsobj,idx=NULL,stat="var", ...){

     T2m.Var.stat<-function(x){       # calculates statistic T2m
        var(rowSums(x[,idx, drop=FALSE]))
     }
     T2m.MAD1.stat<-function(x){       # calculates statistic T2m
        y<-rowSums(x[,idx, drop=FALSE])           # mean absolute deviation
        mean(abs(y-mean(y)))
     }
     T2m.MAD2.stat<-function(x){       # calculates statistic T2m
        mad(rowSums(x[,idx, drop=FALSE]),constant=1) # unscaled median absolute deviation
     }
     T2m.Range.stat<-function(x){     # calculates statistic T2m
        diff(range(rowSums(x[,idx, drop=FALSE])))
     }
     n<-rsobj$n
     n_eff<-rsobj$n_eff
     k<-rsobj$k                      # number of columns of matrices
     if(is.null(idx))
         stop("No item(s) for subscale specified (use idx!)")
     res<-switch(stat,
          "var"=rstats(rsobj,T2m.Var.stat),
          "mad1"=rstats(rsobj,T2m.MAD1.stat),
          "mad2"=rstats(rsobj,T2m.MAD2.stat),
          "range"=rstats(rsobj,T2m.Range.stat),
          stop("stat must be one of \"var\", \"mad1\", \"mad2\", \"range\"")
     )
     res<-unlist(res)
     prop<-sum(res[2:(n_eff+1)]<=res[1])/n_eff
     result<-list(n_eff=n_eff, prop=prop, idx=idx, stat=stat, T2mvec=res) # T2mobj
     class(result)<-"T2mobj"
     result
}


T4<-function(rsobj,idx=NULL,group=NULL,alternative="high", ...){

     T4.stat<-function(x){      # calculates statistic T4
        sign*sum(rowSums(x[gr,idx,drop=FALSE]))
     }
     n_eff<-rsobj$n_eff                         # number of simulated matrices
     n_tot<-rsobj$n_tot                         # number of all matrices
     k<-rsobj$k                                 # number of items
     if(is.null(idx))
         stop("No item(s) for subscale specified (use idx!)")
     if(length(idx)==k)  # rh 2011-03-03
         stop("Subscale containing all items gives meaningless results for T4.")
     if(is.null(group))
         stop("No group specified (use group!)")
     if(!is.logical(group))     # added rh 2011-03-03
         stop("group must be of type \"logical\" (e.g., group = (age==1) )")
     if(alternative=="high")
        sign <- 1
     else if(alternative=="low")
        sign <- -1
     else
        stop("alternative incorrectly specified! (use either \"high\" or \"low\")")

     gr<-as.logical(group)                      # group definition (logical)
     res<-rstats(rsobj,T4.stat)
     res<-unlist(res)
     prop<-sum(res[2:(n_tot)]>=res[1])/n_eff
     gr.nam <- deparse(substitute(group))
     gr.n <- sum(group)
     result<-list(n_eff=n_eff, prop=prop, idx=idx, gr.nam=gr.nam, gr.n=gr.n, T4vec=res, alternative=alternative)   # T4obj
     class(result)<-"T4obj"
     result
}
# removed in version 0.14-5
#T7<-function(rsobj,idx=NULL, ...){
#     T7.stat<-function(x){      # calculates statistic T7
#        calcT7<-function(i,j){  # calculates sum for all items in subscale
#          if(sitscor[i]>sitscor[j]){
#              sum(submat[,j]>submat[,i])   #
#              # t<-table(submat[,i],submat[,j])    # odds ratio gives the same result
#              # OR<-t[1]*t[4]/(t[2]*t[3])
#              # 1/OR
#          } else
#              NA
#        }
#        submat<-x[,idx]
#        submat<-submat[,order(itscor,decreasing=TRUE)]
#        RET<-unlist(lapply(1:(m-1), function(i) lapply((i+1):m, function(j) calcT7(i,j))))
#        RET
#     }
#
#     n_eff<-rsobj$n_eff                         # number of simulated matrices
#     n_tot<-rsobj$n_tot                         # number of all matrices
#     k<-rsobj$k                                 # number of items
#     if(is.null(idx))
#         stop("No items for subscale specified (use idx!)")
#     else if (length(idx)<2)
#         stop("At least 2 items have to be specified with idx!")
#     submat<-rsextrmat(rsobj,1)[,idx]
#     itscor<-colSums(submat)
#     names(itscor)<-colnames(submat)<-idx
#
#     submat<-submat[,order(itscor,decreasing=TRUE)]
#     sitscor<-sort(itscor,decreasing=TRUE)      # sorted itemscore
#     m<-length(itscor)
#
#     resList<-rstats(rsobj,T7.stat)
#     res<-sapply(resList,sum,na.rm=TRUE)
#     prop<-sum(res[2:(n_eff+1)]>=res[1])/n_eff
#     result<-list(n_eff=n_eff, prop=prop, itscor=itscor, T7vec=res)   # T7obj
#     class(result)<-"T7obj"
#     result
#}
#T7a<-function(rsobj,idx=NULL, ...){
#     T7a.stat<-function(x){      # calculates statistic T7a
#        calcT7a<-function(i,j){  # calculates sum for single Itempair
#          if(sitscor[i]>sitscor[j]){
#              sum(submat[,j]>submat[,i])   #
#              # t<-table(submat[,i],submat[,j])    # odds ratio gives the same result
#              # OR<-t[1]*t[4]/(t[2]*t[3])
#              # 1/OR
#          } else
#              NA
#        }
#        submat<-x[,idx]
#        submat<-submat[,order(itscor,decreasing=TRUE)]
#        RET<-unlist(lapply(1:(m-1), function(i) lapply((i+1):m, function(j) calcT7a(i,j))))
#        RET
#     }
#
#     n_eff<-rsobj$n_eff                         # number of simulated matrices
#     n_tot<-rsobj$n_tot                         # number of all matrices
#     k<-rsobj$k                                 # number of items
#     if(is.null(idx))
#         stop("No items for subscale specified (use idx!)")
#     else if (length(idx)<2)
#         stop("At least 2 items have to be specified with idx!")
#     submat<-rsextrmat(rsobj,1)[,idx]
#     itscor<-colSums(submat)
#     names(itscor)<-colnames(submat)<-idx
#     submat<-submat[,order(itscor,decreasing=TRUE)]
#     sitscor<-sort(itscor,decreasing=TRUE)      # sorted itemscore
#     m<-length(itscor)
#
#     res<-rstats(rsobj,T7a.stat)
#     res<-do.call(cbind, lapply(res,as.vector)) # converts result list to matrix
#     T7avec<-apply(res, 1, function(x) sum(x[2:(n_tot)]>=x[1])/n_eff)
#     T7anam<-NULL
#     for (i in 1:(m-1)) for(j in (i+1):m )
#          T7anam<-c(T7anam, paste("(",names(sitscor[i]),">",names(sitscor[j]),")",sep="",collapse=""))
#     names(T7avec)<-T7anam
#     result<-list(n_eff=n_eff, prop=T7avec,itscor=itscor)    # T7aobj
#     class(result)<-"T7aobj"
#     result
#}

T10<-function(rsobj, splitcr="median", ...){
      calc.groups<-function(x,splitcr){
        if (length(splitcr) > 1)  {        # numeric vectors converted to factors
            if (length(splitcr) != nrow(x)) {
                stop("Mismatch between length of split vector and number of persons!")
            }
            splitcr <- as.factor(splitcr)
            if (length(levels(splitcr))>2) {
                stop("Split vector defines more than 2 groups (only two allowed)!")
            }
            spl.lev <- levels(splitcr)
            #spl.gr <- paste(spl.nam, spl.lev, sep = " ")  # not necessary for the time being
            hi <- splitcr==spl.lev[1] # first level is high group
        } else if (!is.numeric(splitcr)) {
            spl.nam <- splitcr
            if (splitcr == "median") {
                spl.gr <- c("Raw Scores <= Median", "Raw Scores > Median")
                rv <- rowSums(x)
                rvsplit <- median(rv)
                hi <- rv > rvsplit
            }
            if (splitcr == "mean") {
                spl.gr <- c("Raw Scores < Mean", "Raw Scores >= Mean")
                rv <- rowSums(x)
                rvsplit <- mean(rv)
                hi <- rv > rvsplit
            }
        }
        list(hi=hi,spl.nam=spl.nam) # spl.nam is returned due to lex scoping even if not defined here
      }
      T10.stat<-function(x){      # calculates statistic T10 for one matrix
        nij.hi<-unlist(lapply(1:k,function(i) lapply(1:k, function(j) sum(x[hi,i]>x[hi,j]))))
        nij.low<-unlist(lapply(1:k,function(i) lapply(1:k, function(j) sum(x[!hi,i]>x[!hi,j]))))
        nji.hi<- unlist(lapply(1:k,function(i) lapply(1:k, function(j) sum(x[hi,i]<x[hi,j]))))
        nji.low<- unlist(lapply(1:k,function(i) lapply(1:k, function(j) sum(x[!hi,i]<x[!hi,j]))))
        RET<-sum(abs(nij.hi*nji.low-nij.low*nji.hi))
        RET
      }
      spl.nam <- deparse(substitute(splitcr))
      n_eff<-rsobj$n_eff                         # number of simulated matrices
      n_tot<-rsobj$n_tot                         # number of all matrices
      k<-rsobj$k                                 # number of columns of matrices
      obj<-rsextrobj(rsobj,1,1)                  # extract first matrix
      x<-matrix(obj$inpmat,obj$n,obj$k)
      ans <- calc.groups(x,splitcr)      # calculate grouping vector (logical)
      hi<-ans$hi
      hi.n<-sum(hi)
      low.n<-sum(!hi)

      res<-rstats(rsobj,T10.stat)                # for each matrix calculate T10

      res<-unlist(res)
      prop<-sum(res[2:(n_eff+1)]>=res[1])/n_eff
      result<-list(n_eff=n_eff, prop=prop,spl.nam=ans$spl.nam,hi.n=hi.n,low.n=low.n,T10vec=res)  # T10obj
      class(result)<-"T10obj"
      result
}


T11<-function(rsobj, ...){
      T11.stat<-function(x){
         as.vector(cor(x))
      }
      calc.T11<-function(x){      # calculates statistic T11 for one matrix
         sum(abs(x-rho))
      }
      n_eff<-rsobj$n_eff                         # number of simulated matrices
      n_tot<-rsobj$n_tot                         # number of all matrices
      k<-rsobj$k                                 # number of columns of matrices
      res<-rstats(rsobj,T11.stat)                # for each matrix calculate all r_ij's

      cormats <- matrix(unlist(res),nrow=k*k)    # k*k x n_tot matrix, each colum contains one corr matrix
      rho<-apply(cormats[,2:n_tot],1,mean)       # vector of estimated "real" rho_ij's
      T11obs<-calc.T11(cormats[,1])              # vector of observed r_ij's
      prop<-sum(apply(cormats[, 2:n_tot],2,calc.T11)>=T11obs)/n_eff
      result<-list(n_eff=n_eff, prop=prop, T11r=cormats[,1], T11rho=rho)   # T11obj
      class(result)<-"T11obj"
      result
}

## The following two functions were included in version 0.16-3
Q3h<-function(rsobj, ...){
    Q3h.stat <- function(x){
      as.vector(x)
    }
  
    calcQ3h.stat <- function(x, exp=exp) { ## Calculates Q3h based on observed matrix and expected values
      # Calculates Q3
	  i <- ncol(x)
      mat <- x - exp
      res <- matrix(nrow=i,ncol=i)
      for(a in 1:(i-1)) {
        for(b in (a+1):i) {
          res[b,a] <- res[a,b] <- -cor(mat[,a],mat[,b])
        }
      }
      return(res)
    }

    n_eff<-rsobj$n_eff                         # number of simulated matrices
    n_tot<-rsobj$n_tot                         # number of simulated matrices
    k <- rsobj$k                               # number of columns of matrices
	n <- rsobj$n                               # number of rows of matrices

	# 1st step of calculating Q3h: Calculate the expected values
    res <- rstats(rsobj,Q3h.stat) # res contains vector with 1st column, 2nd column etc of each matrix as entries
    datmat <- matrix(unlist(res),nrow=k*n) # Contains the entries (columns-wise) of each matrix in each column
    exp <- matrix(apply(datmat, 1, mean), nrow=n, ncol=k)
  
    # 2nd step: Calculate Q3h based on the simulated matrices and the expected values 
    res <- rstats(rsobj, calcQ3h.stat, exp=exp)

	# 3rd step: Calculate p-values (analogous to T1, T1m and T1l)
	res<-do.call(cbind, lapply(res,as.vector)) # converts result list to matrix
    Q3hvec<-apply(res, 1, function(x) sum(x[2:(n_tot)]<x[1])/n_eff)
    Q3hmat<- matrix(Q3hvec, ncol=k)
    Q3hmat[upper.tri(Q3hmat)] <- NA # For consistency with other nonparametric tests
    Q3hvec <- as.vector(Q3hmat)
	  Q3hvec <- Q3hvec[!is.na(Q3hvec)]# For consistency with other nonparametric tests
    result<-list(n_eff=n_eff, prop=Q3hvec, Q3hmat=Q3hmat) # Q3hobj
    class(result)<-"Q3hobj"
    return(result)
}

Q3l<-function(rsobj, ...){
    Q3l.stat <- function(x){
      as.vector(x)
    }
  
    calcQ3l.stat <- function(x, exp=exp) { ## Calculates Q3l based on observed matrix and expected values
      # Calculates Q3
	  i <- ncol(x)
      mat <- x - exp
      res <- matrix(nrow=i,ncol=i)
      for(a in 1:(i-1)) {
        for(b in (a+1):i) {
          res[b,a] <- res[a,b] <- cor(mat[,a],mat[,b])
        }
      }
      return(res)
    }

    n_eff<-rsobj$n_eff                         # number of simulated matrices
    n_tot<-rsobj$n_tot                         # number of simulated matrices
    k <- rsobj$k                               # number of columns of matrices
	n <- rsobj$n                               # number of rows of matrices

	# 1st step of calculating Q3l: Calculate the expected values
    res <- rstats(rsobj,Q3l.stat) # res contains vector with 1st column, 2nd column etc of each matrix as entries
    datmat <- matrix(unlist(res),nrow=k*n) # Contains the entries (columns-wise) of each matrix in each column
    exp <- matrix(apply(datmat, 1, mean), nrow=n, ncol=k)
  
    # 2nd step: Calculate Q3l based on the simulated matrices and the expected values 
    res <- rstats(rsobj, calcQ3l.stat, exp=exp)

	# 3rd step: Calculate p-values (analogous to T1, T1m and T1l)
	res<-do.call(cbind, lapply(res,as.vector)) # converts result list to matrix
	Q3lvec<-apply(res, 1, function(x) sum(x[2:(n_tot)]<x[1])/n_eff)
	Q3lmat<- matrix(Q3lvec, ncol=k)
	Q3lmat[upper.tri(Q3lmat)] <- NA # For consistency with other nonparametric tests
	Q3lvec <- as.vector(Q3lmat)
	Q3lvec <- Q3lvec[!is.na(Q3lvec)]# For consistency with other nonparametric tests
	result<-list(n_eff=n_eff, prop=Q3lvec, Q3lmat=Q3lmat) # Q3lobj
  class(result)<-"Q3lobj"
  return(result)
}

## End of code included in version 0.16-3

print.MLobj<-function(x,...){
  print(x$MLres)
  cat("'exact' p-value =", x$prop, " (based on", x$n_eff, "sampled matrices)\n\n")
}

print.Tmdobj<-function(x,...){
  txt1<-"\nNonparametric RM model test: Tmd (Multidimensionality)"
  writeLines(strwrap(txt1, exdent=4))
  cat("    (correlation of subscale person scores)\n")
  cat("Number of sampled matrices:", x$n_eff,"\n")
  cat("Subscale 1 - Items:", x$idx1,"\n")
  cat("Subscale 2 - Items:", x$idx2,"\n")
  cat("Observed correlation:", x$Tmdvec[1],"\n")
  cat("one-sided p-value:",x$prop,"\n\n")
}

print.Tpbisobj<-function(x,...){
  txt1<-"\nNonparametric RM model test: Tpbis (discrimination)"
  writeLines(strwrap(txt1, exdent=4))
  cat("    (pointbiserial correlation of test item vs. subscale)\n")
  cat("Number of sampled matrices:", x$n_eff,"\n")
  cat("Test Item:", x$idxt,"\n")
  cat("Subscale  - Items:", x$idxs,"\n")
  cat("one-sided p-value (rpbis too low):",x$prop,"\n\n")
}

print.T1obj<-function(x,alpha=0.05,...){
  txt1<-"\nNonparametric RM model test: T1 (local dependence - increased inter-item correlations)\n"
  writeLines(strwrap(txt1, exdent=4))
  cat("    (counting cases with equal responses on both items)\n")
  cat("Number of sampled matrices:", x$n_eff,"\n")
  cat("Number of Item-Pairs tested:", length(x$prop),"\n")
  cat("Item-Pairs with one-sided p <", alpha,"\n")
  T1mat<-x$T1mat
  idx<-which(T1mat<alpha,arr.ind=TRUE)
  val<-T1mat[which(T1mat<alpha)]
  names(val)<-apply(idx,1,function(x) paste("(",x[2],",",x[1],")",sep="",collapse=""))
  if (length(val)>0)
     print(round(val,digits=3))
  else
     cat("none\n\n")
}

print.T1mobj<-function(x,alpha=0.05,...){
  txT1m<-"\nNonparametric RM model test: T1m (multidimensionality - reduced inter-item correlations)\n"
  writeLines(strwrap(txT1m, exdent=4))
  cat("    (counting cases with equal responses on both items)\n")
  cat("Number of sampled matrices:", x$n_eff,"\n")
  cat("Number of Item-Pairs tested:", length(x$prop),"\n")
  cat("Item-Pairs with one-sided p <", alpha,"\n")
  T1mmat<-x$T1mmat
  idx<-which(T1mmat<alpha,arr.ind=TRUE)
  val<-T1mmat[which(T1mmat<alpha)]
  names(val)<-apply(idx,1,function(x) paste("(",x[2],",",x[1],")",sep="",collapse=""))
  if (length(val)>0)
     print(round(val,digits=3))
  else
     cat("none\n\n")
}

print.T1lobj<-function(x,alpha=0.05,...){
  txt1<-"\nNonparametric RM model test: T1 (learning - based on item pairs)\n"
  writeLines(strwrap(txt1, exdent=4))
  cat("    (counting cases with reponsepattern (1,1) for item pair)\n")
  cat("Number of sampled matrices:", x$n_eff,"\n")
  cat("Number of sampled matrices:", x$n_eff,"\n")
  cat("Number of Item-Pairs tested:", length(x$prop),"\n")
  cat("Item-Pairs with one-sided p <", alpha,"\n")
  T1lmat<-x$T1lmat
  idx<-which(T1lmat<alpha,arr.ind=TRUE)
  val<-T1lmat[which(T1lmat<alpha)]
  names(val)<-apply(idx,1,function(x) paste("(",x[2],",",x[1],")",sep="",collapse=""))
  if (length(val)>0)
     print(round(val,digits=3))
  else
     cat("none\n\n")
}

print.T2obj<-function(x,...){
  prop<-x$prop
  idx<-x$idx
  stat<-x$stat
  statnam<-switch(stat,
     "var"="variance",
     "mad1"="mean absolute deviation",
     "mad2"="median absolute deviation",
     "range"="range"
  )
  txt<-"\nNonparametric RM model test: T2 (local dependence - model deviating subscales)\n"
  writeLines(strwrap(txt, exdent=4))
  cat("    (increased dispersion of subscale person rawscores)\n")
  cat("Number of sampled matrices:", x$n_eff,"\n")
  cat("Items in subscale:", idx,"\n")
  cat("Statistic:", statnam,"\n")
  cat("one-sided p-value:",prop,"\n\n")
#  cat("    (proportion of sampled",statnam," GE observed)\n\n")
}

print.T2mobj<-function(x,...){
  prop<-x$prop
  idx<-x$idx
  stat<-x$stat
  statnam<-switch(stat,
     "var"="variance",
     "mad1"="mean absolute deviation",
     "mad2"="median absolute deviation",
     "range"="range"
  )
  txt<-"\nNonparametric RM model test: T2m (multidimensionality - model deviating subscales)\n"
  writeLines(strwrap(txt, exdent=4))
  cat("    (decreased dispersion of subscale person rawscores)\n")
  cat("Number of sampled matrices:", x$n_eff,"\n")
  cat("Items in subscale:", idx,"\n")
  cat("Statistic:", statnam,"\n")
  cat("one-sided p-value:",prop,"\n\n")
#  cat("    (proportion of sampled",statnam," GE observed)\n\n")
}
print.T4obj<-function(x,...){
  prop<-x$prop
  idx<-x$idx
  gr.nam<-x$gr.nam
  gr.n<-x$gr.n
  alternative<-x$alternative
  cat("\nNonparametric RM model test: T4 (Group anomalies - DIF)\n")
  txt<-paste("    (counting", alternative, "raw scores on item(s) for specified group)\n", collapse="")
  writeLines(strwrap(txt, exdent=4))
  cat("Number of sampled matrices:", x$n_eff,"\n")
  cat("Items in Subscale:", idx,"\n")
  cat("Group:",gr.nam,"  n =",gr.n,"\n")
  cat("one-sided p-value:",prop,"\n\n")
#  cat("    (proportion of sampled raw scores GE observed)\n\n")
}

# removed in version 0.14-5
#print.T7obj<-function(x,...){
#  prop<-x$prop
#  cat("\nNonparametric RM model test: T7 (different discrimination - 2PL)\n")
#  txt<-"    (counting cases with response 1 on more difficult and 0 on easier item)\n"
#  writeLines(strwrap(txt, exdent=4))
#  cat("Number of sampled matrices:", x$n_eff,"\n")
#  cat("Item Scores:\n")
#  print(x$itscor)
#  cat("one-sided p-value:",prop,"\n\n")
#}
#print.T7aobj<-function(x,...){
#  prop<-x$prop
#  cat("\nNonparametric RM model test: T7a (different discrimination - 2PL)\n")
#  txt<-"    (counting cases with response 1 on more difficult and 0 on easier item)\n"
#  writeLines(strwrap(txt, exdent=4))
#  cat("Number of sampled matrices:", x$n_eff,"\n")
#  cat("Item Scores:\n")
#  print(x$itscor)
#  cat("\nItem-Pairs: (i>j ... i easier than j)\n\n")
#  print(round(prop,digits=3))
#}
print.T10obj<-function(x,...){
  spl.nam<-x$spl.nam
  prop<-x$prop
  hi.n<-x$hi.n
  low.n<-x$low.n
  txt<-"\nNonparametric RM model test: T10 (global test - subgroup-invariance)\n"
  writeLines(strwrap(txt, exdent=4))
  cat("Number of sampled matrices:", x$n_eff,"\n")
  cat("Split:",spl.nam,"\n")
  cat("Group 1: n = ",hi.n,"  Group 2: n =",low.n,"\n")
  cat("one-sided p-value:",prop,"\n\n")
#  cat("    (proportion of sampled statistics GE observed)\n\n")
}
print.T11obj<-function(x,...){
  prop<-x$prop
  txt<-"\nNonparametric RM model test: T11 (global test - local dependence)\n"
  writeLines(strwrap(txt, exdent=4))
  txt<-"    (sum of deviations between observed and expected inter-item correlations)\n"
  writeLines(strwrap(txt, exdent=4))
  cat("Number of sampled matrices:", x$n_eff,"\n")
  cat("one-sided p-value:",prop,"\n\n")
#  cat("    (proportion of sampled sums GE observed)\n\n")
}

## The following code was in included in version 0.16-3
print.Q3hobj<-function(x,alpha=0.05,...){
  txt1<-"\nNonparametric RM model test: Q3h (local dependence - increased correlation of inter-item residuals)\n"
  writeLines(strwrap(txt1, exdent=4))
  cat("Number of sampled matrices:", x$n_eff,"\n")
  cat("Number of Item-Pairs tested:", length(x$prop),"\n")
  cat("Item-Pairs with one-sided p <", alpha,"\n")
  Q3hmat<-x$Q3hmat
  idx<-which(Q3hmat<alpha,arr.ind=TRUE)
  val<-Q3hmat[which(Q3hmat<alpha)]
  names(val)<-apply(idx,1,function(x) paste("(",x[2],",",x[1],")",sep="",collapse=""))
  if (length(val)>0)
     print(round(val,digits=3))
  else
     cat("none\n\n")
}

print.Q3lobj<-function(x,alpha=0.05,...){
  txt1<-"\nNonparametric RM model test: Q3l (local dependence - decreased correlation of inter-item residuals)\n"
  writeLines(strwrap(txt1, exdent=4))
  cat("Number of sampled matrices:", x$n_eff,"\n")
  cat("Number of Item-Pairs tested:", length(x$prop),"\n")
  cat("Item-Pairs with one-sided p <", alpha,"\n")
  Q3lmat<-x$Q3lmat
  idx<-which(Q3lmat<alpha,arr.ind=TRUE)
  val<-Q3lmat[which(Q3lmat<alpha)]
  names(val)<-apply(idx,1,function(x) paste("(",x[2],",",x[1],")",sep="",collapse=""))
  if (length(val)>0)
    print(round(val,digits=3))
  else
    cat("none\n\n")
}

## End of new code