### --------------------------------------------------------------------------------------### ### Winter Term 2017/2018 ### ### NMSA407 - Linear Regression ### ### ### ### EXERCISE #11 Multiple comparison procedures ### ### ### ### --------------------------------------------------------------------------------------### ### Working directory setwd("H:/nmsa407_LinRegr/") # To be changed to something which exists on your PC. rm(list=ls()) ### Load needed packages library("mffSM") library("lmtest") library("multcomp") ### loading the data file (the working directory has to be set appropriately) load("mana.RData") ### D A T A D E S C R I P T I O N ### Data comes from the project MANA (Maturita na necisto) in 2005. This project was a part ### of bigger project which goal was to prepare a new form of graduation exam ("statni maturita"). ### This dataset contains the results in English language at grammar schools. ### ### VARIABLES ### region: Region of the Czech Republic, ### 1 = Praha ### 2 = Stredocesky ### 3 = Jihocesky ### 4 = Plzensky ### 5 = Karlovarsky ### 6 = Ustecky ### 7 = Liberecky ### 8 = Kralovehradecky ### 9 = Pardubicky ### 10 = Vysocina ### 11 = Jihomoravsky ### 12 = Olomoucky ### 13 = Zlinsky ### 14 = Moravskoslezsky ### ### specialization: The specialization of the university where the student is going to apply. ### 1 = education ### 2 = social science ### 3 = languages ### 4 = law ### 5 = math-physics ### 6 = engineering ### 7 = science ### 8 = medicine ### 9 = economics ### 10 = agriculture ### 11 = art ### ### gender: Gender of the student, ### 0 = female, 1 = male ### ### graduation Is the student going to graduate in English? ### 0 = no, 1 = yes ### ### entry.exam: Is an exam on English knowledge part of an entry exam at the ### university (where the student is going to apply)? ### 0 = no, 1 = yes ### ### score: The score in English language in MANA test. ### ### avg.mark: Average mark (grade) from all subjects at the last report card ### (vysvedceni). ### ### mark.eng: Average mark (grade) from the English language at the last ### 5 report cards (vysvedceni). head(mana) dim(mana) summary(mana) # some necessary transformations with(mana, table(region)) with(mana, table(specialization)) with(mana, table(gender)) with(mana, table(graduation)) with(mana, table(entry.exam)) ### Instead of numbers the following coding of the regions will be used ### ### A - Capital City Prague (Praha) ### B - South Moravia - Jihomoravsky kraj (Brno) ### C - South Bohemia - Jihocesky kraj (Ceske Budejovice) ### E - Pardubice - Pardubicky kraj ### H - Hradec Kralove - Kralovehradecky kraj ### J - Highland (Vysocina) Region - Kraj Vysocina (Jihlava) ### K - Karlovy Vary - Karlovarsky kraj ### L - Liberec - Liberecky kraj ### M - Olomouc - Olomoucky kraj ### P - Plzen - Plzensky kraj ### S - Central Bohemia (Prague) - Stredocesky kraj ### T - Moravia-Silesia - Moravskoslezsky kraj (Ostrava) ### U - Usti nad Labem - Ustecky kraj ### Z - Zlin - Zlinsky kraj mana <- transform(mana, fregion = factor(region, levels = 1:14, labels = c("A", "S", "C", "P", "K", "U", "L", "H", "E", "J", "B", "M", "Z", "T")), fspec = factor(specialization, levels = 1:11, labels = c("Educ", "Social", "Lang", "Law", "Mat-Phys", "Engin", "Science", "Medic", "Econom", "Agric", "Art")), fgender = factor(gender, levels = 0:1, labels = c("Female", "Male")), fgrad = factor(graduation, levels = 0:1, labels = c("No", "Yes")), fentry = factor(entry.exam, levels = 0:1, labels = c("No", "Yes"))) summary(mana) with(mana, table(fregion)) with(mana, table(fspec)) ### QUESTIONS: ### Are there differences among regions in grading (in quality of education?) as expressed ### by average mark on the last certificate [avg.mark] when the prospective specialization is taken into account? ### Some plots lSPEC <- levels(mana[, "fspec"]) colSPEC <- diverge_hcl(length(lSPEC)) names(colSPEC) <- lSPEC lREG <- levels(mana[, "fregion"]) colREG <- diverge_hcl(length(lREG)) names(colREG) <- lREG with(mana, interaction.plot(fregion, fspec, avg.mark, col = colSPEC, type = "b", lwd = 2, lty = 1, xlab = "Region", trace.label = "Specialization", ylab = "Mean of avg.mark")) with(mana, interaction.plot(fspec, fregion, avg.mark, col = colREG, type = "b", lwd = 2, lty = 1, xlab = "Specialization", trace.label = "Region", ylab = "Mean of avg.mark")) ### Not really insightful with so many groups... ### But still, do you think that an additive model will be a "wrong but useful" model? ### Counts for each combination... addmargins(with(mana, table(fregion, fspec))) ### Note, there are some zeros in the Agric column. ### Are there any consequences for the two-way ANOVA linear model with interactions? ### Simple interaction model mInter <- lm(avg.mark ~ fregion*fspec, data = mana) summary(mInter) ### QUESTIONS: ### Do you have any idea why those NA's occur in the output? ### Is this a usefull model in general? ### Is it worth of going for some simplification of the model? ### anova(mInter) Anova(mInter) ### F-tests on submodels are OK... However... ### ### Interaction just not significant, but with our (quite huge) sample size, ### can we assume that the interaction terms are practically unimportant? ### Check again the second spaghetti plot... ### Additive model mAddit <- lm(avg.mark ~ fregion + fspec, data = mana) summary(mAddit) ### QUESTIONS: ### Why there are no NA's observations now? ################### ### HW task - start mAdditSum <- lm(avg.mark ~ fregion + fspec, data = mana, contrasts = list(fregion = "contr.sum", fspec = "contr.sum")) summary(mAdditSum) ### Interpretation of intercept here? ### HW task - end ################# ### Basic diagnostics ### Is model "wrong but useful" for inference? plotLM(mAddit) ### QUESTIONS: ### Do we seem to have serious problems with the model? ### What about the Breusch-Pagan test for homoscedasticity? ncvTest(mAddit, var.formula = ~fitted(mAddit)) bptest(mAddit, varformula = ~fitted(mAddit)) ### Or the Koenker's studentized version instead bptest(mAddit) ### Homoscedasticity rejected only for the second version of Breusch-Pagan test. Do not forget how large is the sample size. ### Shapiro-Wilk test for normality of errors based on raw residuals shapiro.test(residuals(mAddit)) ## Implemented algorithm does not allow for more than 5000 observations. ## Are we actually really interested in testing normality ## with such huge sample size? ### Another commonly used test of normality is D'Agostino test dagostino.test <- function(x) { DNAME <- deparse(substitute(x)) x <- x[complete.cases(x)] n <- length(x) if (n<6) stop("sample size must be at least 6") meanX <- mean(x) s<- sqrt(mean((x-meanX)**2)) a3 <- mean((x-meanX)**3)/s**3 a4 <- mean((x-meanX)**4)/s**4 SD3 <- sqrt(6*(n-2)/((n+1)*(n+3))) SD4 <- sqrt(24*(n-2)*(n-3)*n/((n+1)**2*(n+3)*(n+5))) U3 <- a3/SD3 U4 <- (a4-3+6/(n+1))/SD4 b <-(3*(n**2+27*n-70)*(n+1)*(n+3))/((n-2)*(n+5)*(n+7)*(n+9)) W2 <- sqrt(2*(b-1))-1 delta <- 1/sqrt(log(sqrt(W2))) a <- sqrt(2/(W2-1)) Z3 <- delta*log((U3/a)+sqrt((U3/a)**2+1)) B <- (6*(n*n-5*n+2)/((n+7)*(n+9)))*sqrt((6*(n+3)*(n+5))/(n*(n-2)*(n-3))) A <- 6+(8/B)*((2/B)+sqrt(1+4/(B**2))) jm <- sqrt(2/(9*A)) pos <- ((1-2/A)/(1+U4*sqrt(2/(A-4))))**(1/3) Z4 <- (1-2/(9*A)-pos)/jm omni <- Z3**2+Z4**2 pZ3 <- 2*(1-pnorm(abs(Z3),0,1)) pZ4 <- 2*(1-pnorm(abs(Z4),0,1)) pomni <- 1-pchisq(omni,2) skewness <- c(Z3,pZ3) kurtosis <- c(Z4,pZ4) omnibus <- c(omni,pomni) DA <- cbind(skewness,kurtosis,omnibus) row.names(DA)<-c("statistics","p-value") return(DA) } dagostino.test(residuals(mAddit)) ### Alternatively we can use the R library 'fBasics' and the D'Agostino test ### implemented in this library. ### ### library(fBasics); ### dagoTest(residuals(mAddit)) # for more information see https://ncss-wpengine.netdna-ssl.com/wp-content/themes/ncss/pdf/Procedures/NCSS/Normality_Tests.pdf ### Tests for partial effects of region and specialization Anova(mAddit, type = "II") ### TukeyHSD # see Section 11.3.3 of the lecture notes ### Weighted means of the group means are being compared! ### Is this something we are interested in? ### Reported estimated differences are NOT what we call partial effects! aAddit <- aov(avg.mark ~ fregion + fspec, data = mana) TukeyHSD(aAddit) all.equal(coef(aAddit),coef(mAddit)) ### Which regions are different (given specialization)? (treg <- TukeyHSD(aAddit)[["fregion"]]) (treg <- treg[order(treg[, "p adj"]),]) treg[treg[, "p adj"] < 0.05,] ### Compare the results with ### (treg2 <- TukeyHSD(aov(avg.mark ~ fregion, data = mana))[["fregion"]]) ### (treg2 <- treg2[order(treg2[, "p adj"]),]) ### treg2[treg2[, "p adj"] < 0.05,] sort(reg.means <- tapply(mana$avg.mark, mana$fregion, mean)); (diff.reg.means <- outer(reg.means, reg.means, "-")); diff.reg.means["E","C"] treg[treg[, "p adj"] < 0.05,][1,1] treg2[treg2[, "p adj"] < 0.05,][1,1] # in one-way model G = nrow(diff.reg.means) n = nrow(mana) nG = table(mana$fregion) s2 = sum(aov(avg.mark ~ fregion, data = mana)$res^2)/(n-G) diff.reg.means["E", "C"] + c(-1,1)*qtukey(0.95,G,n-G)*sqrt(1/2*(1/nG["E"]+1/nG["C"])*s2) treg2[treg2[, "p adj"] < 0.05,][1,] # in two-way model H = length(levels(mana$fspec)) s2 = sum(aAddit$res^2)/(n-G-H+1) diff.reg.means["E", "C"] + c(-1,1)*qtukey(0.95,G,n-G-H+1)*sqrt(1/2*(1/nG["E"]+1/nG["C"])*s2) treg[treg[,"p adj"] < 0.05,][1,2:3] ### Hothorn-Bretz-Westfall MCP ### Now, we provide estimates of partial effects ### (differences in means of avg.mark between regions given specialization). ### What we are trying to do? ### Used pseudocontrast matrices ### and the "L" matrix to calculate estimated group means ### (under the assumption of the additive model) coef(mAddit) ### Contrast matrix for the regions CX <- contr.treatment(length(lREG)) rownames(CX) <- lREG colnames(CX) <- names(coef(mAddit))[2:length(lREG)] print(CX) ### Contrast matrix for the specializations CZ <- contr.treatment(length(lSPEC)) rownames(CZ) <- lSPEC colnames(CZ) <- names(coef(mAddit))[(length(lREG) + 1):(length(coef(mAddit)))] print(CZ) ### Estimating group means under the assumption of the additive model LmuAddit <- cbind(1, kronecker(rep(1, length(lSPEC)), CX), kronecker(CZ, rep(1, length(lREG)))) colnames(LmuAddit) <- names(coef(mAddit)) rownames(LmuAddit) <- outer(lREG, lSPEC, paste, sep = ":") print(LmuAddit[1:3,]) print(LmuAddit[15:17,]) LSest(mAddit, L = LmuAddit) ### Estimated group means under the assumption of the additive model muAddit <- matrix(LSest(mAddit, L = LmuAddit)[["Estimate"]], nrow = length(lREG), ncol = length(lSPEC)) rownames(muAddit) <- lREG colnames(muAddit) <- lSPEC print(muAddit) ### Observed group means in a plot with(mana, interaction.plot(fspec, fregion, avg.mark, col = colREG, type = "b", lwd = 2, lty = 1, xlab = "Specialization", trace.label = "Region", ylab = "Mean of avg.mark")) ### Model based group means in a plot dfrAddit <- data.frame(avg.mark = as.numeric(muAddit), fregion = rep(lREG, length(lSPEC)), fspec = rep(lSPEC, each = length(lREG))) print(dfrAddit) dfrAddit$fspec = factor(dfrAddit$fspec, levels=levels(mana$fspec)) with(dfrAddit, interaction.plot(fspec, fregion, avg.mark, col = colREG, type = "b", lwd = 2, lty = 1, xlab = "Specialization", pch=levels(fregion), cex=0.5, trace.label = "Region", ylab = "Model based mean avg.mark")) ### QUESTIONS: ### What is the link between observed and model based estimated group means? ### What are we trying to do when estimating pairwise differences with respect to region ### (partial effects of region given specialization)? ### Do you see ordering of regions with respect to estimated means of the means of the [avg.mark] ### (if averaging is done over specializations)? ### Core part of the "L" matrix to calculate all pairwise differences with respect to region. Lcore <- matrix(nrow = 0, ncol = ncol(CX)) colnames(Lcore) <- colnames(CX) print(Lcore) ## matrix with 0 rows rname <- character(0) for (g1 in 1:(length(lREG) - 1)){ for (g2 in (g1 + 1):length(lREG)){ Lcore <- rbind(Lcore, CX[g2,] - CX[g1,]) rname <- c(rname, paste(lREG[g2], "-", lREG[g1])) } } rownames(Lcore) <- rname dim(Lcore) length(lREG) choose(length(lREG),2) print(Lcore) ### QUESTSIONS: ### If other than the reference group pseudocontrasts matrix ### was used, would the above calculation be different? How? ### Complete "L" matrix to calculate all pairwise differences L <- cbind(0, Lcore, matrix(0, nrow = nrow(Lcore), ncol = length(lSPEC) - 1)) colnames(L) <- names(coef(mAddit)) rownames(L) <- rownames(Lcore) print(L) ### Unadjusted estimation of partial effects temp <- LSest(mAddit, L = L) class(temp) <- c("data.frame") temp <- round(temp,4) temp <- temp[order(temp[["P value"]]),]; temp[temp[["P value"]] <= 0.05,] temp[nrow(temp),] # T-S not significant ### Bretz-Hothorn-Westfall procedure - all functions from package multcomp mcp <- glht(mAddit, linfct = L, rhs = 0) print(mcp) temp["T - Z",] ### The above can be done more easily with the help of the built-in option mcp <- glht(mAddit, linfct = mcp(fregion = "Tukey"), rhs = 0) print(mcp) ### This is the same, "Tukey" in the function call ### is just a keyword to ask for all pairwise differences ### with respect to factor fregion. With respect to underlying ### theory, Tukey's procedure is not involved at all in the ### calculation! ### Using a sandwich estimator to adjust also for possible heteroscedasticity library(sandwich) mcp_sandwich <- glht(mAddit, linfct = mcp(fregion = "Tukey"), rhs = 0, vcov=sandwich) print(mcp_sandwich) ### P-values adjusted for MCP will now be calculated. ### Multivariate integrals now must be computed by a Monte Carlo method. ### This may take some time... (load("nmsa407-mult-comparison.RData")) ## path to this file might need some adjustment ### set.seed(20010911) ### system.time(smcp <- summary(mcp)) # (cca 17 sec. on i5-4570 CPU @ 3.20GH, 2 x 8 GB RAM) ### set.seed(20010911) ### system.time(smcp_sandwich <- summary(mcp_sandwich)) # (cca 17 sec. on i5-4570 CPU @ 3.20GH, 2 x 8 GB RAM) ### smcp_sandwich <- summary(mcp_sandwich) smcp smcp_sandwich # std. errors X = model.matrix(mAddit) V = L %*% solve(t(X) %*% X) %*% t(L) D = diag(V) s = summary(mAddit)$sigma s*sqrt(D) smcp$test$sigma all.equal() ### Simultaneous confidence intervals will now be calculated. ### For this, only one quantile of the max-abs-t distribution ### must be calculated (again by Monte Carlo method). ### Again, some time (but less) is needed... ### set.seed(20010911) ### system.time(cimcp <- confint(mcp, level = 0.95)) # (cca 5 sec. on i5-4570 CPU @ 3.20GH, 2 x 8 GB RAM) ### cimcp <- confint(mcp, level = 0.95) ### system.time(cimcp_sandwich <- confint(mcp_sandwich, level = 0.95)) # (cca 5 sec. on i5-4570 CPU @ 3.20GH, 2 x 8 GB RAM) ### cimcp_sandwich <- confint(mcp_sandwich, level = 0.95) ### print(cimcp) ### print(cimcp_sandwich) ### Put both intervals and P-values into one data.frame MCP <- data.frame(Estimate = cimcp[["confint"]][, "Estimate"], Lower = cimcp[["confint"]][, "lwr"], Upper = cimcp[["confint"]][, "upr"], Pvalue = smcp[["test"]][["pvalues"]]) MCP <- MCP[order(MCP[, "Pvalue"]),] round(MCP[MCP[, "Pvalue"] < 0.05,], digits=4) ### For comparison when sandwich estimator is used MCP_sandwich <- data.frame(Estimate = cimcp_sandwich[["confint"]][, "Estimate"], Lower = cimcp_sandwich[["confint"]][, "lwr"], Upper = cimcp_sandwich[["confint"]][, "upr"], Pvalue = smcp_sandwich[["test"]][["pvalues"]]) MCP_sandwich <- MCP_sandwich[order(MCP_sandwich[, "Pvalue"]),] round(MCP_sandwich[MCP_sandwich[, "Pvalue"] < 0.05,], digits=4) round(sort(muAddit[,1]), 3) ### For comparison Tukey HSD round(treg[treg[, "p adj"] < 0.05,], 4) ### Slightly more complicated model, now with [score] as response m11Sum <- lm(score ~ fgender + avg.mark + fspec*mark.eng, data = mana, contrasts = list(fspec = "contr.sum")) summary(m11Sum) drop1(m11Sum, test = "F") ### Specialization is an effect modifier of the relationship between [score] and [mark.eng]. ### ### Which specializations differ with respect to the size of the effect of [mark.eng] on [score]? ### ### => Pairwise comparisons of slopes pertaining to the [mark.eng] regressor for different specializations. ### Used contrast matrix to parameterize categorical covariate "fspec": CZ <- contr.sum(lSPEC) rownames(CZ) <- lSPEC colnames(CZ) <- names(coef(m11Sum))[4:13] print(CZ) ### "L" matrix to calculate [mark.eng] slopes for different specializations Lslope <- cbind(matrix(0, nrow = nrow(CZ), ncol = 13), 1, CZ) rownames(Lslope) <- lSPEC colnames(Lslope) <- names(coef(m11Sum)) print(Lslope) LSest(m11Sum, L = Lslope) (slps <- LSest(m11Sum, L = Lslope)[["Estimate"]]) ### Is there any specialization for which there is no significant relationship ### between [mark.eng] (English at school) and [score] (English at MANA test) ### given gender and [avg.mark]? mean(slps) coef(m11Sum)["mark.eng"] ## the same as mean(slps), surprised? # contrast sum parametrization ### Core of the "L" matrix to calculate all pairwise differences between the slopes. Lcore2 <- matrix(nrow = 0, ncol = ncol(CZ)) colnames(Lcore2) <- colnames(CZ) print(Lcore2) ## matrix with 0 rows rname <- character(0) for (h1 in 1:(length(lSPEC) - 1)){ for (h2 in (h1 + 1):length(lSPEC)){ Lcore2 <- rbind(Lcore2, CZ[h2,] - CZ[h1,]) rname <- c(rname, paste(lSPEC[h2], "-", lSPEC[h1])) } } rownames(Lcore2) <- rname print(Lcore2) ### Complete "L" matrix to calculate all pairwise differences in slopes... ### Can you explain why is the "L" matrix constructed in this way? coef(m11Sum) LslopeDiff <- cbind(matrix(0, nrow = nrow(Lcore2), ncol = 14), Lcore2) colnames(LslopeDiff) <- names(coef(m11Sum)) rownames(LslopeDiff) <- rownames(Lcore2) print(LslopeDiff) ### Unadjusted inference temp <- LSest(m11Sum, LslopeDiff); class(temp) <- c("data.frame") temp <- round(temp,4) temp <- temp[order(temp[["P value"]]),]; temp[temp[["P value"]] <= 0.05,] ### Bretz-Hothorn-Westfall MCP mcpSlp <- glht(m11Sum, linfct = LslopeDiff, rhs = 0) mcpSlp_sandwich <- glht(m11Sum, linfct = LslopeDiff, rhs = 0, vcov=sandwich) ### set.seed(20010911) ### smcpSlp <- summary(mcpSlp) ## 7 sec. ### cimcpSlp <- confint(mcpSlp, level = 0.95) ## 4 sec. ### set.seed(20010911) ### smcpSlp_sandwich <- summary(mcpSlp_sandwich) ## 7 sec. ### cimcpSlp_sandwich <- confint(mcpSlp_sandwich, level = 0.95) ## 4 sec. print(smcpSlp, 4) ### For comparison with sandwich estimator print(smcpSlp_sandwich, 4) print(cimcpSlp, 2) print(cimcpSlp_sandwich, 2) ### When MCP taken into account, are there any two specializations ### with significantly different [mark.eng] slopes? ### save(list = c("smcp", "cimcp", "smcpSlp", "cimcpSlp", "smcp_sandwich", "cimcp_sandwich", "smcpSlp_sandwich", "cimcpSlp_sandwich"), file="nmsa407-mult-comparison.RData") ### QUESTIONS: ### Which specializations differ with respect to the mean [score] ### when [mark.eng] is equal to 2 (median of [mark.eng]) ### while comparing students of the same gender and the same value of [avg.mark]. ### ### => Pairwise comparisons of specializations for a particular value ### of [mark.eng] adjusted for a possible effect of gender and [avg.mark]. median(mana[, "mark.eng"]) ### Complete "L" matrix to calculate all requested pairwise differences ### from the fitted linear model. ### ### Can you explain why is the "L" matrix constructed in this way? coef(m11Sum) mmark.eng <- 2 LspecDiff <- cbind(0, 0, 0, Lcore2, 0, mmark.eng * Lcore2) colnames(LspecDiff) <- names(coef(m11Sum)) rownames(LspecDiff) <- rownames(Lcore2) print(LspecDiff) print(LspecDiff[1:5,]) ### Unadjusted inference temp <- LSest(m11Sum, LspecDiff) class(temp) <- c("data.frame") temp <- round(temp,4) temp <- temp[order(temp[["P value"]]),]; temp[temp[["P value"]] <= 0.05,] ### Bretz-Hothorn-Westfall MCP mcpSpec <- glht(m11Sum, linfct = LspecDiff, rhs = 0) mcpSpec_sandwich <- glht(m11Sum, linfct = LspecDiff, rhs = 0, vcov=sandwich) set.seed(20010911) smcpSpec <- summary(mcpSpec) cimcpSpec <- confint(mcpSpec, level = 0.95) set.seed(20010911) smcpSpec_sandwich <- summary(mcpSpec_sandwich) cimcpSpec_sandwich <- confint(mcpSpec_sandwich, level = 0.95) ### print(smcpSpec) ### print(cimcpSpec) ### Put both intervals and P-values into one data.frame MCPSpec <- data.frame(Estimate = cimcpSpec[["confint"]][, "Estimate"], Lower = cimcpSpec[["confint"]][, "lwr"], Upper = cimcpSpec[["confint"]][, "upr"], Pvalue = smcpSpec[["test"]][["pvalues"]]) MCPSpec <- MCPSpec[order(MCPSpec[, "Pvalue"]),] MCPSpec[, "PvalueF"] <- format.pval(MCPSpec[, "Pvalue"], digits = 2, eps = 1e-4) MCPSpec[MCPSpec[, "Pvalue"] < 0.05, c("Estimate", "Lower", "Upper", "PvalueF")] MCPSpec_sandwich <- data.frame(Estimate = cimcpSpec_sandwich[["confint"]][, "Estimate"], Lower = cimcpSpec_sandwich[["confint"]][, "lwr"], Upper = cimcpSpec_sandwich[["confint"]][, "upr"], Pvalue = smcpSpec_sandwich[["test"]][["pvalues"]]) MCPSpec_sandwich <- MCPSpec_sandwich[order(MCPSpec_sandwich[, "Pvalue"]),] MCPSpec_sandwich[, "PvalueF"] <- format.pval(MCPSpec_sandwich[, "Pvalue"], digits = 2, eps = 1e-4) MCPSpec_sandwich[MCPSpec_sandwich[, "Pvalue"] < 0.05, c("Estimate", "Lower", "Upper", "PvalueF")]