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