# Linear model fitting with covariance, by group

I need to fit a linear model to a percent change value over time (categorical), for grouped data and need to include age covariance.

I've tried to work with these models:

fit <- lm(Average ~ Event + Age, data=df)

fit2 <- lm(Average ~ Event + Age + group, data=df)


but it doesn't seem to show any difference in the lines slopes between the groups..

Sample data:

structure(list(Event = c(0L, 0L, 0L, 0L, 6L, 12L, 18L, 18L, 0L,
12L, 18L, 6L), group = structure(c(1L, 1L, 2L, 2L, 1L, 3L, 3L,
3L, 3L, 4L, 4L, 1L), .Label = c("A", "B", "C", "D"), class = "factor"),
Age = c(77, 70.2, 69.9, 65.7, 66.2, 66.7, 67.2, 67.7, 66.8,
67.8, 68.3, 68.8), Average = c(0, 5, 0, -2.34, -6, 0, -11,
-3.3, 4.5, 0, -2, 3)), row.names = c(NA, 12L), class = "data.frame")


I ran into the example below in "R-blogger" and was wondering if there's anything like it that would fit my needs

lm(formula = circumference ~ age + Tree + age:Tree, data = orange.df)

• Hi, when you look at the output of summary(fit2), the model fit seems to have improved. with your example data. If you want to create interaction terms, you can do that by Age*group – heck1 Sep 3 '19 at 13:24
• a) I suggest this is a stats question, not a programming question, and it is a better fit at stats.stackexchage. b) Rather than just trying different formulas, you should write out the linear equation you want to estimate. c) I'm not sure what Event is, since it's not in your sampled data. (is it the same as Time?) When you talk about "slope", are you meaning with respect to Age or with respect to Event? Or both? To have another variable influence the slope, you will probably need an interaction term. – Gregor Thomas Sep 3 '19 at 13:26
• Hi @Gregor, thanks for the inputs. The time values are under the 'Event' label. I am trying to fit the 'Average' values over 'Event' (time) – YBB Sep 3 '19 at 13:37
• I suggest to visualize your data before trying to estimate any model. For example, you can use: library(ggplot2) and ggplot(data=df, aes(x=Event, y=Average, color=group)) + geom_line(size=1) + theme_bw(). You can learn a lot from a good plot. – Marco Sandri Sep 3 '19 at 13:37
• Okay, I'd suggest editing your question so that the sample model runs on the sample data. That way, people won't need to read comments and make edits to to run your example. That is, either change Time to Event, or vice versa. And can you respond to my other question: what do you mean by "slope"? coefficient of Age, or coefficient of Time/Event? – Gregor Thomas Sep 3 '19 at 13:41

What you are looking for is an Interaction between your slope variable and the age variable. In R this can be done with the formula syntax, which is somewhat poorly described in the documentation of help(formula), and maybe slightly better described in the comments of my answer here (no promises).

Your formula Average ~ Event + Age + group states "Let average be linearly dependent on event, age and group". For an interaction you would use the :, *, ^ or | descriptors depending on what you want to model.

For your case you likely want Average ~ Event + Age*group which is equivalent to Average ~ Event + Age + group + Age:group, letting each group having a different Intercept and slope. From a statistical point of view this is often also the most sensible, as this allows each group to have a different intercept, and most models don't have a logical reasoning for this not being the case.

Alternatively if age is Nested in group, (each group is a specific age group) you could use Average ~ Event + Age|group with equivalent Average ~ Event + Age + Age:group. This basically translates to removing group dependent intercepts, and often becomes much less reasonable.

For extra confusion, Age*group is also equivalent to (Age+group)^2 in R formula terms. This basically expands the parenthesis to the second order interactions while ^3 would add interactions up to the third order (2 and 3 variable interactions)

# Edit to accommodate comment.

The OP would like to know either the average of Event grouped by group or maybe the %-change of the predicted value. In either case this can be achieved using tapply or a for loop. I will assume the data is called dat for this example.

tapply:

#Average of Average by event given group
tapply(X = dat$Average, INDEX = dat$group, FUN = mean)
#3, 0, 12, 15

#Percentage change (securely)
Percentage_Change <- function(x){
if(!is.numeric(x) || length(x) == 1) return(NA_real_)
rn <- range(x)
#(end - start) / start
diff(rn) / rn
}
fit <- lm(Average ~ Event + Age * group, data = dat)
tapply(X = predict(fit), INDEX = dat$group, FUN = Percentage_Change) #-2.8, -1.0, -1.9, -1.0  This is equivalent to using a for loop as below groups <- unique(dat$group)
nGroups <- length(groups)
predictChange <- eventAverage <- numeric(nGroups)
pred <- predict(fit)
for(i in 1:nGroups){
group <- groups[i]
index <- which(dat\$group == group)
eventAverage[i] <- mean(dat[index , "Event"])
predictChange[i] <- Percentage_Change(pred[index])
}
#Set names
names(predictChange) <- names(eventAverage) <- groups
predictChange
#-2.8, -1.0, -1.9, -1.0
eventAverage
#3, 0, 12, 15

• great answer thank you! Any chance you know how I can get the statics and and percent change from start to finish, to each group? The data is displayed faceted by group – YBB Sep 3 '19 at 13:51
• As in calculate the fitted value at the start and end of each age range, grouped by the grouping factor? – Oliver Sep 3 '19 at 14:09
• The calculated average by event grouped by the 'groups' – YBB Sep 3 '19 at 14:18
• I am still a bit confused, whether you are trying to find the average of the fitted model, the variable Event or the percentage change of either. I will add a short edit to my answer, to illustrate how either could be done. – Oliver Sep 3 '19 at 14:20
• There we go, i hope this helps out. :-) – Oliver Sep 3 '19 at 14:38