This page might help you:

https://stackoverflow.com/questions/45227690/how-to-run-post-hoc-test-pairwise-comparison-with-lsmeans-on-a-polynomial-lmer

Other techniques are Turkey's test and Duncan multiple range test.

If I understand the question correctly, this is a good question.

To clarify: You are conducting an analysis of covariance (ancova) with more than two lines. Assuming equal slopes, it is straightforward to compare among the intercepts of the lines. However, it is not clear how you would compare among the slopes of the multiple lines when the slopes differ.

As a general comment, I would move away from using functions like HSD.test or scheffe.test, and instead transition to glht or emmeans (née lsmeans). These are more flexible, contemporary, and tolerant of irregularities in the data. emmeans is useful because it can make comparisons among interaction terms.

I find interpreting ancova with more than two lines confusing. I will try to include another answer later with a different approach.

The following example creates some data, conducts an ancova, and uses glht to compare among the intercepts. (I think that far is good.) From there, it simply subsets the data, and creates new models for each pair of lines. I don't particularly like this approach, but it might be as good as any.

I would love to hear from those better versed than me in statistics or R programming to find a more rigorous solution.

### --------------------------------------------------------------

### Analysis of covariance, hypothetical data

### --------------------------------------------------------------

if(!require(car)){install.packages("car")}

if(!require(multcomp)){install.packages("multcomp")}

if(!require(ggplot2)){install.packages("ggplot2")}

X = 0:19

Group1 = 0 + X + rnorm(20,0,1)

Group2 = 0 + 2*X + rnorm(20,0,1)

Group3 = 10 + X + rnorm(20,0,1)

Y = c(Group1, Group2, Group3)

Group = c(rep("I", 20), rep("II", 20), rep("III", 20))

X = rep(X,3)

Data = data.frame(Group, X, Y)

library(ggplot2)

ggplot(data = Data,

aes(x = X,

y = Y,

col = Group)) +

geom_point() +

geom_smooth(method=lm) +

scale_color_manual(values = c("red", "darkgreen", "blue"))

model = lm (Y ~ X + Group + X:Group,

data = Data)

library(car)

Anova(model, type="II")

library(multcomp)

mc = glht(model,

mcp(Group = "Tukey"))

summary(mc, test=adjusted("single-step"))

### Above is a test of intercepts, at least assuming equal slopes

##############################

### Now, subset data and make new models

Data.1 = Data[Data$Group=="I" | Data$Group=="II",]

model.1 = lm (Y ~ X + Group + X:Group,

data = Data.1)

Anova(model.1, type="II")

ggplot(data = Data.1,

aes(x = X,

y = Y,

col = Group)) +

geom_point() +

geom_smooth(method=lm) +

scale_color_manual(values = c("red", "darkgreen"))

# # #

Data.2 = Data[Data$Group=="I" | Data$Group=="III",]

model.2 = lm (Y ~ X + Group + X:Group,

data = Data.2)

Anova(model.2, type="II")

ggplot(data = Data.2,

aes(x = X,

y = Y,

col = Group)) +

geom_point() +

geom_smooth(method=lm) +

scale_color_manual(values = c("red", "blue"))

# # #

Data.3 = Data[Data$Group=="II" | Data$Group=="III",]

model.3 = lm (Y ~ X + Group + X:Group,

data = Data.3)

Anova(model.3, type="II")

ggplot(data = Data.3,

aes(x = X,

y = Y,

col = Group)) +

geom_point() +

geom_smooth(method=lm) +

scale_color_manual(values = c("darkgreen", "blue"))

##########################

I think that if the idea of the analysis is to compare slopes and intercepts of linear models, that it makes sense to compare these directly, and not worry so much about the usual ancova framework. To me, the point of ancova is really to see if there is a difference among the categorical treatments or groups, and use the continuous variable as a covariate.

I'm sure there are ways to conduct tests to compare among the intercepts or the slopes of linear models. But here, I'll just propose to construct 95% confidence intervals of the intercepts and slopes, and consider those statistics without overlapping confidence intervals to be significantly different.

The following example uses the same data as my previous example, and uses bootstrap for the 95% confidence intervals. A summary of results are at the bottom of the code.

EDIT: NOTE: because I didn't use set.seed, and there are random values in the data, your results may differ somewhat from mine.

### --------------------------------------------------------------

### Hypothetical data

### --------------------------------------------------------------

if(!require(ggplot2)){install.packages("ggplot2")}

if(!require(boot)){install.packages("boot")}

X = 0:19

Group1 = 0 + X + rnorm(20,0,1)

Group2 = 0 + 2*X + rnorm(20,0,1)

Group3 = 10 + X + rnorm(20,0,1)

Y = c(Group1, Group2, Group3)

Group = c(rep("I", 20), rep("II", 20), rep("III", 20))

X = rep(X,3)

Data = data.frame(Group, X, Y)

library(ggplot2)

ggplot(data = Data,

aes(x = X,

y = Y,

col = Group)) +

geom_point() +

geom_smooth(method=lm) +

scale_color_manual(values = c("red", "darkgreen", "blue"))

############# Intercept #############

Data.i = Data[Data$Group=="I",]

model.i = lm (Y ~ X,

data = Data.i)

summary(model.i)

library(boot)

Function = function(input, index){

Input = input[index,]

Result = lm (Y ~ X,

data = Input)

Stat = summary(Result)$coefficients[1,1]

return(Stat)}

Boot = boot(Data.i,

Function,

R=5000)

boot.ci(Boot,

conf = 0.95,

type = "perc")

############# Slope, coefficient for X #############

Function = function(input, index){

Input = input[index,]

Result = lm (Y ~ X,

data = Input)

Stat = summary(Result)$coefficients[2,1]

return(Stat)}

Boot = boot(Data.i,

Function,

R=5000)

boot.ci(Boot,

conf = 0.95,

type = "perc")

####################################

### Results:

### Model Intercept Lower.CI Upper.CI Separation.letter

### I 0.04919 -1.0195 0.6936 a

### II 0.4130 -0.1723 0.8746 a

### III 9.94016 8.816 11.053 b

###

### Model Slope Lower.CI Upper.CI Separation.letter

### I 0.99337 0.9336 1.0778 a

### II 1.9891 1.926 2.052 b

### III 1.02730 0.936 1.112 a

Dear Salvatore S. Mangiafico,

1. I'm glad you understood what i was asking and agree with this question. Over two days of internet search, I am surprised by the lack of information and statisticians addressing this situation. I guess many merely follow the ANCOVA framework instead of really getting into the details - what each step mean and address. I could be wrong though, - it could be that post-hoc slope comparison is redundant for some reason.

2. Once again, thank you very much for the elaborate answers, explanation, and example script! Using confidence interval is good and creative idea. I just briefly looked through the script and will put it to test later!

And you are being too modest! Your R Handbook is a gem!

By the way, Salvatore,

"From there, it simply subsets the data, and creates new models for each pair of lines... "

I may have misunderstood, but are you saying that the post-hoc creates new models for each factor? I thought it uses the same model produced in ANCOVA but merely compares them using an adjusted p-value, or adjusted r-squared.

Thanks...

"From there, it simply subsets the data, and creates new models for each pair of lines... "

The approach is to create pairwise data frames, each with only two lines, and conduct an ancova on each of these data frames.

So, the first data frame (Data.1) includes only Groups I and II. The anova table from this model addresses only these two models, so the p-value for the interaction tells you if this pair of lines has different slopes. The next data frame (Data.2) includes only Groups I and III, so the p-value for the interaction tells you if this pair of lines has different slopes. And so on.

You might consider the compSlopes() function in my FSA package. It basically does each pairwise slope comparison (as suggested by Salvatore above, but in a different way) and then applies an adjustment to the p-values for the problem of multiple comparisons. It takes the results of the "ANCOVA" model (i.e., with the covariate, the factor, and the covariate and factor interaction) as input and returns a fairly simple (in my opinion) table of results as output.

See the documentation for compSlopes() at https://www.rdocumentation.org/packages/FSA/versions/0.8.17/topics/compSlopes

Note that this function was my solution to this fairly common problem in my line of work as I could not find a better solution. However, I am open to hearing better solutions.

Thank you Waldemar Koczkodaj and Derek H. Ogle!

Will try out those methods!