全文链接：http://tecdat.cn/?p=31010

平衡不完全区组设计（BIBD）是一个很好的研究实验设计，可以从统计的角度看各种所需的特征（点击文末“阅读原文”获取完整代码数据）。

最近我们被客户要求撰写关于BIBD的研究报告，包括一些图形和统计输出。

对于一个BIBD有K个观测，重复r次实验。还有第5参数lamda，记录其中每对实验发生在设计块的数目。

生成一组BIBD设计，设计行列和每块的元素具体数目。如果BIBD(b,v,r,k)存在则：1=v

我们设置区组

BIB(7,7, 4, 2)

##      [,1] [,2] [,3] [,4]  
## [1,]    2    3    5    6  
## [2,]    3    4    6    7  
## [3,]    1    2    4    6  
## [4,]    1    5    6    7  
## [5,]    2    4    5    7  
## [6,]    1    2    3    7  
## [7,]    1    3    4    5

这种设计不是BIBD，因为处理不是所有重复的设计都有相同的次数，我们可以通过isGUID检查。对于本例：

BIB(7,7, 4, 2)

##      [,1] [,2] [,3] [,4]  
## [1,]    2    3    5    6  
## [2,]    1    5    6    7  
## [3,]    2    4    5    7  
## [4,]    1    2    4    6  
## [5,]    1    2    3    7  
## [6,]    3    4    6    7  
## [7,]    1    3    4    5

然后，我们修改参数，来查看该模型是否生产BIBD

my.design

##      [,1] [,2] [,3]  
## [1,]    1    2    6  
## [2,]    2    3    7  
## [3,]    1    4    7  
## [4,]    3    4    6  
## [5,]    1    3    5  
## [6,]    2    4    5  
## [7,]    5    6    7


##  
## [1] The design is a balanced incomplete block design w.r.t. rows.

从结果来看，该设计是一个平衡不完全区组设计。在这种情况下，我们能够生成有效BIBD实验使用指定的参数。

点击标题查阅往期内容

R语言因子实验设计nlme拟合非线性混合模型分析有机农业施氮水平

左右滑动查看更多

分析Box-Behnken设计

Box-Behnken设计的优良在于，可以将其应用于分析2至5个因子的实验。

下面将其扩展到回归模型的实验设计中，比如在下面的一个纸飞机的飞行时间的实验。这是另一个多种因子的实验，在四个变量。这些数据已经被编码。原始的变量是机翼面积A，翼状R，机身宽度W，和长度L ，在数据集中的每个观测代表的10次重复的的纸飞机在每个实验条件下的结果。我们在这里研究平均飞行时间。

使用响应曲面法对变量进行回归模型拟合

相关视频

查看模型结果

summary(heli.rsm)

##  
## Call:  
## rsm(formula = ave ~ block + SO(x1, x2, x3, x4), data = heli)  
##  
##               Estimate Std. Error  t value  Pr(>|t|)     
## (Intercept) 372.800000   1.506375 247.4815 < 2.2e-16 ***  
## block2       -2.950000   1.207787  -2.4425 0.0284522 *   
## x1           -0.083333   0.636560  -0.1309 0.8977075     
## x2            5.083333   0.636560   7.9856 1.398e-06 ***  
## x3            0.250000   0.636560   0.3927 0.7004292     
## x4           -6.083333   0.636560  -9.5566 1.633e-07 ***  
## x1:x2        -2.875000   0.779623  -3.6877 0.0024360 **  
## x1:x3        -3.750000   0.779623  -4.8100 0.0002773 ***  
## x1:x4         4.375000   0.779623   5.6117 6.412e-05 ***  
## x2:x3         4.625000   0.779623   5.9324 3.657e-05 ***  
## x2:x4        -1.500000   0.779623  -1.9240 0.0749257 .   
## x3:x4        -2.125000   0.779623  -2.7257 0.0164099 *   
## x1^2         -2.037500   0.603894  -3.3739 0.0045424 **  
## x2^2         -1.662500   0.603894  -2.7530 0.0155541 *   
## x3^2         -2.537500   0.603894  -4.2019 0.0008873 ***  
## x4^2         -0.162500   0.603894  -0.2691 0.7917877     
## ---  
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1  
##  
## Multiple R-squared:  0.9555, Adjusted R-squared:  0.9078  
## F-statistic: 20.04 on 15 and 14 DF,  p-value: 6.54e-07  
##  
## Analysis of Variance Table  
##  
## Response: ave  
##                     Df  Sum Sq Mean Sq F value    Pr(>F)  
## block                1   16.81   16.81  1.7281  0.209786  
## FO(x1, x2, x3, x4)   4 1510.00  377.50 38.8175 1.965e-07  
## TWI(x1, x2, x3, x4)  6 1114.00  185.67 19.0917 5.355e-06  
## PQ(x1, x2, x3, x4)   4  282.54   70.64  7.2634  0.002201  
## Residuals           14  136.15    9.72                   
## Lack of fit         10  125.40   12.54  4.6660  0.075500  
## Pure error           4   10.75    2.69                   
##  
## Stationary point of response surface:  
##         x1         x2         x3         x4  
##  0.8607107 -0.3307115 -0.8394866 -0.1161465  
##  
## Stationary point in original units:  
##         A         R         W         L  
## 12.916426  2.434015  1.040128  1.941927  
##  
## Eigenanalysis:  
## $values  
## [1]  3.258222 -1.198324 -3.807935 -4.651963  
##  
## $vectors  
##          [,1]       [,2]       [,3]        [,4]  
## x1  0.5177048 0.04099358  0.7608371 -0.38913772  
## x2 -0.4504231 0.58176202  0.5056034  0.45059647  
## x3 -0.4517232 0.37582195 -0.1219894 -0.79988915  
## x4  0.5701289 0.72015994 -0.3880860  0.07557783

绘制拟合值的等高线图

contour(

可视化结果

围绕拟合面，我们可以画出样本拟合点的位置。默认情况下，每个小区显示多个轮廓线的图像。可以看到，图中显示的不一定是等高线图的中心（默认可变范围是从数据中获得）;而是它设置在在坐标轴上的变量对应的值。因此，左上角的图中绘制了在x1和x2对应的拟合值，其中x3 =-0.839和x4=-0.116，在固定的值，最大的就是该坐标X1 =0.861，X2=-0.331。

点击文末“阅读原文”

获取全文完整代码数据资料。

本文选自《R语言统计学DOE实验设计：用平衡不完全区组设计（BIBD）分析纸飞机飞行时间实验数据》。

点击标题查阅往期内容

非线性混合效应 NLME模型对抗哮喘药物茶碱动力学研究

R语言用线性混合效应（多水平/层次/嵌套）模型分析声调高低与礼貌态度的关系

R语言LME4混合效应模型研究教师的受欢迎程度

R语言nlme、nlmer、lme4用（非）线性混合模型non-linear mixed model分析藻类数据实例

R语言混合线性模型、多层次模型、回归模型分析学生平均成绩GPA和可视化

R语言线性混合效应模型（固定效应&随机效应）和交互可视化3案例

R语言用lme4多层次（混合效应）广义线性模型（GLM），逻辑回归分析教育留级调查数据

R语言线性混合效应模型实战案例