从具有 factor/character 级输入的拟合模型模拟 draws/data
Simulating draws/data from a fitted model with factor/character level inputs
我有这样的数据：
library(tidyverse)
set.seed(2017)
df <- tibble(
product = c(rep("A", 50), rep("B", 50)),
sales = round(c(rnorm(50, mean = 55, sd = 10), rnorm(50, mean = 60, sd = 15)))
)
我可以对数据建立线性回归：
mod1 <- lm(sales ~ product, data = df)
并预测产品 "A" 和 "B" 的销售额：
predict(mod1, tibble(product = c("A", "B")))
> 1 2
> 55.78 58.96
但我想从拟合模型中模拟绘图，而不是仅仅预测拟合值。我想要抽奖，这样我就可以捕捉到 point-estimate 周围的不确定性（无需使用 SD、CI 等）。
我通常会使用 simulate() 并更改 model_object$fitted.values。但我不能这样做，因为我的模型的输入是 factor/character 级别（"A" 和 "B"）。
我可以得到分布的形状：
a_mu <- coef(summary(mod1))["(Intercept)", "Estimate"]
a_se <- coef(summary(mod1))["(Intercept)", "Std. Error"]
b_mu <- coef(summary(mod1))["productB", "Estimate"]
b_se <- coef(summary(mod1))["productB", "Std. Error"]
然后像这样模拟开奖：
N <- 100
product_A <- replicate(N,
rnorm(n = 1, mean = a_mu, sd = a_se) + rnorm(n = 1, mean = b_mu, sd = b_se) * 0)
product_B <- replicate(N,
rnorm(n = 1, mean = a_mu, sd = a_se) + rnorm(n = 1, mean = b_mu, sd = b_se) * 1)
并将其全部塞入小标题以供可视化：
pred <- tibble(A = product_A, B = product_B)
但是这个过程看起来超级简陋。如果我的数据增长到例如 5 个输入变量，每个变量有 10 个因子水平，则不会扩展。那么，我怎样才能使它具有普遍性？
我更愿意留在基地 R and/or tidyverse。是的，我知道我在这里与贝叶斯统计调情，我也许可以使用 Stan 从后部绘制...但这不是重点。
对于点估计的不确定性，1）如果选择模拟，我推荐boxplot。 2）如果选择CI，可以手动计算，或者使用predict()比如Webb的评论，plot intervals。在这里，我只是向您展示如何以通用形式进行模拟。您快到了，所以希望这对您有所帮助。
myfactor_pred<-function(factor,N){
if(factor==0){
return(rnorm(N,coef(summary(mod1))[1,1],coef(summary(mod1))[1,2]))
}else{
return(rnorm(N,coef(summary(mod1))[1,1],coef(summary(mod1))[1,2])+
rnorm(N,coef(summary(mod1))[2,1],coef(summary(mod1))[2,2]))
}
}
A<-myfactor_pred(0,100)#call function and get simulation for A
B<-myfactor_pred(1,100)#call function and get simulation for B
boxplot(data.frame(A,B),xlab="product",ylab="sales")
我相信，如果您想显示预测的不确定性，贝叶斯回归比传统回归更适合。
也就是说，你可以通过以下方式得到你想要的（你将不得不重命名SimulatedMat的列）：
# All the possible combinations of factors
modmat<-unique(model.matrix(sales ~.,df))
# Number of simulations
simulations<-100L
# initialise result matrix
SimulatedMat<-matrix(0,nrow=simulations,ncol=0)
# iterate amongst all combinations of factors
for(i in 1:nrow(modmat)){
# columns with value one
selcols<-which(modmat[i,]==1)
# simulation for the factors with value 1
simul<-apply(mapply(rnorm,n=simulations,coef(summary(mod1))[selcols, "Estimate"],
coef(summary(mod1))[selcols, "Std. Error"]),1,sum)
# incorporate result to the matrix
SimulatedMat<-cbind(SimulatedMat,simul)
}
嗯，我们来看看Bayesian linear regression, do some sampling, and then compare to frequentist prediction intervals. We'll try to follow the notation in the linked Wikipedia page, where we'll be working on the posterior distribution。
设置贝叶斯后验模拟
X <- model.matrix(sales~product,data=df)
n <- nrow(X)
k <- ncol(X)
v <- n - k
y <- df$sales
# Take Lambda_0 <- 0 so these simplify
beta.hat <- solve(crossprod(X),crossprod(X,y))
S <- solve(crossprod(X)) # Lambda_n^-1
mu_n <- beta.hat
a_n <- v/2 # I think this is supposed to be v instead of n, the factor with k was dropped?
b_n <- (crossprod(y) - crossprod(mu_n, crossprod(X) %*% mu_n))/2
运行模拟
现在我们从反伽马中得出 sigma^2 a_n, b_n
N <- 10000
sigma.2 <- 1/rgamma(N, a_n, b_n)
并从正常的 u_n、S * sigma.2（刚刚生成）
中绘制 beta
require("MASS")
beta <- sapply(sigma.2, function(s) MASS::mvrnorm(1, mu_n, S * s))
我们将把这些全部放入 data.frame
sim <- data.frame(t(beta),sigma=sqrt(sigma.2))
与 lm 输出相比
统计数据
t(sapply(sim,function(x) c(mean=mean(x),sd=sd(x))))
mean sd
X.Intercept. 55.786709 1.9585332
productB 3.159653 2.7552841
sigma 13.736069 0.9932521
与 lm 比较
mod1 <- lm(sales ~ product, data = df)
summary(mod1)
...
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 55.780 1.929 28.922 <2e-16 ***
productB 3.180 2.728 1.166 0.246
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 13.64 on 98 degrees of freedom
...
虽然您会注意到贝叶斯方法往往更保守（更大的标准误差）。
模拟点数
A 和 B 因子水平的模拟预测为
d <- unique(df$product)
mm <- cbind(1,contrasts(d))
sim.y <- crossprod(beta,t(mm))
head(sim.y)
A B
[1,] 56.09510 61.45903
[2,] 55.43892 57.87281
[3,] 58.49551 60.59784
[4,] 52.55529 62.71117
[5,] 62.18198 59.27573
[6,] 59.50713 57.39560
与 lm 的置信区间比较。
我们可以根据 A 和 B
的模拟值计算分位数
t(apply(sim.y,2,function(col) quantile(col,c(0.025,0.975))))
2.5% 97.5%
A 51.90695 59.62353
B 55.14255 62.78334
并与频率线性回归的置信区间进行比较
predict(mod1, data.frame(product = c("A", "B")), interval="confidence",level=0.95)
fit lwr upr
1 55.78 51.95266 59.60734
2 58.96 55.13266 62.78734
数据
# those without tidyverse, this will suffice
if (!require("tidyverse")) tibble <- data.frame
set.seed(2017)
df <- tibble(
product = c(rep("A", 50), rep("B", 50)),
sales = round(c(rnorm(50, mean = 55, sd = 10), rnorm(50, mean = 60, sd = 15)))
)
Gelman 和 Hill (2007)1 提供了一个 Bayesian-flavored 函数，用于使用模拟估计频率回归中的不确定性。该功能在他们的（恕我直言）文本中从第 142 页底部开始描述，可以是 viewed on google books。
该函数称为 sim，可从 arm 包（Gelman 和 Hill 的文本附带的包）中获得。它使用模型参数（包括考虑系数的协方差和标准误差）来模拟系数的联合分布。该函数自本书出版后发生了变化，现在 returns 一个使用槽访问的 S4 对象，因此实际实现与书中描述的略有不同。
这是一个使用您的数据的示例：
library(ggplot2)
library(ggbeeswarm)
theme_set(theme_classic())
library(arm)
首先，我们将使用 sim 函数对模型系数进行 1000 次模拟：
sim.mod = sim(mod1, 1000)
每个模拟的系数可以在sim.mod@coef中找到，这是一个矩阵。这是前四行：
sim.mod@coef[1:4,]
(Intercept) productB
[1,] 55.25320 3.5782625
[2,] 59.90534 0.4608387
[3,] 55.79126 5.1872595
[4,] 57.97446 1.0012827
现在让我们提取系数模拟，将它们转换为数据框，并缩短列名。这将为我们提供一个数据框 sc，其中一列用于模拟截距，一列用于虚拟变量的模拟系数 product=="B":
sc = setNames(as.data.frame(sim.mod@coef), c("Int","prodB"))
从这里，您可以使用模拟来评估系数和预测销售额的不确定性和可能范围。下面是一些可视化。
让我们为每对模拟系数绘制一条蓝色回归线。我们将得到 1,000 条线，线的密度将向我们展示最可能的系数组合。我们还以黄色显示拟合回归线，以红色显示基础数据点。显然，这些线仅在 x-axis 上的 A 和 B 点有意义。这类似于 Gelman 和 Hill 在他们的书中展示模拟结果的方式。
ggplot() +
geom_abline(data=sc, aes(slope=prodB, intercept=Int), colour="blue", alpha=0.03) +
geom_beeswarm(data=df, aes(product, sales), alpha=1, colour="red", size=0.7) +
geom_abline(slope=coef(mod1)[2], intercept=coef(mod1)[1], colour="yellow", size=0.8)
另一种选择是针对每对模拟系数计算每种产品的预测平均销售额。我们在下面这样做并将结果绘制为小提琴图。此外，我们还包括平均销售额的中值预测，以及平均销售额的 2.5% 到 97.5% 分位数的范围：
pd = data.frame(product=rep(c("A","B"), each=1000), sc)
pd$sales = ifelse(pd$product=="A", pd$Int, pd$Int + pd$prodB)
ggplot(pd, aes(product, sales)) +
geom_violin() +
stat_summary(fun.data=median_hilow, colour="red", geom="errorbar", width=0.05, size=0.8, alpha=0.6) +
stat_summary(fun.y=mean, aes(label=round(..y..,1)), geom="text", size=4, colour="blue")
最后，我们用 50% 和 95% 的椭圆绘制模拟系数值的分布。 coord_equal() 确保一个单元在水平和垂直轴上覆盖相同的物理距离。截距（横轴）是product=="A"时销售额的预测值。当 product=="B":
时，斜率（纵轴）是预测的销售额差异（相对于 product=="A"）
ggplot(sc, aes(Int, prodB)) +
geom_point(alpha=0.5, colour="red", size=1) +
stat_ellipse(level=c(0.5), colour="blue") +
stat_ellipse(level=c(0.95), colour="blue") +
coord_equal() +
scale_x_continuous(breaks=seq(50,62,2)) +
scale_y_continuous(breaks=seq(-6,12,2))
如果你有多个变量，可视化会更复杂，但原理类似于上面说明的 single-predictor 案例。 sim 函数将使用多个预测变量和具有多个级别的分类变量，因此这种方法应该扩展到更复杂的数据集。
- 一个。 Gelman 和 J. Hill，使用回归和 Multilevel/Hierarchical 模型进行数据分析，剑桥大学出版社（2007 年）。
我有这样的数据：
library(tidyverse)
set.seed(2017)
df <- tibble(
product = c(rep("A", 50), rep("B", 50)),
sales = round(c(rnorm(50, mean = 55, sd = 10), rnorm(50, mean = 60, sd = 15)))
)
我可以对数据建立线性回归：
mod1 <- lm(sales ~ product, data = df)
并预测产品 "A" 和 "B" 的销售额：
predict(mod1, tibble(product = c("A", "B")))
> 1 2
> 55.78 58.96
但我想从拟合模型中模拟绘图，而不是仅仅预测拟合值。我想要抽奖，这样我就可以捕捉到 point-estimate 周围的不确定性（无需使用 SD、CI 等）。
我通常会使用 simulate() 并更改 model_object$fitted.values。但我不能这样做，因为我的模型的输入是 factor/character 级别（"A" 和 "B"）。
我可以得到分布的形状：
a_mu <- coef(summary(mod1))["(Intercept)", "Estimate"]
a_se <- coef(summary(mod1))["(Intercept)", "Std. Error"]
b_mu <- coef(summary(mod1))["productB", "Estimate"]
b_se <- coef(summary(mod1))["productB", "Std. Error"]
然后像这样模拟开奖：
N <- 100
product_A <- replicate(N,
rnorm(n = 1, mean = a_mu, sd = a_se) + rnorm(n = 1, mean = b_mu, sd = b_se) * 0)
product_B <- replicate(N,
rnorm(n = 1, mean = a_mu, sd = a_se) + rnorm(n = 1, mean = b_mu, sd = b_se) * 1)
并将其全部塞入小标题以供可视化：
pred <- tibble(A = product_A, B = product_B)
但是这个过程看起来超级简陋。如果我的数据增长到例如 5 个输入变量，每个变量有 10 个因子水平，则不会扩展。那么，我怎样才能使它具有普遍性？
我更愿意留在基地 R and/or tidyverse。是的，我知道我在这里与贝叶斯统计调情，我也许可以使用 Stan 从后部绘制...但这不是重点。
对于点估计的不确定性，1）如果选择模拟，我推荐boxplot。 2）如果选择CI，可以手动计算，或者使用predict()比如Webb的评论，plot intervals。在这里，我只是向您展示如何以通用形式进行模拟。您快到了，所以希望这对您有所帮助。
myfactor_pred<-function(factor,N){
if(factor==0){
return(rnorm(N,coef(summary(mod1))[1,1],coef(summary(mod1))[1,2]))
}else{
return(rnorm(N,coef(summary(mod1))[1,1],coef(summary(mod1))[1,2])+
rnorm(N,coef(summary(mod1))[2,1],coef(summary(mod1))[2,2]))
}
}
A<-myfactor_pred(0,100)#call function and get simulation for A
B<-myfactor_pred(1,100)#call function and get simulation for B
boxplot(data.frame(A,B),xlab="product",ylab="sales")
我相信，如果您想显示预测的不确定性，贝叶斯回归比传统回归更适合。
也就是说，你可以通过以下方式得到你想要的（你将不得不重命名SimulatedMat的列）：
# All the possible combinations of factors
modmat<-unique(model.matrix(sales ~.,df))
# Number of simulations
simulations<-100L
# initialise result matrix
SimulatedMat<-matrix(0,nrow=simulations,ncol=0)
# iterate amongst all combinations of factors
for(i in 1:nrow(modmat)){
# columns with value one
selcols<-which(modmat[i,]==1)
# simulation for the factors with value 1
simul<-apply(mapply(rnorm,n=simulations,coef(summary(mod1))[selcols, "Estimate"],
coef(summary(mod1))[selcols, "Std. Error"]),1,sum)
# incorporate result to the matrix
SimulatedMat<-cbind(SimulatedMat,simul)
}
嗯，我们来看看Bayesian linear regression, do some sampling, and then compare to frequentist prediction intervals. We'll try to follow the notation in the linked Wikipedia page, where we'll be working on the posterior distribution。
设置贝叶斯后验模拟
X <- model.matrix(sales~product,data=df)
n <- nrow(X)
k <- ncol(X)
v <- n - k
y <- df$sales
# Take Lambda_0 <- 0 so these simplify
beta.hat <- solve(crossprod(X),crossprod(X,y))
S <- solve(crossprod(X)) # Lambda_n^-1
mu_n <- beta.hat
a_n <- v/2 # I think this is supposed to be v instead of n, the factor with k was dropped?
b_n <- (crossprod(y) - crossprod(mu_n, crossprod(X) %*% mu_n))/2
运行模拟
现在我们从反伽马中得出 sigma^2 a_n, b_n
N <- 10000
sigma.2 <- 1/rgamma(N, a_n, b_n)
并从正常的 u_n、S * sigma.2（刚刚生成）
中绘制 betarequire("MASS")
beta <- sapply(sigma.2, function(s) MASS::mvrnorm(1, mu_n, S * s))
我们将把这些全部放入 data.frame
sim <- data.frame(t(beta),sigma=sqrt(sigma.2))
与 lm 输出相比
统计数据
t(sapply(sim,function(x) c(mean=mean(x),sd=sd(x))))
mean sd
X.Intercept. 55.786709 1.9585332
productB 3.159653 2.7552841
sigma 13.736069 0.9932521
与 lm 比较
mod1 <- lm(sales ~ product, data = df)
summary(mod1)
...
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 55.780 1.929 28.922 <2e-16 ***
productB 3.180 2.728 1.166 0.246
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 13.64 on 98 degrees of freedom
...
虽然您会注意到贝叶斯方法往往更保守（更大的标准误差）。
模拟点数
A 和 B 因子水平的模拟预测为
d <- unique(df$product)
mm <- cbind(1,contrasts(d))
sim.y <- crossprod(beta,t(mm))
head(sim.y)
A B
[1,] 56.09510 61.45903
[2,] 55.43892 57.87281
[3,] 58.49551 60.59784
[4,] 52.55529 62.71117
[5,] 62.18198 59.27573
[6,] 59.50713 57.39560
与 lm 的置信区间比较。
我们可以根据 A 和 B
t(apply(sim.y,2,function(col) quantile(col,c(0.025,0.975))))
2.5% 97.5%
A 51.90695 59.62353
B 55.14255 62.78334
并与频率线性回归的置信区间进行比较
predict(mod1, data.frame(product = c("A", "B")), interval="confidence",level=0.95)
fit lwr upr
1 55.78 51.95266 59.60734
2 58.96 55.13266 62.78734
数据
# those without tidyverse, this will suffice
if (!require("tidyverse")) tibble <- data.frame
set.seed(2017)
df <- tibble(
product = c(rep("A", 50), rep("B", 50)),
sales = round(c(rnorm(50, mean = 55, sd = 10), rnorm(50, mean = 60, sd = 15)))
)
Gelman 和 Hill (2007)1 提供了一个 Bayesian-flavored 函数，用于使用模拟估计频率回归中的不确定性。该功能在他们的（恕我直言）文本中从第 142 页底部开始描述，可以是 viewed on google books。
该函数称为 sim，可从 arm 包（Gelman 和 Hill 的文本附带的包）中获得。它使用模型参数（包括考虑系数的协方差和标准误差）来模拟系数的联合分布。该函数自本书出版后发生了变化，现在 returns 一个使用槽访问的 S4 对象，因此实际实现与书中描述的略有不同。
这是一个使用您的数据的示例：
library(ggplot2)
library(ggbeeswarm)
theme_set(theme_classic())
library(arm)
首先，我们将使用 sim 函数对模型系数进行 1000 次模拟：
sim.mod = sim(mod1, 1000)
每个模拟的系数可以在sim.mod@coef中找到，这是一个矩阵。这是前四行：
sim.mod@coef[1:4,]
(Intercept) productB [1,] 55.25320 3.5782625 [2,] 59.90534 0.4608387 [3,] 55.79126 5.1872595 [4,] 57.97446 1.0012827
现在让我们提取系数模拟，将它们转换为数据框，并缩短列名。这将为我们提供一个数据框 sc，其中一列用于模拟截距，一列用于虚拟变量的模拟系数 product=="B":
sc = setNames(as.data.frame(sim.mod@coef), c("Int","prodB"))
从这里，您可以使用模拟来评估系数和预测销售额的不确定性和可能范围。下面是一些可视化。
让我们为每对模拟系数绘制一条蓝色回归线。我们将得到 1,000 条线，线的密度将向我们展示最可能的系数组合。我们还以黄色显示拟合回归线，以红色显示基础数据点。显然，这些线仅在 x-axis 上的 A 和 B 点有意义。这类似于 Gelman 和 Hill 在他们的书中展示模拟结果的方式。
ggplot() +
geom_abline(data=sc, aes(slope=prodB, intercept=Int), colour="blue", alpha=0.03) +
geom_beeswarm(data=df, aes(product, sales), alpha=1, colour="red", size=0.7) +
geom_abline(slope=coef(mod1)[2], intercept=coef(mod1)[1], colour="yellow", size=0.8)
另一种选择是针对每对模拟系数计算每种产品的预测平均销售额。我们在下面这样做并将结果绘制为小提琴图。此外，我们还包括平均销售额的中值预测，以及平均销售额的 2.5% 到 97.5% 分位数的范围：
pd = data.frame(product=rep(c("A","B"), each=1000), sc)
pd$sales = ifelse(pd$product=="A", pd$Int, pd$Int + pd$prodB)
ggplot(pd, aes(product, sales)) +
geom_violin() +
stat_summary(fun.data=median_hilow, colour="red", geom="errorbar", width=0.05, size=0.8, alpha=0.6) +
stat_summary(fun.y=mean, aes(label=round(..y..,1)), geom="text", size=4, colour="blue")
最后，我们用 50% 和 95% 的椭圆绘制模拟系数值的分布。 coord_equal() 确保一个单元在水平和垂直轴上覆盖相同的物理距离。截距（横轴）是product=="A"时销售额的预测值。当 product=="B":
product=="A"）
ggplot(sc, aes(Int, prodB)) +
geom_point(alpha=0.5, colour="red", size=1) +
stat_ellipse(level=c(0.5), colour="blue") +
stat_ellipse(level=c(0.95), colour="blue") +
coord_equal() +
scale_x_continuous(breaks=seq(50,62,2)) +
scale_y_continuous(breaks=seq(-6,12,2))
如果你有多个变量，可视化会更复杂，但原理类似于上面说明的 single-predictor 案例。 sim 函数将使用多个预测变量和具有多个级别的分类变量，因此这种方法应该扩展到更复杂的数据集。
- 一个。 Gelman 和 J. Hill，使用回归和 Multilevel/Hierarchical 模型进行数据分析，剑桥大学出版社（2007 年）。