R中的統(tǒng)計(jì)模型這一部分假定讀者已經(jīng)對統(tǒng)計(jì)方法，特別是回歸分析和方差分析有一定的了解。后面我們還會(huì)假定讀者對廣義線性模型和非線性模型也有所了解。R已經(jīng)很好地定義了統(tǒng)計(jì)模型擬合中的一些前提條件，因此我們能構(gòu)建出一些通用的方法以用于各種問題。R提供了一系列緊密聯(lián)系的統(tǒng)計(jì)模型擬合的工具，使得擬合工作變得簡單。正如我們在緒論中提到的一樣，基本的屏幕輸出是簡潔的，因此用戶需要調(diào)用一些輔助函數(shù)來提取細(xì)節(jié)的結(jié)果信息。1定義統(tǒng)計(jì)模型的公式下面統(tǒng)計(jì)模型的模板是一個(gè)基于獨(dú)立的方差齊性數(shù)據(jù)的線性模型yx+ee~NID(0,g2)ijijiij二0用矩陣術(shù)語表示，它可以寫成y=XP+e其中y是響應(yīng)向量，X是模型矩陣(modelmatrix)或者設(shè)計(jì)矩陣(designmatrix)oX的列x,xx是決定變量(determiningvariable)。通常，x列都是1,01p0用來定義截距(intercept)項(xiàng)。例子在給予正式的定義前，舉一些的例子可能更容易了解全貌。假定y,x,x0,X],x2,...是數(shù)值變量，X是一個(gè)矩陣，而A,B,C,...是因子。下面的例子中，左邊給出公式，右邊給出該公式的統(tǒng)計(jì)模型的描述。y~xy~1+x二者都反映了y對x的簡單線性模型。第一個(gè)公式包含了一個(gè)隱式的截距項(xiàng)，而第二個(gè)則是一個(gè)顯式的截距項(xiàng)。y~0+xy~-1+xy~x-1log(y)~x1+x2y對x過原點(diǎn)的簡單線性模型(也就是說，沒有截距項(xiàng))。y的變換形式log(y)對X]和x2進(jìn)行的多重回歸(有一個(gè)隱式的截距項(xiàng))。y~poly(x,2)y~1+x+I(xA2)y對x的二次多項(xiàng)式回歸。第一種是正交多項(xiàng)式(orthogonalpolynomial)，第二種則顯式地注明各項(xiàng)的冪次。y~X+poly(x,2)y~Ay~A+xy利用模型矩陣X和二次多項(xiàng)式項(xiàng)x進(jìn)行多重回歸。y的單因素方差分析模型，類別由A決定。y的單因素協(xié)方差分析模型，類別由A決定，協(xié)方差項(xiàng)為x。

y~A*By~A+B+A:By?B%in%Ay?A/By對A和B的非可加兩因子方差分析模型(twofactornon-additivemodel)。前兩個(gè)公式表示相同的父叉分類設(shè)計(jì)(crossedclassification),后兩個(gè)公式表示相同的嵌套分類設(shè)計(jì)(nestedclassification)o抽象一點(diǎn)說，這四個(gè)公式指明同一個(gè)模型子空間。y?(A+B+C)人2y?A*B*C-A:B:C二因子實(shí)驗(yàn)。該模型包括個(gè)主效應(yīng)(maineffects)和兩個(gè)因子的父互效應(yīng)(interactions)。這兩個(gè)公式等價(jià)。y~A*xy?A/xy?A/(1+x)-1在A的各個(gè)水平獨(dú)立擬合y對x的簡單線性回歸。三個(gè)公式的編碼不一樣。最后一個(gè)公式會(huì)對A各個(gè)水平分別估計(jì)截距項(xiàng)和斜率項(xiàng)的。y?A*B+Error(C)一個(gè)實(shí)驗(yàn)設(shè)計(jì)有兩個(gè)處理因素A和B以及因子C決定的誤差分層(errorstrata)o如在裂區(qū)實(shí)驗(yàn)設(shè)計(jì)(splitplotexperiment)中，所有區(qū)組(還包括子區(qū)組)都由因子C決定的。操作符?用來定義R的模型公式(modelformula)。一個(gè)普通的線性模型式可以表示為response?op1term1op2term2op3term3...其中response是一個(gè)作為響應(yīng)變量的向量或者矩陣，或者是一個(gè)值為向量/矩陣的表達(dá)式。opi是一個(gè)操作符。它要么是+要么是-,分別表示在一個(gè)模型中加入或者去掉某一項(xiàng)(公式第一項(xiàng)的操作符可選)。termi可以(1)是一個(gè)向量，矩陣表達(dá)式或者1,(2)因子,(3)是一個(gè)由因子,向量或矩陣通過公式操作符連接產(chǎn)生的公式表達(dá)式(formulaexpression)?；旧?，公式中的項(xiàng)決定了模型矩陣中的列要么被加入要么被去除。1表示截距項(xiàng)，并且默認(rèn)就已加入模型矩陣，除非顯式地去除這一選項(xiàng)。公式操作符(formulaoperators)在效果上和用于程序Glim和Genstat中的Wilkinson&Rogers標(biāo)記符(notation湘似。一個(gè)不可避免的改變是操作符.在R里面變成了:，因?yàn)辄c(diǎn)號(hào)在R里面是合法的命名字符。這些符號(hào)總結(jié)如下(參考Chambers&Hastie,1992,p.29):Y?MY由模型M解釋。M1+M2同時(shí)包括M1和M2項(xiàng)。M1-M2包括M1但排除M2項(xiàng)。M1:M2M1和M2的張量積(tensorproduct)o如果兩項(xiàng)都是因子，那么將產(chǎn)生“子類”因子(subclassesfactor)。M1%in%M2和M1:M2類似，但編碼方式不一樣。M1*M2M1+M2+M1:M2.M1/M2M1+M2%in%M1.MAnM的所有各項(xiàng)以及所有到n階為止的“父互作用”項(xiàng)I(M)隔離M。M內(nèi)所有操作符當(dāng)一般的運(yùn)算符處理。并且該項(xiàng)出現(xiàn)在模型矩陣中。注意，在常常用來封裝函數(shù)參數(shù)的括弧中的操作符按普通的四則運(yùn)算法則解釋。1()是一個(gè)恒等函數(shù)(identityfunction),它使得常規(guī)的算術(shù)運(yùn)算符可以用在模型公式中。還要特別注意模型公式僅僅指定了模型矩陣的列項(xiàng)，暗含了對參數(shù)項(xiàng)的指定。在某些情況下可能不是這樣，如非線性模型的參數(shù)指定。1對照我們至少要知道模型公式是如何指定模型矩陣的列項(xiàng)的。對于連續(xù)變量這是比較簡單的，因?yàn)槊恳粋€(gè)變量對應(yīng)于模型矩陣的一個(gè)列(如果模型中包含截距，會(huì)在矩陣中列出值都是1的一列)。對于一個(gè)k-水平的因子A該如何處理呢？無序和有序因子給出的結(jié)論是不一樣的。對于無序因子，因子第2,...，第k不同水平的指標(biāo)產(chǎn)生k?l列。(因此隱含的參數(shù)設(shè)置就是把其他水平和第一個(gè)水平的響應(yīng)程度進(jìn)行比較)。對于有序因子，k-1列是在1,...,k上的正交項(xiàng)(orthogonalpolynomial),并且忽略常數(shù)項(xiàng)。盡管這里的回答有點(diǎn)復(fù)雜，但這不是事情的全部。首先在含有一個(gè)因子項(xiàng)的模型中忽略截距項(xiàng)，這一項(xiàng)將會(huì)被編入指示所有因子水平的k列中。其次整個(gè)行為可以通過options設(shè)置參數(shù)contrasts而改變。R的默認(rèn)設(shè)置為options(contrasts=c("contr.treatment","contr.poly"))提這些內(nèi)容的主要原因是R和S對無序因子采用不同的默認(rèn)值。S采用Helmert對照。因此，當(dāng)你需要比較你的結(jié)果和某本書上或論文上用SPLUS代碼的結(jié)果時(shí)，你必須設(shè)置options(contrasts=c("contr.helmert","contr.poly"))這是一個(gè)經(jīng)過認(rèn)真考慮的改變。因?yàn)樘幚韺φ?treatmentcontrast)(R默認(rèn))對于新手是比較容易理解的。這還沒有結(jié)束，因?yàn)樵诟鱾€(gè)模型的各個(gè)項(xiàng)中對照方式可以用函數(shù)contrasts和C重新設(shè)置。我們還沒有考慮交互作用項(xiàng)：這些交互作用項(xiàng)將會(huì)產(chǎn)生各分量項(xiàng)的乘積。盡管細(xì)節(jié)是復(fù)雜的，R里面的模型公式在要求不是太離譜的情況下可以產(chǎn)生統(tǒng)計(jì)專家所期望的各種模型。提供模型公式的各種擴(kuò)展特性是讓R更靈活。例如，利用關(guān)聯(lián)項(xiàng)而非主要效應(yīng)的模型擬合常常會(huì)產(chǎn)生令人驚訝的結(jié)果，不過這些僅僅為統(tǒng)計(jì)專家們設(shè)計(jì)的。2線性模型對于常規(guī)的多重模型(multiplemodel)擬合，最基本的函數(shù)是lm()。下面是調(diào)用它的方式的一種改進(jìn)版：>fitted.model<-lm(formula,data=data.frame)例如>fm2<-lm(y~x1+x2,data=production)將會(huì)擬合y對x1和x2的多重回歸模型(和一個(gè)隱式的截距項(xiàng))。一個(gè)重要的(技術(shù)上可選)參數(shù)是data=production。它指定任何構(gòu)建這個(gè)模型的變量首先必須來自數(shù)據(jù)框production。這里不需要考慮數(shù)據(jù)框production是否被綁定在搜索路徑中。

3提取模型信息的泛型函數(shù)lm()的返回值是一個(gè)模型擬合結(jié)果對象；技術(shù)上就是屬于類"lm"的一個(gè)結(jié)果列表。關(guān)于擬合模型的信息可以用適合對象類,'lm"的泛型函數(shù)顯示，提取，圖示等等。這包括add1coefeffectskappapredictresidualsaliasdeviancefamilylabelsprintstepanovadrop1formulaplotprojsummary其中一些常用的泛型函數(shù)可以簡潔描述如下。anova(object1,object2)coef(object)deviance(object)formula(object)plot(object)比較一個(gè)子模型和外部模型，并且產(chǎn)生方差分析表。提取回歸系數(shù)(矩陣)。全稱：coefficients(object).殘差平方和，若有權(quán)重可加權(quán)。提取模型公式信息。產(chǎn)生四個(gè)圖，顯式殘差，擬合值和一些診斷圖。predict(object,newdata=data.frame)提供的數(shù)據(jù)框必須有同原始變量一樣標(biāo)簽的變量。結(jié)果是對應(yīng)于data.frame中決定變量預(yù)測值的向量或矩陣。predict.gam(object,newdata=data.frame)predict.gam()是安全模式的predict。。它可以用于lm,glm和gam擬合對象。在正交多項(xiàng)式作為原始的基本函數(shù)并且增加新數(shù)據(jù)意味著必須使用不同的原始基本函數(shù)。print(object)residuals(object)簡要打印一個(gè)對象的內(nèi)容。常常隱式使用。提取殘差(矩陣)，有權(quán)重時(shí)可加權(quán)，省略方式：resid(object)。step(object)通過增加或者減少模型中的項(xiàng)并且保留層次來選擇合適的模型。在逐步搜索過程中，AIC(Akaike信息規(guī)范)值最大的模型將會(huì)被返回。summary(object)顯示較詳細(xì)的模型擬合結(jié)果。4方差分析和模型比較aov(formula,data=data.frame)和函數(shù)lm()非常的相似，在泛型函數(shù)提取模型信息部分列出的泛型函數(shù)同樣適用。需要注意的是aov()還允許分析多誤差層次(multipleerrorstrata)的模型，女口裂區(qū)實(shí)驗(yàn)設(shè)計(jì)(splitplotexperiments)，利用區(qū)組內(nèi)信息進(jìn)行的平衡不完全區(qū)組設(shè)計(jì)(balancedincompleteblockdesign，等。模型公式responsemean.formula+Error(strata.formula)指定了一個(gè)多層次實(shí)驗(yàn)設(shè)計(jì)，誤差層由strata.formula定義。最簡單的情況是，strata.formula是單因素的。它定義了一個(gè)雙層次的實(shí)驗(yàn)，也就是研究在這些因子的水平內(nèi)或者水平間的實(shí)驗(yàn)響應(yīng)。例如，已知所有的決定變量因子，模型公式可以設(shè)計(jì)如下：>fm<-aov(yield~v+n*p*k+Error(farms/blocks),data=farm.data)這常常用來描述一個(gè)同時(shí)含有均值模型v+n*p*k和三個(gè)誤差層次“農(nóng)田之間”農(nóng)田內(nèi)但在區(qū)組之間”和“區(qū)組內(nèi)”)的實(shí)驗(yàn)。方差表的分析實(shí)際上是為擬合模型序列(sequence)進(jìn)行的。在模型序列的特定地方增加特定的項(xiàng)會(huì)使殘差平方和降低。因此僅僅在正交實(shí)驗(yàn)中，模型中增加項(xiàng)的次序是沒有影響的。在多層實(shí)驗(yàn)設(shè)計(jì)(multistratumexperiments)中，程序首先把響應(yīng)值依次投射到各個(gè)誤差層次上，并且用均值模型去擬合各個(gè)投射。細(xì)節(jié)內(nèi)容可以參考Chambers&Hastie(1992)。除了使用常規(guī)的方差分析表(ANOVAtable),你還可以直接用函數(shù)anova()來比較兩個(gè)模型。這種方法更為靈活。>anova(fitted.model.1,fitted.model.2,...)結(jié)果將是一個(gè)方差分析表以顯示依次加入的擬合模型的差異。需要比較的擬合模型常常是等級(jí)序列(hierarchicalsequence)。這個(gè)和默認(rèn)的設(shè)置實(shí)際上沒有差別，只是使它更容易理解和控制。5更新擬合模型函數(shù)update()是一個(gè)非常便利的函數(shù)。它允許擬合一個(gè)比原先模型增加或減少一個(gè)項(xiàng)的模型。它的形式是>new.model<-update(old.model,new.formula)在new.formula中，公式中包含的句點(diǎn)，.，僅僅表示“舊的公式模型中的對應(yīng)部分”。例如>fm05<-lm(y~x1+x2+x3+x4+x5,data=production)>fm6<-update(fm05,.~.+x6)>smf6<-update(fm6,sqrt(.)~.)這將分別擬合從數(shù)據(jù)框production中得到的五個(gè)變量的多重回歸模型，擬合額外增加一個(gè)變量的六個(gè)回歸量的模型，和進(jìn)一步對響應(yīng)值進(jìn)行平方根變換后的模型擬合。注意參數(shù)data-在最開始調(diào)用模型擬合函數(shù)的時(shí)候指定。這個(gè)信息將會(huì)通過擬合模型對象傳遞給函數(shù)update。及其相關(guān)者。符號(hào).同樣可以用在其他情況下，不過含義有點(diǎn)不同。如>fmfull<-lm(y~.,data=production)它將會(huì)擬合響應(yīng)量y對數(shù)據(jù)框production中所有變量回歸的模型。其他研究逐步回歸的函數(shù)是addl(),dropl()和step()。從字面上就可以看出這些函數(shù)的意義，更細(xì)節(jié)的內(nèi)容可以參考在線幫助文檔。6廣義線性模型廣義線性建模是線性建模的一種發(fā)展，它通過一種簡潔而又直接的方式使得線性模型既適合非正態(tài)分布的響應(yīng)值又可以進(jìn)行線性變換。廣義線性模型是基于下面一系列假設(shè)前提的：有一個(gè)有意思的響應(yīng)變量y和一系列刺激變量(stimulusvariable)x1,x2,...。這些刺激變量決定響應(yīng)變量的最終分布。刺激變量僅僅通過一個(gè)線性函數(shù)影響響應(yīng)值y的分布。該線性函數(shù)稱為線性預(yù)測器(linearpredictor)，常常寫成n邙占+卩2電++卩pXp,因此xi當(dāng)且僅當(dāng)0i=0時(shí)對y的分布沒有影響。(3)y分布的形式為fY(y；卩,8)=exp[y九(卩)-Y(九(卩))}+T(y8)]Y8其中8是尺度參數(shù)(scaleparameter)(可能已知)，對所有觀測恒定，A是一個(gè)先驗(yàn)的權(quán)重，假定知道但可能隨觀測不同有所不同，M是y的均值。也就是說假定y的分布是由均值和一個(gè)可能的尺度參數(shù)決定的。⑷均值》是線性預(yù)測器的平滑可逆函數(shù)(smoothinvertiblefunction)：M=m(n),n=m-i(卩)=1(卩)其中的反函數(shù)(inversefunction)l()被稱為關(guān)聯(lián)函數(shù)(linkfunction))這些假定比較寬松，足以包括統(tǒng)計(jì)實(shí)踐中大多數(shù)有用的統(tǒng)計(jì)模型，同時(shí)也足夠嚴(yán)謹(jǐn)，使得可以發(fā)展參數(shù)估計(jì)和統(tǒng)計(jì)推論(estimationandinference)中一致的方法(至少可以近似一致)。讀者如果想了解這方面最新的進(jìn)展，可以參考McCullagh&Nelder(1989)或者Dobson(1990)。6.1族R提供了一系列廣義線性建模工具，從類型上來說包括高斯(gaussian),二項(xiàng)式，泊松(poisson),逆高斯(inversegaussian)和伽馬(gamma)模型的響應(yīng)變量分布以及響應(yīng)變量分布無須明確給定的擬似然(quasi-likelihood)模型。在后者，方差函數(shù)(variancefunction)可以由均值的函數(shù)指定，但在其它情況下，該函數(shù)可以由響應(yīng)變量的分布得到。每一種響應(yīng)分布允許各種關(guān)聯(lián)函數(shù)將均值和線性預(yù)測器關(guān)聯(lián)起來。這些自動(dòng)可用的關(guān)聯(lián)函數(shù)如下表所示：族名字BinomialGaussianGammainverse.aussianpoissonquasi關(guān)聯(lián)函數(shù)logit,probit,log,cloglogidentity,log,inverseidentity,inverse,log1/muA2,identity,inverse,logidentity,log,sqrtlogit,probit,cloglog,identity,inverse,log,1/muA2,sqrt這些用于模型構(gòu)建過程中的響應(yīng)分布，關(guān)聯(lián)函數(shù)和各種其他必要的信息統(tǒng)稱為廣義線性模型的族(family)。6.2glm()函數(shù)既然響應(yīng)的分布僅僅通過單一的一個(gè)線性函數(shù)依賴于刺激變量，那么用于線性模型的機(jī)制同樣可以用于指定一個(gè)廣義模型的線性部分。但是族必須以一種不同的方式指定。R用于廣義線性回歸的函數(shù)是glm()，它的使用形式為>fitted.model<-glm(formula,family=family.generator,data=data.frame)和lm()相比，唯一的一個(gè)新特性就是描述族的參數(shù)family.generator。它其實(shí)是一個(gè)函數(shù)的名字，這個(gè)函數(shù)將產(chǎn)生一個(gè)函數(shù)和表達(dá)式列表用于定義和控制模型的構(gòu)建與估計(jì)過程。盡管這些內(nèi)容開始看起來有點(diǎn)復(fù)雜，但它們非常容易使用。這些名字是標(biāo)準(zhǔn)的。程序給定的族生成器可以參見族部分表格中的“族名”。當(dāng)選擇一個(gè)關(guān)聯(lián)函數(shù)時(shí)，該關(guān)聯(lián)函數(shù)名和族名可以同時(shí)在括弧里面作為參數(shù)設(shè)定。在擬(quasi)家族里面，方差函數(shù)也是以這種方式設(shè)定。一些例子可能會(huì)使這個(gè)過程更清楚。gaussian族命令>fm<-glm(y~x1+x2,family=gaussian,data=sales)和下面的命令結(jié)果一致>fm<-lm(y~x1+x2,data=sales)但是效率上，前者差一點(diǎn)。注意，高斯族沒有自動(dòng)提供關(guān)聯(lián)函數(shù)設(shè)定的選項(xiàng)，因此不允許設(shè)置參數(shù)。如一個(gè)問題需要用非標(biāo)準(zhǔn)關(guān)聯(lián)函數(shù)的高斯族，那么只能采用我們后面討論的擬族。二項(xiàng)式族考慮Silvey(1970)提供的一個(gè)人造的小例子。在Kalythos的Aegean島上，男性居民常?；加幸环N先天的眼科疾病，并且隨著年齡的增長而變的愈明顯?，F(xiàn)在搜集了各種年齡段島上男性居民的樣本，同時(shí)記錄了盲眼的數(shù)目。數(shù)據(jù)展示如下：Age:2035455570No.:tested:5050505050No.:blind:617263744我們想知道的是這些數(shù)據(jù)是否吻合logistic和probit模型，并且分別估計(jì)各個(gè)模型的LD50,也就是一個(gè)男性居民盲眼的概率為50%時(shí)候的年齡。如果y和n分別是年齡為x時(shí)的盲眼數(shù)目和檢測樣本數(shù)目，兩種模型的形式都為y?B(n,F(P°+01x))其中在probit模型中，F(xiàn)(z)=e⑵是標(biāo)準(zhǔn)的正態(tài)分布函數(shù)，而在logit模型(默認(rèn))中，F(xiàn)(z)=ez/(1+ez)。這兩種模型中，LD50=-B0/B1，即分布函數(shù)的參數(shù)為0時(shí)所在的點(diǎn)。第一步是把數(shù)據(jù)轉(zhuǎn)換成數(shù)據(jù)框,>kalythos<-data.frame(x=c(20,35,45,55,70),n=rep(50,5),y=c(6,17,26,37,44))在glm()擬合二項(xiàng)式模型時(shí)，響應(yīng)變量有三種可能性：如果響應(yīng)變量是向量，則假定操作二元(binary)數(shù)據(jù)，因此要求是0/1向量。如果響應(yīng)變量是雙列矩陣，則假定第一列為試驗(yàn)成功的次數(shù)第二列為試驗(yàn)失敗的次數(shù)。如果響應(yīng)變量是因子，則第一水平作為失敗(0)考慮而其他的作為‘成功'(1)考慮。這里，我們采用的是第二種慣例。我們在數(shù)據(jù)框中增加了一個(gè)矩陣：>kalythos$Ymat<-cbind(kalythos$y,kalythos$n-kalythos$y)為了擬合這些模型，我們采用>fmp<-glm(Ymat~x,family=binomial(link=probit),data=kalythos)>fml<-glm(Ymat~x,family=binomial,data=kalythos)既然logit的關(guān)聯(lián)函數(shù)是默認(rèn)的，因此我們可以在第二條命令中省看擬合結(jié)果，我們使用>summary(fmp)>summary(fml)兩種模型都擬合的很好。為了計(jì)算LD50，我們可以利用一個(gè)簡單>ld50<-function(b)-b[1]/b[2]>ldp<-ld50(coef(fmp));ldl<-ld50(coef(fml));c(ldp,ldl)從這些數(shù)據(jù)中得到的年齡分別是43.663年和43.601年。Poisson模型Poisson族默認(rèn)的關(guān)聯(lián)函數(shù)是log。在實(shí)際操作中，這一族常常用于擬合計(jì)數(shù)資料的Poisson對數(shù)線性模型。這些計(jì)數(shù)資料的實(shí)際分布往往符合二項(xiàng)式分布。這是一個(gè)非常重要而又龐大的話題，我們不想在這里深入展開。它甚至是非-高斯廣義模型內(nèi)容的主要部分。有時(shí)候，實(shí)踐中產(chǎn)生的Poisson數(shù)據(jù)在對數(shù)或者平方根轉(zhuǎn)化后可當(dāng)作正態(tài)數(shù)據(jù)處理。作為后者的另一種選擇是，一個(gè)Poisson廣義線性模型可以通過下面的方式擬合：>fmod<-glm(y~A+B+x,family=poisson(link=sqrt),data=worm.counts)擬似然模型對于所有的族，響應(yīng)變量的方差依賴于均值并且擁有作為乘數(shù)(multiplier)的尺度參數(shù)。方差對均值的依賴方式是響應(yīng)分布的一個(gè)特性；例如對于poisson分布Var[y]=“。對于擬似然估計(jì)和推斷，我們不是設(shè)定精確的響應(yīng)分布而是設(shè)定關(guān)聯(lián)函數(shù)和方差函數(shù)的形式，因?yàn)殛P(guān)聯(lián)函數(shù)和方差函數(shù)都依賴于均值。既然擬似然估計(jì)和gaussian分布使用的技術(shù)非常相似，因此這一族順帶提供了一種用非標(biāo)準(zhǔn)關(guān)聯(lián)函數(shù)或者方差函數(shù)擬合gaussian模型的方法。例如，考慮非線性回歸的擬合y=%zi/(z2-°2)+e同樣還可以寫成y=1/(P1x1+P2x2)+e其中x1=z2/z1,x2=-1/x1,P1=1/01andP2=02/01O假如有適合的數(shù)據(jù)框，我們可以如下進(jìn)行非線性擬合>nlfit<-glm(y~x1+x2-1,family=quasi(link=inverse,variance=constant),data=biochem)如果需要的話，讀者可以從其他手冊或者幫助文檔中得到更多的信息。7非線性最小二乘法和最大似然法模型特定形式的非線性模型可以通過廣義線性模型(glm())擬合。但是許多時(shí)候，我們必須把非線性擬合的問題作為一個(gè)非線性優(yōu)化的問題解決。R的非線性優(yōu)化程序是optim()，nlm()和nlminb()(自R2.2.0開始)。二者分別替換SPLUS的ms()和nlminb()但功能更強(qiáng)。我們通過搜尋參數(shù)值使得缺乏度(lack-of-fit)指標(biāo)最低，如nlm()就是通過循環(huán)調(diào)試各種參數(shù)值得到最優(yōu)值。和線性回歸不同，程序不一定會(huì)收斂到一個(gè)穩(wěn)定值。nlm()需要設(shè)定參數(shù)搜索的初始值，而參數(shù)估計(jì)是否收斂在很大程度上依賴于初始值設(shè)置的質(zhì)量。7.1最小二乘法擬合非線性模型的一種辦法就是使誤差平方和(SSE)或殘差平方和最小。如果觀測到的誤差極似正態(tài)分布，這種方法是非常有效的。下面是例子來自Bates&Watts(1988),51頁。具體數(shù)據(jù)是：>x<-c(0.02,0.02,0.06,0.06,0.11,0.11,0.22,0.22,0.56,0.56,1.10,1.10)>y<-c(76,47,97,107,123,139,159,152,191,201,207,200)被擬合的模型是：>fnv-function(p)sum((y-(p[1]*x)/(p[2]+x))人2)為了進(jìn)行擬合,我們需要估計(jì)參數(shù)初始值。一種尋找合理初始值的辦法把數(shù)據(jù)圖形化,然后估計(jì)一些參數(shù)值,并且利用這些值初步添加模型曲線。>plot(x,y)>xfit<-seq(.02,1.1,.05)>yfit<-200*xfit/(0.1+xfit)>lines(spline(xfit,yfit))當(dāng)然,我們可以做的更好,但是初始值200和0.1應(yīng)該足夠了?，F(xiàn)在做擬合：>out<-nlm(fn,p=c(200,0.1),hessian=TRUE)擬合后，out$minimum是誤差的平方和(SSE),out$estimate是參數(shù)的最小二乘估計(jì)值。為了得到參數(shù)估計(jì)過程中近似的標(biāo)準(zhǔn)誤(SE)，我們可以：>sqrt(diag(2*out$minimum/(length(y)-2)*solve(out$hessian)))上述命令中的2表示參數(shù)的個(gè)數(shù)。一個(gè)95%的信度區(qū)間可以通過±1.96SE計(jì)算得到。我們可以把這個(gè)最小二乘擬合曲線畫在一個(gè)新的圖上：>plot(x,y)>xfit<-seq(.02,1.1,.05)>yfit<-212.68384222*xfit/(0.06412146+xfit)>lines(spline(xfit,yfit))標(biāo)準(zhǔn)包stats提供了許多用最小二乘法擬合非線性模型的擴(kuò)充工具。我們剛剛擬合過的模型是Michaelis-Menten模型，因此可以利用下面的命令得到類似的結(jié)論。>df<-data.frame(x=x,y=y)>fit<-nls(y~SSmicmen(x,Vm,K),df)>fitNonlinearregressionmodelmodel:y~SSmicmen(x,Vm,K)data:dfVm212.68370711K0.06412123residualsum-of-squares:1195.449>summary(fit)Formula:y~SSmicmen(x,Vm,K)Parameters:EstimateStd.ErrortvaluePr(>|t|)Vm2.127e+026.947e+0030.6153.24e-11K6.412e-028.281e-037.7431.57e-05Residualstandarderror:10.93on10degreesoffreedomCorrelationofParameterEstimates:VmK0.76517.2最大似然法最大似然法(Maximumlikelihood)也是一種非線性擬合方法。它甚至可以用在誤差非正態(tài)的數(shù)據(jù)中。這種方法估計(jì)的參數(shù)將會(huì)使得對數(shù)似然值最大或者負(fù)的對數(shù)似然值最小。下面的例子來自Dobson(1990),pp.:108-111。這個(gè)例子對劑量一響應(yīng)數(shù)據(jù)擬合logistic模型(當(dāng)然也可以用glm()擬合)。數(shù)據(jù)是：>x<-c(1.6907,1.7242,1.7552,1.7842,1.8113,1.8369,1.8610,1.8839)>y<-c(6,13,18,28,52,53,61,60)>n<-c(59,60,62,56,63,59,6