线性混合模型（NLME 或 LMER）- 在 coefficients/estimates 上设置界限
Linear Mixed Model (NLME or LMER) - Putting bounds on coefficients/estimates
我正在尝试在混合效果模型上建立放置界限。我通常为我的模型使用 LMER 函数，但我找不到任何方法来限制系数。我尝试使用 LME，但即使这样也无济于事。有人可以帮忙吗？
library(nlme)
library(lmerTest)
myDat = structure(list(Score = c(1.62, 2.18, 2.3, 3.46, 3.85, 4.7, 1.41,
2.21, 3.32, 2.73, 3.34, 3.27, 2.14, 2.73, 2.74, 3.39, 3.59, 4.01,
1.81, 1.83, 3.22, 3.64, 3.51, 4.26), Subject = structure(c(1L,
1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L, 5L, 5L, 5L, 6L, 6L,
6L, 7L, 7L, 7L, 8L, 8L, 8L), .Label = c("A", "B", "C", "D", "E",
"F", "G", "H"), class = "factor"), Condition = structure(c(1L,
1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L,
2L, 1L, 1L, 1L, 2L, 2L, 2L), .Label = c("No", "Yes"), class = "factor"),
Time = structure(c(1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L,
2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L), .Label = c("1PM",
"2PM", "3PM"), class = "factor")), .Names = c("Score", "Subject",
"Condition", "Time"), class = "data.frame", row.names = c(NA,
-24L))
m1 <- lmer(Score ~ Condition + Time + Condition*Time + (1 | Subject),
data = myDat)
fixef(m1)
ranef(m1)
# This is a lame attempt to get positive coeffiecients but I am not sure what objective function, I need to spcify
m2 <- lme(Score ~ Condition + Time + Condition*Time,
data = myDat, random=c(~1 | Subject)
# ,
# nlminb(start = c(0,0,0,0,0),lower = c(0,0,0,0,0))
)
fixef(m2)
ranef(m2)
如果您只想接受正参数的可能性，我认为 Ben Bolker 已经写过关于在 lme4 中实现它的文章。
虽然，在贝叶斯方法中，您可以根据需要的支持将先验放在参数上。下面的示例使用来自 brms 包（具有类似的公式符号）的下限参数 lb 来编码折叠法线先验。当然，在正确的支持下还有很多可能的先验。
library(lme4)
library(brms)
myDat <- structure(list(Score = c(1.62, 2.18, 2.3, 3.46, 3.85, 4.7, 1.41,
2.21, 3.32, 2.73, 3.34, 3.27, 2.14, 2.73, 2.74, 3.39, 3.59, 4.01,
1.81, 1.83, 3.22, 3.64, 3.51, 4.26), Subject = structure(c(1L,
1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L, 5L, 5L, 5L, 6L, 6L,
6L, 7L, 7L, 7L, 8L, 8L, 8L), .Label = c("A", "B", "C", "D", "E",
"F", "G", "H"), class = "factor"), Condition = structure(c(1L,
1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L,
2L, 1L, 1L, 1L, 2L, 2L, 2L), .Label = c("No", "Yes"), class = "factor"),
Time = structure(c(1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L,
2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L), .Label = c("1PM",
"2PM", "3PM"), class = "factor")), .Names = c("Score", "Subject",
"Condition", "Time"), class = "data.frame", row.names = c(NA,
## your first try -24L))
m1 <- lmer(Score ~ Condition + Time + Condition*Time + (1 | Subject),
data = myDat)
summary(m1)
#
# Linear mixed model fit by REML ['lmerMod']
# Formula: Score ~ Condition + Time + Condition * Time + (1 | Subject)
# Data: myDat
#
# REML criterion at convergence: 25.7
#
# Scaled residuals:
# Min 1Q Median 3Q Max
# -1.4009 -0.5397 0.1392 0.4146 1.2662
# Random effects:
# Groups Name Variance Std.Dev.
# Subject (Intercept) 0.05619 0.2370
# Residual 0.11329 0.3366
# Number of obs: 24, groups: Subject, 8
#
# Fixed effects:
# Estimate Std. Error t value
# (Intercept) 1.7450 0.2058 8.477
# ConditionYes 1.5600 0.2911 5.359
# Time2PM 0.4925 0.2380 2.069
# Time3PM 1.1500 0.2380 4.832
# ConditionYes:Time2PM -0.2250 0.3366 -0.668
# ConditionYes:Time3PM -0.3950 0.3366 -1.174
#
# Correlation of Fixed Effects:
# (Intr) CndtnY Tim2PM Tim3PM CY:T2P
# ConditionYs -0.707
# Time2PM -0.578 0.409
# Time3PM -0.578 0.409 0.500
# CndtnY:T2PM 0.409 -0.578 -0.707 -0.354
# CndtnY:T3PM 0.409 -0.578 -0.354 -0.707 0.500
## brms - get priors
get_prior(Score ~ Condition + Time + Condition*Time + (1 | Subject),
data = myDat, family= gaussian())
prior_positive <- c(
prior(normal(0, 1), class = Intercept),
prior(normal(0, 1), class = b, lb = 0), # specify lower bound with "lb"
#prior(cauchy(0,1), class = shape),
prior(exponential(1), class = sd),
prior(exponential(1), class = sigma))
mod1_positive <- brm(Score ~ Condition + Time + Condition*Time + (1 | Subject),
data = myDat, family= gaussian(), prior = prior_positive)
summary(mod1_positive)
# Family: gaussian
# Links: mu = identity; sigma = identity
# Formula: Score ~ Condition + Time + Condition * Time + (1 | Subject)
# Data: myDat (Number of observations: 24)
# Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
# total post-warmup samples = 4000
#
# Group-Level Effects:
# ~Subject (Number of levels: 8)
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept) 0.26 0.16 0.02 0.65 1.00 954 1653
#
# Population-Level Effects:
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# Intercept 1.91 0.21 1.48 2.32 1.00 1645 2280
# ConditionYes 1.11 0.27 0.55 1.65 1.00 1431 1742
# Time2PM 0.27 0.17 0.02 0.65 1.00 1861 1176
# Time3PM 0.82 0.22 0.37 1.23 1.00 1793 1129
# ConditionYes:Time2PM 0.28 0.21 0.01 0.76 1.00 2183 1490
# ConditionYes:Time3PM 0.24 0.20 0.01 0.75 1.00 1624 1475
#
# Family Specific Parameters:
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sigma 0.40 0.08 0.27 0.58 1.00 1379 2100
#
# Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
# and Tail_ESS are effective sample size measures, and Rhat is the potential
# scale reduction factor on split chains (at convergence, Rhat = 1).
我正在尝试在混合效果模型上建立放置界限。我通常为我的模型使用 LMER 函数，但我找不到任何方法来限制系数。我尝试使用 LME，但即使这样也无济于事。有人可以帮忙吗？
library(nlme)
library(lmerTest)
myDat = structure(list(Score = c(1.62, 2.18, 2.3, 3.46, 3.85, 4.7, 1.41,
2.21, 3.32, 2.73, 3.34, 3.27, 2.14, 2.73, 2.74, 3.39, 3.59, 4.01,
1.81, 1.83, 3.22, 3.64, 3.51, 4.26), Subject = structure(c(1L,
1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L, 5L, 5L, 5L, 6L, 6L,
6L, 7L, 7L, 7L, 8L, 8L, 8L), .Label = c("A", "B", "C", "D", "E",
"F", "G", "H"), class = "factor"), Condition = structure(c(1L,
1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L,
2L, 1L, 1L, 1L, 2L, 2L, 2L), .Label = c("No", "Yes"), class = "factor"),
Time = structure(c(1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L,
2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L), .Label = c("1PM",
"2PM", "3PM"), class = "factor")), .Names = c("Score", "Subject",
"Condition", "Time"), class = "data.frame", row.names = c(NA,
-24L))
m1 <- lmer(Score ~ Condition + Time + Condition*Time + (1 | Subject),
data = myDat)
fixef(m1)
ranef(m1)
# This is a lame attempt to get positive coeffiecients but I am not sure what objective function, I need to spcify
m2 <- lme(Score ~ Condition + Time + Condition*Time,
data = myDat, random=c(~1 | Subject)
# ,
# nlminb(start = c(0,0,0,0,0),lower = c(0,0,0,0,0))
)
fixef(m2)
ranef(m2)
如果您只想接受正参数的可能性，我认为 Ben Bolker 已经写过关于在 lme4 中实现它的文章。
虽然，在贝叶斯方法中，您可以根据需要的支持将先验放在参数上。下面的示例使用来自 brms 包（具有类似的公式符号）的下限参数 lb 来编码折叠法线先验。当然，在正确的支持下还有很多可能的先验。
library(lme4)
library(brms)
myDat <- structure(list(Score = c(1.62, 2.18, 2.3, 3.46, 3.85, 4.7, 1.41,
2.21, 3.32, 2.73, 3.34, 3.27, 2.14, 2.73, 2.74, 3.39, 3.59, 4.01,
1.81, 1.83, 3.22, 3.64, 3.51, 4.26), Subject = structure(c(1L,
1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L, 5L, 5L, 5L, 6L, 6L,
6L, 7L, 7L, 7L, 8L, 8L, 8L), .Label = c("A", "B", "C", "D", "E",
"F", "G", "H"), class = "factor"), Condition = structure(c(1L,
1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L,
2L, 1L, 1L, 1L, 2L, 2L, 2L), .Label = c("No", "Yes"), class = "factor"),
Time = structure(c(1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L,
2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L), .Label = c("1PM",
"2PM", "3PM"), class = "factor")), .Names = c("Score", "Subject",
"Condition", "Time"), class = "data.frame", row.names = c(NA,
## your first try -24L))
m1 <- lmer(Score ~ Condition + Time + Condition*Time + (1 | Subject),
data = myDat)
summary(m1)
#
# Linear mixed model fit by REML ['lmerMod']
# Formula: Score ~ Condition + Time + Condition * Time + (1 | Subject)
# Data: myDat
#
# REML criterion at convergence: 25.7
#
# Scaled residuals:
# Min 1Q Median 3Q Max
# -1.4009 -0.5397 0.1392 0.4146 1.2662
# Random effects:
# Groups Name Variance Std.Dev.
# Subject (Intercept) 0.05619 0.2370
# Residual 0.11329 0.3366
# Number of obs: 24, groups: Subject, 8
#
# Fixed effects:
# Estimate Std. Error t value
# (Intercept) 1.7450 0.2058 8.477
# ConditionYes 1.5600 0.2911 5.359
# Time2PM 0.4925 0.2380 2.069
# Time3PM 1.1500 0.2380 4.832
# ConditionYes:Time2PM -0.2250 0.3366 -0.668
# ConditionYes:Time3PM -0.3950 0.3366 -1.174
#
# Correlation of Fixed Effects:
# (Intr) CndtnY Tim2PM Tim3PM CY:T2P
# ConditionYs -0.707
# Time2PM -0.578 0.409
# Time3PM -0.578 0.409 0.500
# CndtnY:T2PM 0.409 -0.578 -0.707 -0.354
# CndtnY:T3PM 0.409 -0.578 -0.354 -0.707 0.500
## brms - get priors
get_prior(Score ~ Condition + Time + Condition*Time + (1 | Subject),
data = myDat, family= gaussian())
prior_positive <- c(
prior(normal(0, 1), class = Intercept),
prior(normal(0, 1), class = b, lb = 0), # specify lower bound with "lb"
#prior(cauchy(0,1), class = shape),
prior(exponential(1), class = sd),
prior(exponential(1), class = sigma))
mod1_positive <- brm(Score ~ Condition + Time + Condition*Time + (1 | Subject),
data = myDat, family= gaussian(), prior = prior_positive)
summary(mod1_positive)
# Family: gaussian
# Links: mu = identity; sigma = identity
# Formula: Score ~ Condition + Time + Condition * Time + (1 | Subject)
# Data: myDat (Number of observations: 24)
# Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
# total post-warmup samples = 4000
#
# Group-Level Effects:
# ~Subject (Number of levels: 8)
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept) 0.26 0.16 0.02 0.65 1.00 954 1653
#
# Population-Level Effects:
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# Intercept 1.91 0.21 1.48 2.32 1.00 1645 2280
# ConditionYes 1.11 0.27 0.55 1.65 1.00 1431 1742
# Time2PM 0.27 0.17 0.02 0.65 1.00 1861 1176
# Time3PM 0.82 0.22 0.37 1.23 1.00 1793 1129
# ConditionYes:Time2PM 0.28 0.21 0.01 0.76 1.00 2183 1490
# ConditionYes:Time3PM 0.24 0.20 0.01 0.75 1.00 1624 1475
#
# Family Specific Parameters:
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sigma 0.40 0.08 0.27 0.58 1.00 1379 2100
#
# Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
# and Tail_ESS are effective sample size measures, and Rhat is the potential
# scale reduction factor on split chains (at convergence, Rhat = 1).