R中的生长曲线分析——两条生长曲线的比较
Growth curve analysis in R - comparison of two growth curves
我有一个小任务要做,不幸的是我对统计这个领域不熟悉...实际上我确实需要计算(我不是在寻找现成的解决方案),但我不知道是否他们是正确的,也是我的思维方式,因此如果你看一下并指出我的错误,我将不胜感激。
这里是假数据,呈现狗和猫的增长率(完全虚构):
time <- c(1,2,3,4,5,6,7,8,9,10,1,2,3,4,5,6,7,8,9,10)
a <- rep('dog', 10)
b <- rep('cat', 10)
animal <- c(a,b)
val <- c(2.00,8.00,17.00,21.00,29.00,37.00,41.00,56.00,67.00,82.00,1.00,3.00,6.00,8.00,11.00,15.00,21.00,26.00,31.00,37.00)
data <- data.frame(time,animal,val)
仔细看看:
require(ggplot2)
ggplot(data, aes(time, val, color=animal)) +
stat_summary(fun.data=mean_se, geom="pointrange") +
geom_point()
如您所见,狗比猫长得快 - 这可能是我的假设。但是我需要做一些统计来符合它。
所以我决定进行 生长曲线分析 (GCA)。
我是基于 this 教程。在我的结果下面有简要的解释。
所以首先我做了一个基础模型,为每只动物随机截取:
m.base <- lmer(val ~ time + (1 | animal), data=data, REML = F)
这里我遇到了问题,其实我这里没有任何固定效应,我的数据集很简单,我只想知道我的两组(狗和猫)的时间增长率在统计上显着不同。
换一种说法。
动物在这段时间内的生长速度是否不同?
因此我把我的动物作为一个额外的固定效果:
m.1 <- lmer(val ~ time * animal + (1 | animal), data=data, REML = F)
现在,为了检查是否存在统计显着差异,我使用方差分析比较了两个模型。
> anova(m.base,m.1)
Data: data
Models:
m.base: val ~ time + (1 | animal)
m.1: val ~ time * animal + (1 | animal)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
m.base 4 151.43 155.41 -71.714 143.43
m.1 6 116.29 122.26 -52.145 104.29 39.138 2 3.171e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
现在我很困惑,我不完全理解所有这些分析,几个问题...
这个值3.171e-09表示我组的增长率在统计上有显着差异吗?
要不要我再做一个模型:
m.0 <- lmer(val ~ time + animal + (1 | animal), data=data, REML = F)
然后进行模型测试?
> anova(m.base,m.0,m.1)
Data: data
Models:
m.base: val ~ time + (1 | animal)
m.0: val ~ time + animal + (1 | animal)
m.1: val ~ time * animal + (1 | animal)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
m.base 4 151.43 155.41 -71.714 143.43
m.0 5 145.58 150.56 -67.789 135.58 7.8499 1 0.005082 **
m.1 6 116.29 122.26 -52.145 104.29 31.2884 1 2.224e-08 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
哪个值可以让我证实我的假设?
非常感谢指出我的错误,任何线索和解释!
我假设每个数据点都来自不同的动物。如果你只有两只动物的数据,你只能比较这两种动物,而不能推断出关于这两个种群的任何信息。如果你有几只动物的数据,但每只动物都被重复测量,你确实需要一个混合效应模型。但是根据我上面的假设,你不需要它。
您现在可以使用来自特定领域理论的参数模型并使用 nlme::gnls
。该函数基本上适合非线性模型,其中参数是一些其他变量的线性模型(在您的情况下是动物的类型)。然后可以测试这些线性模型的参数的显着性,summary
方法可以为您完成。如果您有重复测量,nlme::nlme
将其扩展到混合效应模型。
另一种方法是非参数模型:
library(mgcv)
mod1 <- gam(val ~ s(time, k = 4), data = data, select = TRUE)
mod2 <- gam(val ~ animal + s(time, k = 4, by = animal), data = data, select = TRUE)
#we need the parametric effect because smoothers are centered
#compare both models, not sure which test is more appropriate,
#let's just do both Chisq and F
anova(mod1, mod2, test = "Chisq")
anova(mod1, mod2, test = "F")
#significant difference between animal types
#plots show which one grows faster
gam.check(mod2)
plot(mod2)
summary(mod2)
如有必要,这也可以扩展到混合效应模型。
我有一个小任务要做,不幸的是我对统计这个领域不熟悉...实际上我确实需要计算(我不是在寻找现成的解决方案),但我不知道是否他们是正确的,也是我的思维方式,因此如果你看一下并指出我的错误,我将不胜感激。
这里是假数据,呈现狗和猫的增长率(完全虚构):
time <- c(1,2,3,4,5,6,7,8,9,10,1,2,3,4,5,6,7,8,9,10)
a <- rep('dog', 10)
b <- rep('cat', 10)
animal <- c(a,b)
val <- c(2.00,8.00,17.00,21.00,29.00,37.00,41.00,56.00,67.00,82.00,1.00,3.00,6.00,8.00,11.00,15.00,21.00,26.00,31.00,37.00)
data <- data.frame(time,animal,val)
仔细看看:
require(ggplot2)
ggplot(data, aes(time, val, color=animal)) +
stat_summary(fun.data=mean_se, geom="pointrange") +
geom_point()
如您所见,狗比猫长得快 - 这可能是我的假设。但是我需要做一些统计来符合它。
所以我决定进行 生长曲线分析 (GCA)。 我是基于 this 教程。在我的结果下面有简要的解释。
所以首先我做了一个基础模型,为每只动物随机截取:
m.base <- lmer(val ~ time + (1 | animal), data=data, REML = F)
这里我遇到了问题,其实我这里没有任何固定效应,我的数据集很简单,我只想知道我的两组(狗和猫)的时间增长率在统计上显着不同。 换一种说法。 动物在这段时间内的生长速度是否不同?
因此我把我的动物作为一个额外的固定效果:
m.1 <- lmer(val ~ time * animal + (1 | animal), data=data, REML = F)
现在,为了检查是否存在统计显着差异,我使用方差分析比较了两个模型。
> anova(m.base,m.1)
Data: data
Models:
m.base: val ~ time + (1 | animal)
m.1: val ~ time * animal + (1 | animal)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
m.base 4 151.43 155.41 -71.714 143.43
m.1 6 116.29 122.26 -52.145 104.29 39.138 2 3.171e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
现在我很困惑,我不完全理解所有这些分析,几个问题...
这个值3.171e-09表示我组的增长率在统计上有显着差异吗?
要不要我再做一个模型:
m.0 <- lmer(val ~ time + animal + (1 | animal), data=data, REML = F)
然后进行模型测试?
> anova(m.base,m.0,m.1)
Data: data
Models:
m.base: val ~ time + (1 | animal)
m.0: val ~ time + animal + (1 | animal)
m.1: val ~ time * animal + (1 | animal)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
m.base 4 151.43 155.41 -71.714 143.43
m.0 5 145.58 150.56 -67.789 135.58 7.8499 1 0.005082 **
m.1 6 116.29 122.26 -52.145 104.29 31.2884 1 2.224e-08 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
哪个值可以让我证实我的假设?
非常感谢指出我的错误,任何线索和解释!
我假设每个数据点都来自不同的动物。如果你只有两只动物的数据,你只能比较这两种动物,而不能推断出关于这两个种群的任何信息。如果你有几只动物的数据,但每只动物都被重复测量,你确实需要一个混合效应模型。但是根据我上面的假设,你不需要它。
您现在可以使用来自特定领域理论的参数模型并使用 nlme::gnls
。该函数基本上适合非线性模型,其中参数是一些其他变量的线性模型(在您的情况下是动物的类型)。然后可以测试这些线性模型的参数的显着性,summary
方法可以为您完成。如果您有重复测量,nlme::nlme
将其扩展到混合效应模型。
另一种方法是非参数模型:
library(mgcv)
mod1 <- gam(val ~ s(time, k = 4), data = data, select = TRUE)
mod2 <- gam(val ~ animal + s(time, k = 4, by = animal), data = data, select = TRUE)
#we need the parametric effect because smoothers are centered
#compare both models, not sure which test is more appropriate,
#let's just do both Chisq and F
anova(mod1, mod2, test = "Chisq")
anova(mod1, mod2, test = "F")
#significant difference between animal types
#plots show which one grows faster
gam.check(mod2)
plot(mod2)
summary(mod2)
如有必要,这也可以扩展到混合效应模型。