dta <- read.csv("ncku_prof_V6.csv", h=T, stringsAsFactors = TRUE)
dta <- dta %>%
  select(H.id, Gender, Degree, Rank, College, Dept, Grads, FPY, Articles) %>%
  mutate(researchy = 2022 - FPY)

Degree

## 檢驗學位國籍是否是研究生涯年與論文引數的調節變項
fullmod <- lm(H.id ~ researchy + Degree + Degree*researchy, data = dta)
reducemod <- lm(H.id ~ researchy + Degree, data = dta)
anova(fullmod, reducemod)
## Analysis of Variance Table
## 
## Model 1: H.id ~ researchy + Degree + Degree * researchy
## Model 2: H.id ~ researchy + Degree
##   Res.Df   RSS Df Sum of Sq      F Pr(>F)
## 1    454 49809                           
## 2    455 49863 -1   -54.377 0.4956 0.4818

根據統計檢定結果，我們可以看見兩個模型間的差異並未達顯著，因此學位國籍可能與研究生涯年與論文引數之間沒有交互作用(沒有調節效應)。

summary(fullmod)$r.sq - summary(reducemod)$r.sq
## [1] 0.0008148254

兩個模型的解釋量差異為0.0008148254

## 調節效果圖
ggplot(aes(y = H.id, x = researchy, color = Degree), data = dta) +
 geom_smooth(method = "lm", se = F, fullrange = T) +
  theme_bw()
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 2 rows containing non-finite values (stat_smooth).

從這個圖我們可以看見，雖然在研究生涯年較少時，本國籍博士畢業的教職員，其論文指數稍微高一點，但當研究生涯年落在35年左右時，兩者的論文指數會很接近，而國外畢業的教職員在研究生涯較久時，其論文指數可能會比本國籍博士畢業的教職員來的高一點點。

而且這兩條線雖有交點，但因為其斜率看起來很接近，因此可能學位國籍對此兩個變項真的沒有什麼調節效果。

## 檢驗簡單斜率是否顯著
dta$nDegree <- as.numeric(dta$Degree)
modmod <- lmres(H.id ~ researchy + nDegree + researchy*nDegree, dta)
summary(modmod)
## Formula:
## H.id ~ researchy + nDegree + researchy.XX.nDegree
## <environment: 0x000001f7d584b868>
## 
## Models
##          R     R^2   Adj. R^2    F     df1  df2  p.value    
## Model  0.504  0.254     0.249 51.426  3.000  454  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residuals
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -23.5372  -6.2005  -0.7767   0.0000   4.9457  70.9566 
## 
## Coefficients
##                      Estimate   StdErr  t.value    beta p.value  
## (Intercept)           4.87999  4.81851  1.01276         0.31171  
## researchy             0.52852  0.22973  2.30068  0.3899 0.02186 *
## nDegree              -3.60111  2.72920 -1.31947 -0.1369 0.18768  
## researchy.XX.nDegree  0.08976  0.12750  0.70402  0.1430 0.48178  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Collinearity
##                          VIF Tolerance
## researchy            17.4680    0.0572
## nDegree               6.5461    0.1528
## researchy.XX.nDegree 25.0876    0.0399
ss <- simpleSlope(modmod, pred = "researchy", mod1 = "nDegree", coded = "nDegree")
summary(ss)
## 
## ** Estimated points of H.id  **
## 
##                    Low researchy (-1 SD) High researchy (+1 SD)
## Low nDegree ( 1 )                  8.526                 19.549
## High nDegree ( 2 )                 5.977                 18.600
## 
## 
## 
## ** Simple Slopes analysis ( df= 454 ) **
## 
##                    simple slope standard error t-value p.value    
## Low nDegree ( 1 )        0.6183         0.1103     5.6  <2e-16 ***
## High nDegree ( 2 )       0.7080         0.0639    11.1  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## 
## ** Bauer & Curran 95% CI **
## 
##         lower CI upper CI
## nDegree  -0.2349   5.9949

於此統計檢定結果中，我們可以看見不論是本國籍畢業的教職員抑或是國外畢業的教職員，其斜率相較於0皆為顯著。

## 調節效果圖
PlotSlope(ss)