分類
R

在R中執行 SAS的glm lsmeans

在R執行SAS的glm lsmeans
使用的資料為R內建的dataset- airquality
可以先看一下該資料的描述(紐約某一年的空氣品質指標數據)

?airquality
head(airquality)

輸出為csv檔供SAS使用

airquality$Month<- as.factor(airquality$Month)
write.csv(airquality, file = airquality.csv)

假設我們的實驗目的是要比較各月份的溫度是否有差異

建立線性模型
formula用於將Temp(反應變數)迴歸於Month(預測變數),使用的資料為airquality
day為類別變數# 把-1放在formula裏是不要讓截距項涵蓋在模型裏;
類別變數將自動地被設定成每個level都會有一個迴歸係數

Temp_lm<- lm(Temp~ Month -1 , data= airquality)
Temp_lm #結果顯示Intercept(截距項)和Month各類別的係數
Call:
lm(formula = Temp ~ Month - 1, data = airquality)

Coefficients:
Month5 Month6 Month7 Month8 Month9 
 65.55 79.10 83.90 83.97 76.90 

summary(Temp_lm)
#summary顯示更多模型的結果,其中包括標準誤差(Std.Error)、t 檢定值(t value)、針對回歸係數的 p 值(Pr(>|t|))、自由度(degrees of freedom)、殘差(Residual)的摘要統計和F檢定結果

Call:
lm(formula = Temp ~ Month - 1, data = airquality)

Residuals:
 Min 1Q Median 3Q Max 
-14.100 -4.548 -0.900 3.900 16.100

Coefficients:
 Estimate Std. Error t value Pr(>|t|) 
Month5 65.548 1.195 54.83 <2e-16 ***
Month6 79.100 1.215 65.09 <2e-16 ***
Month7 83.903 1.195 70.19 <2e-16 ***
Month8 83.968 1.195 70.24 <2e-16 ***
Month9 76.900 1.215 63.28 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.656 on 148 degrees of freedom
Multiple R-squared: 0.993, Adjusted R-squared: 0.9928 
F-statistic: 4221 on 5 and 148 DF, p-value: < 2.2e-16

其實會ANOVA和回歸結果都是一樣的線性模型推導出來的

執行lsmenas比較每個月溫度是否有差異
安裝執行lsmeans所需套件

install.packages(emmeans)
library(emmeans)
Temp_lsm<- lsmeans(Temp_lm, Month)
Temp_lsm

Month lsmean SE df lower.CL upper.CL
 5 65.54839 1.195454 148 63.18602 67.91075
 6 79.10000 1.215215 148 76.69859 81.50141
 7 83.90323 1.195454 148 81.54086 86.26559
 8 83.96774 1.195454 148 81.60538 86.33010
 9 76.90000 1.215215 148 74.49859 79.30141

Confidence level used: 0.95

以contrast函數執行兩兩比較各月份溫度差異,引數的interaction= T

contrast(Temp_lsm, method = pairwise, interaction = T)

 Month_pairwise estimate SE df t.ratio p.value
 5 - 6 -13.55161290 1.704657 148 -7.950 <.0001
 5 - 7 -18.35483871 1.690627 148 -10.857 <.0001
 5 - 8 -18.41935484 1.690627 148 -10.895 <.0001
 5 - 9 -11.35161290 1.704657 148 -6.659 <.0001
 6 - 7 -4.80322581 1.704657 148 -2.818 0.0055
 6 - 8 -4.86774194 1.704657 148 -2.856 0.0049
 6 - 9 2.20000000 1.718573 148 1.280 0.2025
 7 - 8 -0.06451613 1.690627 148 -0.038 0.9696
 7 - 9 7.00322581 1.704657 148 4.108 0.0001
 8 - 9 7.06774194 1.704657 148 4.146 0.0001

將Temp_lm製成圖表

library(plyr)
library(ggplot2)
TempInfo<- summary(Temp_lm)
TempCoef <- as.data.frame(TempInfo$coefficients[, 1:2])
TempCoef <- within(TempCoef, {
 Lower <- Estimate - qt(p=0.90, df=TempInfo$df[2]) * `Std. Error`
 Upper <- Estimate + qt(p=0.90, df=TempInfo$df[2]) * `Std. Error`
 Month <- rownames(TempCoef)
})
ggplot(TempCoef, aes(x=Estimate, y=Month)) + geom_point() +
 geom_errorbarh(aes(xmin=Lower, xmax=Upper), height=.3) +
 ggtitle(Temp by Month calculated from regression model)
![SAS_glm_lsmeans.png](https://pic.pimg.tw/alex59638/1538829973-658346672.png "SAS_glm_lsmeans.png")

現在換成在SAS中執行glm中的lsmeans語法

PROC IMPORT
DATAFILE= C:\\R\\airquality.csv /* 要匯入資料的資料位置 */
OUT= Course.airquality /* 儲存到SAS的哪個位置*/ 
DBMS= CSV /* 匯入檔案的格式*/ 
REPLACE; /*宣告匯入的資料覆寫現有的資料集,系統預設是不覆寫,如不宣告則可能導致匯入資料重複*/
GETNAMES=YES; /*宣告是否將資料的第一行當作變項名稱*/
run;
/*匯入資料方式的參考來源 https://blog.timshan.idv.tw/2013/06/howtoimportsas.html */

proc print data = Course.airquality;
run;

proc glm data = Course.airquality;
class Month;
model Temp = Month ;
lsmeans Month/ stderr pdiff;
run;

下面這個網誌簡單的說明lsmeans的功能,使各組間means的權重一致,達到修正模型的結果

[SAS] 多重比較 based on “Least Squares Means (lsmeans)

SAS得到的結果如下

1538829823571.jpg

1538829837687.jpg

BoxPlot.png

後記: 其實我是在先學SAS的glm語法,所以一直以為這種統計方法是glm(廣義線性回歸模型)

後來才發現這只是SAS的語法名稱而已,實際上的統計方式還是簡單線性回歸模型而已

想說明明比較的只有各類別中的數值而已,怎麼會需要用到glm

附上紀錄檔

R script

SAS

發表迴響