R語言實現 廣義加性模型 Generalized Additive Models(GAM) 入門

轉載請說明。

R語言官網：http://www.r-project.org/

R語言軟件下載：http://ftp.ctex.org/mirrors/CRAN/         注：下載時點擊 install R for the first time

下面進行一個簡單的入門程序學習。

先新建一個txt，叫做 Rice_insect.txt 點我下載，內容為：（用制表符Tab）

Year	Adult	Day	Precipitation
1973	27285	15	387.3
1974	239	14	126.3
1975	6164	11	165.9
1976	2535	24	184.9
1977	4875	30	166.9
1978	9564	24	146.0
1979	263	3	24.0
1980	3600	21	23.0
1981	21225	13	167.0
1982	915	12	67.0
1983	225	17	307.0
1984	240	40	295.0
1985	5055	25	266.0
1986	4095	15	115.0
1987	1875	21	140.0
1988	12810	32	369.0
1989	5850	21	167.0
1990	4260	39	270.8

　　Adult為累計蛾量，Day為降雨持續天數，Precipitation為降雨量。

輸入代碼：

library(mgcv)      #加載mgcv軟件包，因為gam函數在這個包里
Data <- read.delim("Rice_insect.txt")     #讀取txt數據，存到Data變量中
Data <- as.matrix(Data)     #轉為矩陣形式
#查看Data數據：Data，查看第2列：Data[,2]，第2行：Data[2,]

Adult<-Data[,2]
Day<-Data[,3]
Precipitation<-Data[,4]

result1 <- gam(log(Adult) ~ s(Day))     #此時，Adult為相應變量，Day為解釋變量
summary(result1)      #輸出計算結果

　　此時可以看到：

 

Family: gaussian
Link function: identity

Formula:
log(Adult) ~ s(Day)

Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.9013 0.3562 22.18 4.83e-13 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
edf Ref.df F p-value
s(Day) 1.713 2.139 0.797 0.473

R-sq.(adj) = 0.0471 Deviance explained = 14.3%
GCV score = 2.6898 Scale est. = 2.2844 n = 18

Day的影響水平p-value=0.473，解釋能力為14.3%，說明影響不明顯。

其中，edf為自由度，理論上，當自由度接近1時，表示是線性關系；當自由度比1大，則表示為曲線關系。

接下來使用另一個解釋變量Precipitation。輸入代碼：

result2 <- gam(log(Adult) ~ s(Precipitation)) 
summary(result2) 

　　發現：

Family: gaussian
Link function: identity

Formula:
log(Adult) ~ s(Precipitation)

Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.9013 0.2731 28.94 1.87e-12 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
edf Ref.df F p-value
s(Precipitation) 5.022 6.032 2.561 0.0774 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) = 0.44 Deviance explained = 60.6%
GCV score = 2.0168 Scale est. = 1.342 n = 18

此時p-value為0.0774，說明該因子在P<0.1水平下影響顯著。（一般情況下p<0.05為顯著。）

接下來畫圖：

plot(result2,se=T,resid=T,pch=16)

　　

pch=16這個是圖標的序號，比如改成17就是三角形了。

 

log(Adult)中的log是什么意思呢?

log是數據變換，取對數可以把大范圍的數變成小范圍的數，這在將幾組相差太大的數據畫在同一個坐標軸時特別有用，比如一組數據范圍是1～10，第二組數據范圍是10～100000000，要是不對第二組取常用對數，第一組在坐標軸上只是一點點，都看不到，對第二組取常用對數后，第二組范圍變成1～8了，這樣兩組數據都能看到了。

 

下面嘗試將兩個變量同時作為解釋變量。

result3<-gam(log(Adult)~s(Precipitation)+s(Day))

　　出錯：Model has more coefficients than data

解決辦法：改成：

result3<-gam(log(Adult)~s(Precipitation,k=9)+s(Day,k=9))

　　k是什么？

　　k is the dimension of the basis used to represent the smooth term. If k is not specified then basis specific defaults are used.

K的最小值是3，最大值是17，為3、4的時候都是直線，說明太小了體現不出來，在不斷增大的過程中發現，K越大，曲線原來越平滑，再大時，曲線就出現了一些彎曲，說明更精准了，再大時，圖形就基本不變了，說明也沒必要設那么大了。

再

summary(result3)

　　結果：

Family: gaussian
Link function: identity

Formula:
log(Adult) ~ s(Precipitation, k = 9) + s(Day, k = 9)

Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.9013 0.2831 27.91 8.16e-12 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
edf Ref.df F p-value
s(Precipitation) 4.653 5.572 2.546 0.0848 .
s(Day) 1.000 1.000 0.500 0.4939
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) = 0.398 Deviance explained = 59.8%
GCV score = 2.288 Scale est. = 1.4423 n = 18

 

 

 

 

 

 

 

 

 

 

R語言基礎包函數中文幫助文檔(中英文對照v10)   http://www.docin.com/p-585638419.html

R語言:formula.gam()函數中文幫助文檔    http://www.biostatistic.net/thread-6589-1-45.html

R語言 mgcv包 gam()函數中文幫助文檔    http://www.biostatistic.net/thread-56086-1-1.html