试验设计与数据处理第六章续(第七,八节)

第七节二次响应面分析法当试验中考察的指标宜于用多元二次回归方程来拟合因素与指标的函数关系，就可以分析回归方程所反应的曲面形状，如果得到的曲面是凸面(像山丘)或凹面(像山谷)这类简单曲面，那么预测的最佳指标值(极大值或极小值)可以从所估计的曲面上获得；如果曲面很复杂，或者预测的最佳点远离所考察因素的试验范围，那么可以通过岭嵴分析来确定重新进行试验的方向.这就是应用较广，颇有实用价值的响应面分析法(ResponseSurfaceAnalysis).响应面分析的三个任务:任务1.模型拟合和方差分析;任务2.典型分析(寻找最优值点);任务3.岭嵴分析(若最优值点在试验范围内达到,则不必进行岭嵴分析;否则就要进行).例6.7.2求面包酵母高产和高发酵活力的最优工艺参数.指标:y1---面包酵母产量[指最终发酵液的菌体浓度(g/L)]y2---发酵活力[指面团从水中浮起的时间(min)]因素:z1---还原糖浓度,2%~6%;z2---氮源浓度,0.4%~1.2%;z3---磷源浓度,0.2%~0.4%.取三水平,z1:2%,4%,6%;z2:0.4%,0.8%,1.2%;z3:0.2%,0.3%,0.4%.令编码变量xi=[zi-(zi1+zi2)/2]/(zi2-zi1)/2故可作变换:X1=(z1-4)/2X2=(z2-0.8)/0.4X3=(z3-0.3)/0.1响应面分析的SAS程序如下:DataE62;Inputz1-z3x1-x3y1y2@@;Cards;

40.40.20-1-124.513.9

40.40.40-1135.314.1………40.80.300037.58.1

;Proc

RsReg

data=E62;/*响应面分析*/Modely1y2=z1-z3;Run;

TheRSREGProcedureCodingCoefficientsfortheIndependentVariablesFactorSubtractedoffDividedbyz14.0000002.000000z20.8000000.400000z30.3000000.100000ResponseSurfaceforVariabley1ResponseMean32.400000RootMSE1.329474R-Square0.9775CoefficientofVariation4.1033下面给出程序运行得到的部分结果:

TypeISumRegressionDFofSquaresR-SquareFValuePr>FLinear314.5575000.03712.750.1524Quadratic3116.5325000.297221.980.0026Crossproduct3252.1725000.643147.560.0004TotalModel9383.2625000.977524.090.0013SumofResidualDFSquaresMeanSquareTotalError58.8375001.767500

TheRSREGProcedureCanonicalAnalysisofResponseSurfaceBasedonCodedDataCriticalValueFactorCodedUncodedz1-0.0988443.802311z20.2648400.905936z3-0.3158000.268420Predictedvalueatstationarypoint:38.040882EigenvectorsEigenvaluesz1z2z3-0.5185070.184252-0.7614870.621441-1.7733340.894929-0.131446-0.426408-8.0206590.4063900.6347120.657258Stationarypointisamaximum.如果稳定点不是理想点就要进一步作岭嵴分析,请看示意图和例子演示(例见datarubber或datapectic)例6.7.3橡胶树幼苗对硫35的吸收实验(2×3析因试验)因素:温度t(℃)---7℃,20℃

取样时间d---2d,7d,12d响应面分析SAS程序如下:datarubber;inputftd@@;cards;数(略);procsort;bytd;/*对自变量td进行sort由小到大排序*/procrsreg;modelf=td/lackfit;/*选项lackfit要求对有重复试验自由度的回归模型执行不适合度检定(lack-of-fittest)(误差参照为重复试验的均方误差),预先应先对自变量进行sort由小到大排序*/run;可以作二个因素的响应面图(固定其它因素),E62的响应面图如下:回归方程作响应面,在对数据集E62进行响应面回归(Rsreg)后得回归方程(固定水平x3=0,本例自变量均已经编码处理)y1=37.9-0.5625x1-0.6x2-2.875x12-3.575x1x2-3.45x22SAS程序如下:

DataE62G;Dox1=-1to1by0.1;/*当自变量是原始变量时取类似的二元网格点,即do变量=左端点

to右端点

by步长*/Dox2=-1to1by0.1;/*当自变量是原始变量时取类似的二元网格点,即do变量=左端点

to右端点

by步长*/Y1=37.9-0.5625*x1-0.6*x2-2.875*x1**2-3.575*x1*x2-3.45*x2**2;Output;End;End;Procprint;Run;/*以上为由回归方程建立作图数据集*/等值线图(contourplot)37.95/*以下利用程序作三维图*/procg3ddata=E62G;plotX2*X1=Y1;/*作x1(纵)-x2(横)-y1(竖)三维图/run;第八节非线性回归非线性回归是指模型从参数角度为非线性的且无法线性化的模型,例如y=b0+b1x1+b2x2+b11x12+b12x1x2+b22x22是诸参数bi及bij的线性模型(从自变量x1,x2角度模型是非线性的),y=aebx是可以线性化的非线性模型(线性化后lny=bx+lna),而y=aexp(-exp(b-cx))是非线性的且无法线性化的模型.此时SAS采用nlin非线性回归程序程序形式:procnlindata=…;model因变量

=非线性函数表达式;

parms

参数1=初值参数2=初值

…;

der.参数1=偏导表达式;der.参数2=偏导表达式;……………;runs;方法:采用迭代求参数,要给出参数的一阶偏导表达式der.参数=…;,要赋参数初值parms参数=…;请看例子演示例6.8.1已知牧草重量y与生长天数x的关系是y=aexp(-exp(b-cx))+ε.9次观察的数据如表6.8.1,求a,b,c的值.datahw;inputxy;cards;数据(略)

；procnlin

data=hw;/*调用nlin进行非线性回归*/zz=exp(b-c*x);/*简化记号*/term=exp(-zz);/*简化记号*/model

y=a*term;

/*因变量=非线性函数表达式,即y=aexp(-exp(b-cx))*/parmsa=70b=1.1c=0.2;/*赋参数初值*/der.a=term;/*给出参数的一阶偏导表达式*/der.b=-a*zz*term;der.c=a*x*zz*term;run;

DependentVariableyMethod:Gauss-NewtonSumofIterbcaSquares01.10000.200070.00009161.21-0.1190