3

我的数据如下所示：

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

但我想模拟拟合模型的绘制，而不是仅仅预测拟合值。我想要绘制，以便我可以捕捉点估计周围的不确定性（无需使用 SD、CI 等）。

我通常会使用simulate()和更改model_object$fitted.values. 但我不能这样做，因为我的模型的输入是因子/字符级别（“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 和/或tidyverse. 是的，我知道我在这里与贝叶斯统计调情，我也许可以使用 Stan 从后验中提取……但这不是重点。

4

4 回答 4

4

Gelman 和 Hill (2007) 1提供了一个贝叶斯风格的函数，用于使用模拟估计频率回归中的不确定性。该功能在其（恕我直言）文本中从第 142 页底部开始描述，可以在 google 图书上查看。

该函数被调用sim并且可以从arm包中获得（这是伴随着 Gelman 和 Hill 的文本的包）。它使用模型参数（包括考虑系数的协方差和标准误差）来模拟系数的联合分布。自从本书出版以来，该函数发生了变化，现在返回一个使用插槽访问的 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 轴上的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=="A"），当product=="B"：

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))

在此处输入图像描述

如果您有多个变量，可视化会更加复杂，但原理类似于上面说明的单预测器情况。该sim函数将使用多个预测变量和具有多个级别的分类变量，因此这种方法应该扩展到更复杂的数据集。

  1. A. Gelman 和 J. Hill，使用回归和多级/分层模型进行数据分析，剑桥大学出版社（2007 年）。
于 2017-02-10T04:03:03.417 回答
2

好吧，让我们看一下贝叶斯线性回归，做一些抽样，然后与常客预测区间进行比较。我们将尝试遵循链接的 Wikipedia 页面中的符号，我们将在其中研究后验分布。

设置贝叶斯后验模拟

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

运行模拟

现在我们从反伽玛 a_n, b_n 中绘制 sigma^2

N <- 10000
sigma.2 <- 1/rgamma(N, a_n, b_n)

并从正常的 u_n 中绘制 beta，S * sigma.2（刚刚生成）

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)))
)
于 2017-02-09T21:38:33.753 回答
0

对于点估计的不确定性，1）如果您选择模拟，我会推荐箱线图。2）如果选择CI，可以手动计算，或者使用Webb评论中的predict()，绘制区间。在这里，我只是向您展示如何以通用形式进行模拟。你快到了，所以希望这会有所帮助。

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")

在此处输入图像描述

于 2017-02-07T00:31:21.707 回答
0

我相信，如果您想显示预测的不确定性，贝叶斯回归比传统回归更适合。

也就是说，您可以通过以下方式获得所需的内容（您必须重命名 的列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)
}
于 2017-02-07T10:47:28.323 回答