PROSAGA码农传奇-模型可视化-使用因子/字符级输入模拟拟合模型的绘图/数据

0# 易米烊光 | 2019-08-31 10-32

<div class =“post-text”itemprop =“text”>
  <P>
    好吧，我们来看看吧
    <a href="https://en.wikipedia.org/wiki/Bayesian_linear_regression" rel="nofollow noreferrer">
      贝叶斯线性回归
    </A>
    ，做一些采样，然后与频率预测间隔进行比较。我们将尝试按照链接的维基百科页面中的符号进行操作，我们将在其中处理
    <a href="https://en.wikipedia.org/wiki/Bayesian_linear_regression#Posterior_distribution" rel="nofollow noreferrer">
      后验分布
    </A>
    。
  </p>
  <H2>
    建立贝叶斯后验模拟
  </H2>
   <pre>
    <code>
      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

</code>
  </pre>
  <H2>
    运行模拟
  </H2>
  <P>
    现在我们从反伽马a_n，b_n中绘制sigma ^ 2
  </p>
   <pre>
    <code>
      N <- 10000
sigma.2 <- 1/rgamma(N, a_n, b_n)

</code>
  </pre>
  <P>
    并从正常的u_n，S * sigma.2（刚生成）中绘制beta
  </p>
   <pre>
    <code>
      require("MASS")
beta <- sapply(sigma.2, function(s) MASS::mvrnorm(1, mu_n, S * s))

</code>
  </pre>
  <P>
    我们将把这一切都放到data.frame中
  </p>
   <pre>
    <code>
      sim <- data.frame(t(beta),sigma=sqrt(sigma.2))

</code>
  </pre>
  <H2>
    比较
     <code>
      lm
    </code>
     产量
  </H2>
  <P>
    统计数据
  </p>
   <pre>
    <code>
      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

</code>
  </pre>
  <P>
    比较密切
     <code>
      lm
    </code>
  </p>
   <pre>
    <code>
      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
...

</code>
  </pre>
  <P>
    虽然你会注意到贝叶斯方法往往更保守（更大的标准误差）。
  </p>
  <H2>
    模拟点
  </H2>
  <P>
    模拟预测
     <code>
      A
    </code>
     和
     <code>
      B
    </code>
     因子水平是
  </p>
   <pre>
    <code>
      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

</code>
  </pre>
  <H2>
    比较来自的置信区间
     <code>
      lm
    </code>
    。
  </H2>
  <P>
    我们可以计算模拟值的分位数
     <code>
      A
    </code>
     和
     <code>
      B
    </code>
  </p>
   <pre>
    <code>
      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

</code>
  </pre>
  <P>
    并与频率线性回归的置信区间进行比较
  </p>
   <pre>
    <code>
      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

</code>
  </pre>
  <H2>
    数据
  </H2>
   <pre>
    <code>
      # 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)))
)

</code>
  </pre>
</DIV>

1# 淋了一整夜的雨 | 2019-08-31 10-32

<div class =“post-text”itemprop =“text”>
  <P>
    格尔曼和希尔（2007）
    <SUP>
      1
    </SUP>
     提供贝叶斯风格的函数，用于使用模拟估计频率回归中的不确定性。该功能在其（IMHO优秀）文本的第142页底部开始描述，可以是
    <a href =“https://books.google.com/books?id=lV3DIdV0F9AC&pg=PA142&lpg=PA142&dq=%22to+get+n.sims+simulation+draws%22&source=bl& ots = 6klLHfSyU0＆amp; sig = XqPRDzyNQS4LCF7yqtx5Gcd4clk＆amp; hl = en＆amp; sa = X＆amp; ved = 0ahUKEwjq0tnDxITSAhUDz2MKHWP0ChQQ6AEIHDAA #v = onepage＆amp; q =％22to％20get％20n.sims％20simulation％20draws％22＆amp; f = false“rel =”noreferrer“ >
      在谷歌书籍上查看
    </A>
    。
  </p>
  <P>
    该函数被调用
     <code>
      sim
    </code>
     并可从
     <code>
      arm
    </code>
     包（这是Gelman和Hill的文本附带的包）。它使用模型参数（包括考虑系数的协方差和标准误差）来模拟系数的联合分布。自本书出版以来，该函数已发生变化，现在返回一个使用插槽访问的S4对象，因此实际实现与本书中描述的略有不同。
  </p>
  <P>
    以下是使用您的数据的示例：
  </p>
   <pre>
    <code>
      library(ggplot2)
library(ggbeeswarm)
theme_set(theme_classic())
library(arm)

</code>
  </pre>
  <P>
    首先，我们将使用。生成1000个模型系数的模拟
     <code>
      sim
    </code>
     功能：
  </p>
   <pre>
    <code>
      sim.mod = sim(mod1, 1000)

</code>
  </pre>
  <P>
    每个模拟的系数可以在中找到
     <code>
      sim.mod@coef
    </code>
    ，这是一个矩阵。以下是前四行：
  </p>
   <pre>
    <code>
      sim.mod@coef[1:4,]

</code>
  </pre>
  <BLOCKQUOTE>
     <pre>
      <code>
             (Intercept)  productB
[1,]    55.25320 3.5782625
[2,]    59.90534 0.4608387
[3,]    55.79126 5.1872595
[4,]    57.97446 1.0012827

</code>
    </pre>
  </BLOCKQUOTE>
  <P>
    现在让我们提取系数模拟，将它们转换为数据框，并缩短列名。这将为我们提供数据框架
     <code>
      sc
    </code>
     一列用于模拟截距，一列用于虚拟变量的模拟系数
     <code>
      product=="B"
    </code>
    ：
  </p>
   <pre>
    <code>
      sc = setNames(as.data.frame(sim.mod@coef), c("Int","prodB"))

</code>
  </pre>
  <P>
    从这里，您可以使用模拟来评估系数和预测销售的不确定性和可能的范围。以下是一些可视化。
  </p>
  <P>
    让我们为每个模拟系数对绘制蓝色回归线。我们将获得1,000行，并且线的密度将向我们显示最可能的系数组合。我们还显示黄色的拟合回归线和红色的基础数据点。显然这些线条只对它有意义
     <code>
      A
    </code>
     和
     <code>
      B
    </code>
     x轴上的点。这类似于Gelman和Hill在他们的书中展示模拟结果的方式。
  </p>
   <pre>
    <code>
      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)

</code>
  </pre>
  <P>
    <a href="https://i.stack.imgur.com/Xp4Rw.png" rel="noreferrer">
      <img src =“https://i.stack.imgur.com/Xp4Rw.png”alt =“在此输入图像说明”/>
    </A>
  </p>
  <P>
    另一种选择是计算每对模拟系数的每种产品的预测平均销售额。我们在下面这样做，并将结果绘制为小提琴情节。此外，我们还包括平均销售额的中位数预测，以及2.5％到97.5％的平均销售额分布范围：
  </p>
   <pre>
    <code>
      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")

</code>
  </pre>
  <P>
    <a href="https://i.stack.imgur.com/YJBda.png" rel="noreferrer">
      <img src =“https://i.stack.imgur.com/YJBda.png”alt =“在此处输入图片说明”/>
    </A>
  </p>
  <P>
    最后，我们用50％和95％的椭圆绘制模拟系数值的分布。
     <code>
      coord_equal()
    </code>
     确保一个单元在水平和垂直轴上覆盖相同的物理距离。截距（横轴）是销售时的预测值
     <code>
      product=="A"
    </code>
    。斜率（纵轴）是预测的销售差异（相对于
     <code>
      product=="A"
    </code>
    ） 什么时候
     <code>
      product=="B"
    </code>
    ：
  </p>
   <pre>
    <code>
      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))

</code>
  </pre>
  <P>
    <a href="https://i.stack.imgur.com/RK98D.png" rel="noreferrer">
      <img src =“https://i.stack.imgur.com/RK98D.png”alt =“在此处输入图片说明”/>
    </A>
  </p>
  <P>
    如果您有多个变量，可视化将会更复杂，但原则类似于上面说明的单预测器情况。该
     <code>
      sim
    </code>
     函数将与多个预测变量和具有多个级别的分类变量一起使用，因此这种方法应该扩展到更复杂的数据集。
  </p>
  <OL>
    <LI>
      A. Gelman和J. Hill，
      <EM>
        使用回归和多级/分层模型的数据分析
      </EM>
      ，剑桥大学出版社（2007年）。
    </LI>
  </醇>
</DIV>

2# Let us fly | 2019-08-31 10-32

3# v-star*위위 | 2019-08-31 10-32