在研究導(dǎo)師類型和性別這個兩個自變量學(xué)生GPA的影響時，我們想看一看導(dǎo)師類型和性別的交互作用，這個例子中，導(dǎo)師類型和性別都是分類變量。

數(shù)據(jù)模擬

N <- 250

Q <- sample(rep(c(-1,0,1),N),N,replace = FALSE)

G <- sample(rep(c(0,1),N*3/2),N,replace = FALSE)

Y <- .5*Q + .25*G + 2.5*Q*G+ 1 + rnorm(N, sd=2)

Y = (Y - min(Y)) / (max(Y) - min(Y))*4

GPA.Data.3<-data.frame(GPA=Y,Tutor=Q,Gender=G)

上面的代碼就模擬出了個數(shù)據(jù)框，其中Q為導(dǎo)師類型，有3類，模擬的時候用的-1，0，1代替，性別用的0，1

接下來，我們得將分類變量轉(zhuǎn)化為啞變量，R會自動將因子處理為啞變量，所以我們直接轉(zhuǎn)為因子就行：

GPA.Data.3$Tutor.F <- factor(GPA.Data.3$Tutor,

level=c(-1,0,1),

labels=c("No Tutor", "Group Tutor", "Private Tutor"))

GPA.Data.3$Gender.F <- factor(GPA.Data.3$Gender,

level=c(0,1),

labels=c("Male", "Female"))

回歸分析

然后我們進(jìn)行回歸分析，一個是有交互作用的，一個是不含交互作用的

GPA.3.Model.1<-lm(GPA ~ Tutor.F+Gender.F, data = GPA.Data.3)

GPA.3.Model.2<-lm(GPA ~ Tutor.F*Gender.F, data = GPA.Data.3)

stargazer(GPA.3.Model.1, GPA.3.Model.2,type="html",

column.labels = c("Main Effects", "Interaction"),

intercept.bottom = FALSE,

single.row=TRUE,

notes.append = FALSE,

omit.stat=c("ser"),

star.cutoffs = c(0.05, 0.01, 0.001),

header=FALSE，

out="test.html")

運(yùn)行代碼得出回歸分析的結(jié)果如下：

畫簡單斜率圖

畫斜率圖的第一步依然是跑出每個子類的系數(shù)：

Inter.GPA.3 <- effect('Tutor.F*Gender.F', GPA.3.Model.2,

se=TRUE)

Inter.GPA.3.DF<-as.data.frame(Inter.GPA.3)

我們可以看一看每個子類的擬合系數(shù)長啥樣：

上圖就顯示了，在兩個自變量的不同水平fit值是不一樣的，而且通過回歸的交互項(xiàng)顯著性分析，也是知道這個fit是有顯著性差異的

接下來就開始畫圖：

Inter.GPA.3.DF$Tutor.F <- factor(Inter.GPA.3.DF$Tutor,

level=c("No Tutor", "Group Tutor", "Private Tutor"),

labels=c("No Tutor", "Group Tutor", "Private Tutor"))

Inter.GPA.3.DF$Gender.F <- factor(Inter.GPA.3.DF$Gender,

level=c("Male", "Female"),

labels=c("Male", "Female"))

Plot.GPA.3<-ggplot(data=Inter.GPA.3.DF, aes(x=Tutor.F, y=fit, group=Gender.F))+

geom_line(size=2, aes(color=Gender.F))+

geom_ribbon(aes(ymin=fit-se, ymax=fit+se,fill=Gender.F),alpha=.2)+

ylab("GPA")+

xlab("Tutor")+

ggtitle("Tutors and Gender as GPA Predictors")+ theme_bw()+

theme(text = element_text(size=12),

legend.text = element_text(size=12),

legend.direction = "horizontal",

panel.grid.major = element_blank(),

panel.grid.minor = element_blank(),

legend.position="top")

Plot.GPA.3

到這兒，我們的簡單斜率圖就畫好了，注意兩自變量都是分類變量的時候，x軸上的自變量有3個水平，所以這個斜率圖一般都不是直的，一般斜率在中間都會打個折。