CODEWITHSUNDEEP
rlsmeans
RNB
RNB

Reputation: 477

Planned Contrasts using Multcomp in R

I have an analysis i'd like to perform and I have planned out contrasts (not posthoc comparisons!) I would like to make between treatment groups. The treatment group variable has (k =) 4 levels. I plan to make a total of 3 different comparisons, and therefore,—if I understand correctly—I do not need to make any adjustments to the p-values that are calculated since the comparisons are k-1.

I would like to use the multcomp or lsmeans package in R to do this. My question is: Does anyone know if it is possible to do this planned comparison WITHOUT any adjustment made to the confidence intervals (and p-value)? As far as I can tell from the vignettes i've looked at and examples i've seen, the summary.glht() function makes an adjustment as the default and it is not clear to me what option would undo this.

If someone requires a reproducible example, they could use this example that I found on http://www.ats.ucla.edu/stat/r/faq/testing_contrasts.htm:

    library(multcomp)

    hsb <- read.csv("http://www.ats.ucla.edu/stat/data/hsb2.csv")

    m1 <- lm(read ~ socst + factor(ses) * factor(female), data = hsb)
    summary(m1)

    K <- matrix(c(0, 0, 1, -1, 0, 0, 0), 1)
    t <- glht(m1, linfct = K)
    summary(t)

Upvotes: 1

Views: 1590

Answers (1)

cuttlefish44
cuttlefish44

Reputation: 6786

As far as I see your example, your question is a little strange. At least IMO, If you needn't adjust, you needn't use multcomp package (but in some situation, it saves us some time).

library(multcomp)
hsb <- read.csv("http://www.ats.ucla.edu/stat/data/hsb2.csv")
hsb$ses <- as.factor(hsb$ses)

m3 <- lm(read ~ socst + ses, data = hsb)
l3 <- glht(m3, linfct = mcp(ses = "Tukey"))

# mcp(~) doesn't run with some type of model. If so, you'll give the matrix directly.
# k3 <- matrix(c(0, 0, 1, 0,
#                0, 0, 0, 1,
#                0, 0, -1, 1), byrow = T, ncol = 4)
# rownames(k3) <- c("2-1", "3-1", "3-2")
# l3 <- glht(m3, linfct = k1)

summary(l3, test=adjusted("none"))    # this is the result without adjustment
#            Estimate Std. Error t value Pr(>|t|)
# 2 - 1 == 0   0.6531     1.4562   0.448    0.654
# 3 - 1 == 0   2.7034     1.6697   1.619    0.107
# 3 - 2 == 0   2.0503     1.3685   1.498    0.136

hsb$ses <- relevel(hsb$ses, ref="2")  # change of the order of levels
m3.2 <- lm(read ~ socst + ses, data = hsb)
summary(m3)     # "Without adjustment" means it's equivalent to original model's statistics.
#             Estimate Std. Error t value Pr(>|t|)
#  :
# ses2         0.65309    1.45624   0.448    0.654    
# ses3         2.70342    1.66973   1.619    0.107 

summary(m3.2)
#             Estimate Std. Error t value Pr(>|t|)
#  :
# ses3         2.05033    1.36846   1.498    0.136

# When argument is lmer.obj, summary(~, adjusted("none")) returns p.value by using z value with N(0, 1).

Upvotes: 1

Related Questions