CODEWITHSUNDEEP
r
Mac
Mac

Reputation: 367

Testing differences in coefficients including interactions from piecewise linear model

I'm running a piecewise linear random coefficient model testing the influence of a covariate on the second piece. Thereby, I want to test whether the coefficient of the second piece under the influence of the covariate (piece2 + piece2:covariate) differs from the coefficient of the first piece (piece1), hence whether the growth rate differs.

I set up some exemplary data:

set.seed(100)

# set up dependent variable   
temp <- rep(seq(0,23),50)
y <- c(rep(seq(0,23),50)+rnorm(24*50), ifelse(temp <= 11, temp + runif(1200), temp + rnorm(1200) + (temp/sqrt(temp))))

# set up ID variable, variables indicating pieces and the covariate
id <- sort(rep(seq(1,100),24))
piece1 <- rep(c(seq(0,11), rep(11,12)),100)
piece2 <- rep(c(rep(0,12), seq(1,12)),100)
covariate <- c(rep(0,24*50), rep(c(rep(0,12), rep(1,12)), 50))

# data frame
example.data <- data.frame(id, y, piece1, piece2, covariate)

# run piecewise linear random effects model and show results
library(lme4)
lmer.results <- lmer(y ~ piece1 + piece2*covariate + (1|id) , example.data)  
summary(lmer.results)

I came across the linearHypothesis() command from the car package to test differences in coefficients. However, I could not find an example on how to use it when including interactions.

Can I even use linearHypothesis() to test this or am I aiming for the wrong test?

I appreciate your help. Many thanks in advance! Mac

Upvotes: 2

Views: 800

Answers (1)

Dieter Menne
Dieter Menne

Reputation: 10215

Assuming your output looks like this

                 Estimate Std. Error t value
(Intercept)       0.26293    0.04997     5.3
piece1            0.99582    0.00677   147.2
piece2            0.98083    0.00716   137.0
covariate         2.98265    0.09042    33.0
piece2:covariate  0.15287    0.01286    11.9

if I understand correctly what you want, you are looking for the contrast: piece1-(piece2+piece2:covariate)

or

c(0,1,-1,0,-1)

My preferred tool for this is function estimable in gmodels; you could also do it by hand or with one of the functions in Frank Harrel's packages.

library(gmodels)
estimable(lmer.results,c(0,1,-1,0,-1),conf.int=TRUE)

giving 

              Estimate Std. Error p value Lower.CI Upper.CI
(0 1 -1 0 -1)   -0.138     0.0127       0   -0.182  -0.0928

Upvotes: 1

Related Questions