1ME.Lenth <- function(b, simulated=TRUE, alpha=NULL){
2    s0 <- 1.5 * median(abs(b))
3    cj <- as.numeric(b[abs(b) < 2.5 * s0])
4    PSE <- 1.5 * median(abs(cj))
5    m <- length(b)
6    alphas <- as.numeric(colnames(crit.ME))
7      if (is.null(alpha)) alpha <- alphas
8      else{
9         if (!is.numeric(alpha)) stop("alpha must be numeric")
10         if (min(alpha)<0 || max(alpha)>1) stop("alpha must be between 0 and 1")
11      }
12      if (simulated){
13        if (length(setdiff(alpha, alphas))>0 || m<7 || m>143){
14        simulated <- FALSE
15        message("simulated critical values not available for all requests, used conservative ones")
16        }
17      }
18    ## cover simulated critical values, where available
19    if (simulated){
20      ME <- crit.ME[as.character(m),as.character(alpha)]*PSE
21      SME <- crit.SME[as.character(m),as.character(alpha)]*PSE
22    }
23    else{
24      gamma <- (1 + (1 - alpha)^(1/m))/2
25      ME <- qt(1-alpha/2, m/3)*PSE
26      SME <- qt(gamma, m/3)*PSE
27      names(ME) <- names(SME) <- alpha
28    }
29    list(s0=s0, PSE=PSE, ME=ME, SME=SME)
30  }
31
32critsimul <- function(m, dfe, nsimul=10000, alpha=round(seq(0.25,0.01,-0.01),2),
33                      type="EM08", weight0=5){
34## m the number of effects
35## dfe the number of error df
36## type "EM08" or "LW98"
37## weight0 relevant for EM08 only
38   ergs.CME <- rep(NA, m*dfe)
39   ergs.CSME <- rep(NA, dfe)
40
41   for (i in 1:nsimul){
42      x <- rnorm(m)                   ## all effects inactive
43      se <- sqrt(rchisq(1,dfe)/dfe)   ## estimated standard error
44      if (type=="EM08") hilf <- CME.EM08(x, se, dfe, weight0=weight0, alpha=alpha,
45                    simulated=FALSE)
46      else hilf <- CME.LW98(x, se, dfe, alpha=alpha,
47                    simulated=FALSE)
48      hilf <- x/hilf$CPSE
49      ergs.CME[((i-1)*m+1):(i*m)] <- hilf
50      ergs.CSME[i] <- max(hilf)
51   }
52   CME.mult <- quantile(abs(ergs.CME), 1-alpha)
53   CSME.mult <- quantile(abs(ergs.CSME), 1-alpha)
54   names(CME.mult) <- alpha
55   names(CSME.mult) <- alpha
56   list(CME.mult = CME.mult, CSME.mult = CSME.mult)
57}
58
59CME.LW98 <- function(b, sterr, dfe, simulated=TRUE, alpha=NULL){
60  ## sterr = sqrt(K*MSE), obtainable from linear model analysis
61  ## dfe = residual df
62  s0 <- 1.5 * median(abs(b))
63  cj <- as.numeric(b[abs(b) < 2.5 * s0])
64  m <- length(b)
65  d <- m/3
66  wPSE <- d / (d + dfe)
67  CPSE <- sqrt((1.5 * median(abs(cj)))^2*wPSE + (1-wPSE)*sterr^2)
68    alphas <- as.numeric(colnames(crit.ME))
69      if (is.null(alpha)) alpha <- alphas
70      else{
71         if (!is.numeric(alpha)) stop("alpha must be numeric")
72         if (min(alpha)<0 || max(alpha)>1) stop("alpha must be between 0 and 1")
73      }
74  if (simulated){
75    hilf <- critsimul(m, dfe, type="LW98", alpha=alpha)
76    CME <- hilf$CME.mult*CPSE
77    CSME <- hilf$CSME.mult*CPSE
78  }
79  else{
80    gamma <- (1 + (1 - alpha)^(1/m))/2
81    CME <- qt(1-alpha/2, d+dfe)*CPSE
82    CSME <- qt(gamma, d+dfe)*CPSE
83    names(CME) <- names(CSME) <- alpha
84  }
85  list(s0=s0, CPSE=CPSE, CME=CME, CSME=CSME)
86}
87
88CME.EM08 <- function(b, sterr, dfe, simulated=TRUE, weight0=5, alpha=NULL){
89  ## sterr = sqrt(K*MSE), obtainable from linear model analysis
90  ## dfe = residual df
91  ## weight0 <- multiplier for dfe in weighting sterr^2 for s0
92  m <- length(b)
93  d <- length(b)/3
94  ws0 <- d / (d + weight0 * dfe)
95  wPSE <- d / (d + dfe)
96  s0 <- 1.5 * median(abs(b))
97  s0 <- sqrt(ws0*s0^2 + (1-ws0)*sterr^2)
98  cj <- as.numeric(b[abs(b) < 2.5 * s0])
99  CPSE <- sqrt((1.5 * median(abs(cj)))^2*wPSE + (1-wPSE)*sterr^2)
100    alphas <- as.numeric(colnames(crit.ME))
101      if (is.null(alpha)) alpha <- alphas
102      else{
103         if (!is.numeric(alpha)) stop("alpha must be numeric")
104         if (min(alpha)<0 || max(alpha)>1) stop("alpha must be between 0 and 1")
105      }
106  if (simulated){
107    hilf <- critsimul(m, dfe, type="EM08", weight0=weight0, alpha=alpha)
108    CME <- hilf$CME.mult*CPSE
109    CSME <- hilf$CSME.mult*CPSE
110  }
111  else{
112    gamma <- (1 + (1 - alpha)^(1/m))/2
113    CME <- qt(1-alpha/2, d+dfe)*CPSE
114    CSME <- qt(gamma, d+dfe)*CPSE
115    names(CME) <- names(CSME) <- alpha
116  }
117  list(Cs0=s0, CPSE=CPSE, CME=CME, CSME=CSME)
118}
119