乐趣区

关于数据挖掘:R语言两层2^k析因试验设计因子设计分析工厂产量数据和Lenth方法检验显著性可视化

jiezi

原文链接：http://tecdat.cn/?p=25921 

假如考察人员有趣味查看减肥干涉办法的三个组成部分。这三个组成部分是：

  1. 记录食物日记（是 / 否）
  2. 减少流动（是 / 否）
  3. 家访（是 / 否）

调查员打算考察所有 ，试验条件的组合。试验条件为

  • 要执行因子设计，您须要为多个因子（变量）中的每一个抉择固定数量的程度，而后以所有可能的组合运行试验。
  • 这些因素能够是定量的或定性的。
  • 定量变量的两个程度能够是两个不同的温度或两个不同的浓度。
  • 定性因素可能是两种类型的催化剂或某些实体的存在和不存在。

符号 ：– 因子数 (3) – 每个因子的程度数 (2) – 设计中有多少试验条件 ()

因子试验能够波及具备不同程度数量的因子。

测试：

思考一个  设计。

  1. 有多少因素？
  2. 每个因素有多少个程度？
  3. 多少试验条件（运行）？

答案：

(a) 有 2+2+1 = 5 个因数。

(b) 两个因素有 4 个程度，2 个因素有 3 个程度，1 个因素有 2 个程度。

(c) 有 288 个试验条件。

方差分析和因子设计之间的区别

在 ANOVA 中，指标是比拟各个试验条件。

让咱们考虑一下下面的食物日记钻研。

咱们能够通过比拟食物日记设置为 NO（条件 1-4）的所有条件的平均值和食物日记设置为 YES（条件 5-8）的所有条件的平均值来预计食物日记的成果。这也被称为食物日记的  主效应 ，形容词 _次要_ 是揭示这个平均值超过了其余因素的程度。

食物日记的次要作用是：

体育锻炼的次要作用是：

家访的次要作用是：

应用了所有试验对象，但重新排列以进行每次比拟。受试者被回收以测量不同的成果。这是析因实验更无效的起因之一。

执行  因子设计

要执行因子设计：

  1. 为每个因子抉择固定数量的程度。
  2. 以所有可能的组合运行试验。

咱们将探讨每个因子只有两个程度的设计。因素能够是定量的或定性的。两个程度的定量变量能够是两个不同的温度或浓度。定量变量的两个级别能够是两种不同类型的催化剂或某些实体的存在 / 不存在。

一项试验采纳 2^3 因子设计，具备两个定量因素 – 温度 (T) 和浓度 (C) – 以及一个定性因素 – 催化剂 K 类型。

温度 T（C∘）有两个等级：160C∘和 180C∘。它们别离编码为 -1 和 +1。

浓度 C (%) 有两个级别：20 和 40。它们别离编码为 -1 和 +1。

催化剂 K 有两个级别：A 和 B。它们别离编码为 -1 和 +1。

记录的每个数据值都是针对两次反复运行的均匀因变量产量 y。

立方图

下图显示了立方体角处因子 T、C 和 K 的各种组合的 y 值。例如，当 T=-1、C = 1 和 K=-1 时，从运行 3 取得 y=54。

  • 立方体展现了这种设计如何沿着立方体的 12 个边缘进行 12 次比拟：温度变动影响的四个测量值；浓度变动影响的四种测量方法；催化剂变动成果的四种测量方法。
  • 在立方体的每条边上，只有一个因子发生变化，而其余两个因子放弃不变。
bh4 <- lm
Plot

因子效应

次要影响

运行 1 和 2 的影响仅因温度而不同，因为浓度为 20%，催化剂类型为 A。差别 72-60 = 12 提供了一种温度影响的测量值，而其余因素放弃不变。对于浓度和催化剂的四种组合中的每一种，有四种这样的温度效应测量方法。

T 的次要（均匀）影响是

有一组相似的浓度 C 测量值。在这些测量值中的每一个中，程度 T 和 K 都放弃不变。浓度 C 的次要影响是：

C 的次要（均匀）影响是

K 的次要影响是

K 的次要（均匀）影响是

所有 8 次运行都用于预计每个主效应。这就是因子设计比一次查看一个因子更无效的起因。

一般来说，次要影响是两个平均值之间的差别：

其中 ¯y+ 是对应于因子 +1 程度的均匀响应，而 ¯y− 是对应于因子 -1 程度的均匀响应。

交互成果

两因素相互作用

当催化剂 K 为 A 时，温度效应为：

当催化剂 K 为 B 时，温度效应为：

这两个均匀差别之间的均匀差别称为温度和催化剂之间的  相互作用 ，用 TK 示意。这就是温度和催化剂两个因素之间的相互作用——温度和催化剂之间的两个因素相互作用。

这也能够在立方图上看到：与立方体侧面 (13) 相比，立方体 (33) 反面的平均温度影响更大。

三因素相互作用

当催化剂为 B（在其 +1 程度）时，浓度相互作用的温度为：

当催化剂为 A（在其 -1 程度）时，浓度相互作用的温度为：

这两种相互作用之间的差别掂量了两种催化剂的温度 - 浓度相互作用的一致性。这种差别的一半被定义为温度、浓度和催化剂的三因素相互作用，用 TCK 示意。

因子设计中的复制

工厂试验的后果 y 是两次反复运行的平均值。两个独自的运行如下表所示。运行程序是随机的。例如，运行 6 和 13 是 T、C 和 K（T=-1、C=-1、K=-1）的雷同设置下的两个仿行。

复制运行并不总是可行的。工厂试验运行包含清洁反应器，插入适当的催化剂装料，并在给定的进料浓度下在给定的温度下运行设施 3 小时，以使过程在所选的试验条件下稳定下来，以及 (4) 取样在运行的最初几个小时内每 15 分钟输入一次。

假如每次测量的方差为 σ2。每组条件下的预计方差为：

其中 yi1 是第 i 次运行的第一个后果。在上表中 diffi=(yi1−yi2)。σ2 的汇总预计是 

对于反复运行，具备一个自由度的方差预计为 . 这些产生单自由度预计的平均值产生具备 8 个自由度的合并预计 s2=8。

反复运行效应的误差方差和标准误差的预计

每个预计的效应，例如 T、C、K、TC 等，都是 8 个观测值的两个平均值之间的差别。反复运行的因子效应方差为

因而，任何因子效应的规范误为：

后果解释

哪些影响是实在的，哪些能够偶尔解释？一个粗略的教训法令是，任何 2-3 倍于其标准误差的效应都不容易仅靠必然性来解释。

如果咱们假如观测值是独立且正态分布的，那么

因而，95% 的置信区间能够计算为：

其中 t8,.05/2 是 t8t 的第 97.5 个百分位数。这是通过 qt() 函数在 R 中取得的。

qt(p = 1-.025,df = 8)

因而，因子效应的 95% 置信区间为

T 的 95% 置信区间为

K 的 95% 置信区间为 

1.5-3.2 #上限 
## \[1\] -1.7
1.5+3.2 #下限 
## \[1\] 4.7

温度的影响可能不是偶尔的，但偶尔不能成为催化剂的影响的规定。

只有在没有证据表明该因素与其余因素相互作用时，才应独自解释一个因素的主效应。

交互图

下图显示了每对因子 TC、TK、CK（即这些因子的每个因子程度组合）的均匀产量。这些图通常称为交互图。如果两条线平行，则表明没有相互作用，如果两条线穿插或靠近穿插，则表明可能存在相互作用。

下图显示了催化剂和温度之间的双向相互作用。

plot(tabT,taC,ty, type = "l")

2k 因子设计的线性模型

yi 是第 i 次运行的后果，

2^3 因子设计的线性模型是：

变量  是温度和浓度之间的相互作用，xi1xi3xi1xi3 是温度和催化剂之间的相互作用等。

参数估计是通过 lm() R 中的函数取得的。

fact.mod <-lm(y~T\*K\*C,data = tab0503)
round(summary(fact.mod)$coefficients,2)

设计矩阵  设计是：

这个模型矩阵  设计是：

model.matr

下表显示了具备因变量的模型矩阵：

如果将 T 的列乘以均匀收益率并除以 4，则失去 T 的主效应。

预计的最小二乘系数是因子预计的二分之一，截距 β0 是样本均值。因而，因子预计是最小二乘系数的两倍。例如，

  • 4 的除数将对比度转换为两个平均值之间的差别。
  • 通过乘以各自因素的符号取得交互作用比照的符号。
  • 每列绝对于其余列齐全均衡（负数和正数相等）。
  • 均衡（正交）设计确保每个预计的成果不受其余成果的大小和符号的影响。

最小二乘预计能够在 R 中乘以 2。

fad <-lm
round(2*coeffits,2)

当有反复运行时，咱们还从回归模型中取得因子效应的 p 值和置信区间。例如，β1 的 p 值对应于温度的阶乘效应

如果原假如为真，那么   

为了取得因子效应的 95% 置信区间，咱们将回归参数的 95% 置信区间乘以 2。这在 R 中应用函数 很容易做到 confint()。

2*confint.lm

浓度主效应的 95% 置信区间为 (-8.0,-1.5)，温度和浓度之间的双向交互作用具备 95% 置信区间 (-1.46,4.96)。

因子设计绝对于一次一个因子设计的劣势

假如一次只钻研一个因素。例如，在将浓度放弃在 20% (-1) 并将催化剂放弃在 B (+1) 时钻研温度。

为了使成果具备更广泛的相关性，有必要使成果在所有其余浓度和催化剂程度上都雷同。换句话说，因素（例如，温度和催化剂）之间没有相互作用。如果成果雷同，则因子设计更无效，因为成果的预计须要更少的察看来达到雷同的精度。

如果在其余浓度和催化剂程度下成果不同，则阶乘能够检测和预计相互作用。

非反复因子设计中的正态图

回顾正态分位数图

一组数据的正态性能够通过以下办法来评估。让  示意的有序值 . 例如，r(1) 是 r1,…,rN 的最小值，r(N) 是 r1,…,rN 的最大值。所以，如果数据是：-1, 2, -10, 20，那么 。

N(0,1) 的累积散布函数 (CDF) 具备 S 形。

x <- seq
plot(x,pnorm)

因而，一组数据的正态性测验是绘制数据的有序值 r(i) 与 pi=(i-0.5)/N 的关系。如果该图与正态 CDF 具备雷同的 S 形，则这表明数据来自正态分布。

上面是从图中模仿的 1000 个随机样本的 r(i) 与 pi=(i−0.5)/N,i=1,…,N 的关系图 

N <- 1000
x <- rnorm(N)
p <- ((1:N)-0.5)/N
plot

咱们还能够构建一个正态的分位数 - 分位数图。能够证实 Φ(r(i))Φ(r(i)) 在 [0,1] 上具备均匀分布。这意味着 E(Φ(r(i)))=i/(N+1)（这是来自 [0,1] 上的均匀分布的第 i 阶统计量的期望值。

 这意味着 N 点 (pi,Φ(r(i))) 应该落在一条直线上。当初将 Φ−1 变换利用于程度和垂直尺度。N 个点

造成正态概率图 . 如果  是从正态分布生成的，而后是点图  应该是一条直线。

在 R qnorm() 中是 Φ-1。

plot(qnorm(p),sort(x))

咱们通常应用内置函数 qqnorm()（并 qqline() 增加一条直线进行比拟）来生成 QQ 图。请留神，R 应用略微更通用的分位数版本 (pi=(1−a)/(N+(1−a)−a)，其中 a=3/8，如果 N≤10，a=1/2，如果 N >10。

qqnorm(x);qqline(x)

该图与直线的显着（系统性）偏差表明：

  1. 正态假如不成立。
  2. 方差不是恒定的。

一个次要利用是在因子设计中，其中 r(i) 被有序因子效应代替。设 ^θ(1)<^θ(2)<⋯<^θ(N) 为 N 个有序因子预计。如果咱们绘制

那么靠近 0 的阶乘效应 ^θi 将沿直线降落。因而，偏离直线的点将被发表为重要点。

基本原理如下：1. 假如预计效应 ^θi 为 N(θ,σ)（预计效应波及 N 个观测值的平均值，CLT 确保 N 小至 8 的平均值靠近正态）。2. 如果 H0:θi=0,i=1,…,N 为真，那么所有预计的影响都将为零。3. 预计效应的后果正态概率图将是一条直线。4. 因而，正态概率图是测验所有预计的效应是否具备雷同的散布（即雷同的均值）。

  • 当一些效应不为零时，相应的预计效应将趋于更大并偏离直线。
  • 对于侧面影响，预计的影响落在该线之上，而负面影响落在该线之下。

示例 –  钻研化学反应的设计

一个工艺开发试验钻研了四个因素  因子设计：催化剂装料量 1、温度 2、压力 3 和其中一种反应物的浓度 4。因变量 y 是 16 个运行条件中每个条件下的转化百分比。该设计如下图所示。

该设计未反复，因而无奈预计因子效应的标准误差。

fct1 <- lm

能够取得因子效应的正态图。

Plot(fac

对应的成果 x1, x4, x2:x4, x2 不会沿着直线降落。

半正态图

相干的图形办法称为半正态概率图。让

示意无符号因子效应预计的有序值。

依据半正态分布的坐标绘制它们 – 正态随机变量的绝对值具备半正态分布。

半正态概率图由点组成

该图的一个长处是所有较大的预计效应都呈现在右上角并落在该线之上。

能够取得过程开发示例中成果的半正态图 half = TRUE。

Lenth 办法：测验没有方差预计的试验的显着性

半正态图和正态图是波及视觉判断的非正式图形办法。最好依据正式的显着性测验来定量地判断与直线的偏差。长 (1989) 提出了一种计算简略且性能良好的办法。（第 205 页，盒子、猎人和猎人（2005))

让

在 2k2k 设计 N=2k-1 中预计 θ1,θ2,…,θNθ1,θ2,…,θN 的因子效应。假如所有因子效应具备雷同的标准差。

伪标准误差 (PSE) 定义为

其中中位数是在 ∣∣^θi∣ 中计算的  和

预计的因子效应为：

ef <- 2*fat1$coeffic

s0=1.5⋅median∣∣^θi∣∣的预计是 

s0 <- 1.5*median(abs(eff))
s0

修整常数 2.5s0 是 

2.5*s0

∣∣^θi∣∣≥2.5s0 的成果 ^θi 将被修剪。上面是标记为 TRUE (x1,x2,x4,x2:x4) 的成果 

abs(eff)<2.5*s0

而后将 PSE 计算为这些值中位数的 1.5 倍。

PE <- 1.5*median
PE

ME 和 SME 是 

ME <- PE*qt
ME

PE*qt(p =(1+.95^{1/15})/2,df=(16-1)/3)

因而，成果的 95% 置信区间为：

lor <- round(ef-ME,2)
uper <- round(ef+ME,2)
kable(cbind)

具备 MEME 和 SME 的效果图通常称为 Lenth 图。PSE,ME,SMEPSE,ME,SME 的值是输入的一部分。下图中的尖峰用于显示因子效应。

Plot(fat1,cex.fac = 0.5)

该选项 cex.fac = 0.5 调整用于因子标签的字符大小。

最受欢迎的见解

1.Matlab 马尔可夫链蒙特卡罗法（MCMC）预计随机稳定率（SV，Stochastic Volatility）模型

2. 基于 R 语言的疾病制图中自适应核密度估计的阈值抉择办法

3.WinBUGS 对多元随机稳定率模型：贝叶斯预计与模型比拟

4.R 语言回归中的 hosmer-lemeshow 拟合优度测验

5.matlab 实现 MCMC 的马尔可夫切换 ARMA – GARCH 模型预计

6.R 语言区间数据回归剖析

7.R 语言 WALD 测验 VS 似然比测验

8.python 用线性回归预测股票价格

9.R 语言如何在生存剖析与 Cox 回归中计算 IDI，NRI 指标

退出移动版