We find that the residual standard deviations of the groups that are not group A are about double the residual standard deviation for group A. And that there is negative autocorrelation - positive-negative residual switching pattern by time.
Heteroskedastic Jim
- 5.2k
- 1
- 17
- 40
Heteroskedastic Jim
- 5.2k
- 1
- 17
- 40
However, it appears Group A is on a different scale, so it makes sense to run a heteroskedastic model. Additionally, as @whuber pointed out in the comments, it makes sense to model the autocorrelation. I use the simple autocorrelation of order 1:
library(nlme)
(m1 <- gls(Size ~ Time * Group, dat, correlation = corAR1(form = ~ Time | Group),
weights = varIdent(form = ~ 1 | I(Group == "A"))))
# Coefficients:
# (Intercept) Time GroupB GroupC Time:GroupB Time:GroupC
# 0.17000001768985 0.39300003900313 -0.14821431543012 -0.10321431146352 -0.28895242860587 -0.29561902912242
#
# Correlation Structure: AR(1)
# Formula: ~Time | Group
# Parameter estimate(s):
# Phi
# -0.5295663
# Variance function:
# Structure: Different standard deviations per stratum
# Formula: ~1 | I(Group == "A")
# Parameter estimates:
# TRUE FALSE
# 1.000000 12.676464152732
We find that the residual standard deviations of the groups that are not group A are 68% greater thanabout double the residual standard deviation for group A. And that there is negative autocorrelation
library(emmeans)
pairs(emtrends(m1, ~ Group, var = "Time"))
# contrast estimate SE df t.ratio p.value
# A - B 0.288952381286058700 0.009169459005098842 14 3156.512103 <.0001
# A - C 0.295619048291224187 0.009169459005098842 14 3257.240116 <.0001
# B - C 0.006666667005165488 0.006492917003857697 14 1.027339 0.57293979
#
# P value adjustment: tukey method for comparing a family of 3 estimates
m.pars <- gls(Size ~ Time * I(Group == "A"), dat,
correlation = corAR1(form = ~ Time | Group),
weights = varIdent(form = ~ 1 | I(Group == "A")))
However, it appears Group A is on a different scale, so it makes sense to run a heteroskedastic model:
library(nlme)
(m1 <- gls(Size ~ Time * Group, dat,
weights = varIdent(form = ~ 1 | I(Group == "A"))))
# Coefficients:
# (Intercept) Time GroupB GroupC Time:GroupB Time:GroupC
# 0.1700000 0.3930000 -0.1482143 -0.1032143 -0.2889524 -0.2956190
#
# Variance function:
# Structure: Different standard deviations per stratum
# Formula: ~1 | I(Group == "A")
# Parameter estimates:
# TRUE FALSE
# 1.000000 1.676464
We find that the residual standard deviations of the groups that are not group A are 68% greater than the residual standard deviation for group A.
library(emmeans)
pairs(emtrends(m1, ~ Group, var = "Time"))
# contrast estimate SE df t.ratio p.value
# A - B 0.288952381 0.009169459 14 31.512 <.0001
# A - C 0.295619048 0.009169459 14 32.240 <.0001
# B - C 0.006666667 0.006492917 14 1.027 0.5729
#
# P value adjustment: tukey method for comparing a family of 3 estimates
m.pars <- gls(Size ~ Time * I(Group == "A"), dat,
weights = varIdent(form = ~ 1 | I(Group == "A")))
However, it appears Group A is on a different scale, so it makes sense to run a heteroskedastic model. Additionally, as @whuber pointed out in the comments, it makes sense to model the autocorrelation. I use the simple autocorrelation of order 1:
library(nlme)
(m1 <- gls(Size ~ Time * Group, dat, correlation = corAR1(form = ~ Time | Group),
weights = varIdent(form = ~ 1 | I(Group == "A"))))
# Coefficients:
# (Intercept) Time GroupB GroupC Time:GroupB Time:GroupC
# 0.1768985 0.3900313 -0.1543012 -0.1146352 -0.2860587 -0.2912242
#
# Correlation Structure: AR(1)
# Formula: ~Time | Group
# Parameter estimate(s):
# Phi
# -0.5295663
# Variance function:
# Structure: Different standard deviations per stratum
# Formula: ~1 | I(Group == "A")
# Parameter estimates:
# TRUE FALSE
# 1.000000 2.152732
We find that the residual standard deviations of the groups that are not group A are about double the residual standard deviation for group A. And that there is negative autocorrelation
library(emmeans)
pairs(emtrends(m1, ~ Group, var = "Time"))
# contrast estimate SE df t.ratio p.value
# A - B 0.286058700 0.005098842 14 56.103 <.0001
# A - C 0.291224187 0.005098842 14 57.116 <.0001
# B - C 0.005165488 0.003857697 14 1.339 0.3979
#
# P value adjustment: tukey method for comparing a family of 3 estimates
m.pars <- gls(Size ~ Time * I(Group == "A"), dat,
correlation = corAR1(form = ~ Time | Group),
weights = varIdent(form = ~ 1 | I(Group == "A")))
This is a relatively simple problem. The basic model to test your question about differences in slope is:
(m0 <- lm(Size ~ Time * Group, dat))
# Coefficients:
# (Intercept) Time GroupB GroupC Time:GroupB Time:GroupC
# 0.1700 0.3930 -0.1482 -0.1032 -0.2890 -0.2956
I have ignored the question about the intercepts. More on this at the end. Also, the basic model you ran does not permit testing of differences in slopes. If you perform the diagnostic tests you performed on the model m0 here, they do not confirm misspecification.
However, it appears Group A is on a different scale, so it makes sense to run a heteroskedastic model:
library(nlme)
(m1 <- gls(Size ~ Time * Group, dat,
weights = varIdent(form = ~ 1 | I(Group == "A"))))
# Coefficients:
# (Intercept) Time GroupB GroupC Time:GroupB Time:GroupC
# 0.1700000 0.3930000 -0.1482143 -0.1032143 -0.2889524 -0.2956190
#
# Variance function:
# Structure: Different standard deviations per stratum
# Formula: ~1 | I(Group == "A")
# Parameter estimates:
# TRUE FALSE
# 1.000000 1.676464
We find that the residual standard deviations of the groups that are not group A are 68% greater than the residual standard deviation for group A.
To address your primary research questions, we can go:
library(emmeans)
pairs(emtrends(m1, ~ Group, var = "Time"))
# contrast estimate SE df t.ratio p.value
# A - B 0.288952381 0.009169459 14 31.512 <.0001
# A - C 0.295619048 0.009169459 14 32.240 <.0001
# B - C 0.006666667 0.006492917 14 1.027 0.5729
#
# P value adjustment: tukey method for comparing a family of 3 estimates
We find that there is not much statistical evidence to conclude that the slopes for Group B and C are different from each other. While there is the evidence to differentiate A from B, and A from C.
Since we have an interaction, it is difficult to consider differences in the intercept. Given the current analysis, the intercept relates to group differences at Time 0 which does not exist in the data, minimum Time is 1. The emmeans package provide an option to view differences between the groups at different values of time:
emmip(m1, Time ~ Group, cov.reduce = FALSE)
We find that as time increases, the group differences between A and B, and A and C increase. But B and C continue to be relatively similar. Be careful because there are no time point beyond time 4 for Group A, these are extrapolated values.
Given what we have learned, a parsimonous model would be:
m.pars <- gls(Size ~ Time * I(Group == "A"), dat,
weights = varIdent(form = ~ 1 | I(Group == "A")))