---
title: "การประเมินลักษณะมนุษย์: ตัวอย่าง R Code สำหรับคาบเรียนที่ 5"
author: "สันทัด พรประเสริฐมานิต"
output: html_document
date: "`r format(Sys.time(), '%d/%m/%Y')`"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r xx, include=FALSE}
library(vegan)
library(psych)
library(mvtnorm)
```

## การวิเคราะห์ความซ้ำซ้อน (Redundancy Analysis; RDA)

อ่านข้อมูลด้วยคำสั่ง `read.table` แล้วใส่ไฟล์ที่ต้องการอ่านลงไป ในที่นี้คือ ` leadership.csv` 

```{r openfile}
dat <- read.table("leadership.csv", sep=",", header=TRUE)
```

ข้อมูลนี้เป็นข้อมูลจากลูกน้องจำนวน 201 คน ประเมินหัวหน้าตนเองด้วยมาตร Multifactor Leadership Questionnaire (MLQ) โดยคะแนนจากมาตรนี้จะแบ่งออกเป็นลักษณะของผู้นำ 9 ด้าน ซึ่งจะใช้เป็นตัวแปรทำนาย และตัวแปรผลลัพธ์จะเป็นการประเมินตนเองของลูกน้อง โดยจะประเมินเรื่องประสิทธิภาพในการทำงาน การใช้พยายามพิเศษในการทำงาน และความพึงพอใจในงาน ดูค่าสถิติพื้นฐานของข้อมูลได้โดยคำสั่ง `describe` จาก `psych` package 

```{r describe}
library(psych)
describe(dat)
```

ในการทำ RDA จะใช้คำสั่ง `rda` ในแพ๊คเกจ `vegan` แต่การวิเคราะห์ RDA จะเกี่ยวข้องกับตัวแปรจำนวนมาก เพื่อให้ง่ายต่อการวิเคราะห์ ให้สร้างข้อมูลใหม่ ที่ข้อมูลหนึ่งมีเฉพาะตัวแปรทำนาย และอีกข้อมูลหนึ่งมีเฉพาะตัวแปรผลลัพธ์

```{r sepdata}
predictors <- dat[,c("idealattributed", "idealbehavior", 
                     "individualconsideration", "inspiremotive", 
                     "intelstimulation", "contingentreward", 
                     "managebyexceptactive", "managebyexceptpassive", 
                     "laissezfaire")]
outcomes <- dat[,c("effectiveness", "extraeffort", "satisfaction")]
```

จากนั้น ใช้ฟังก์ชั่น `rda` โดยใส่ข้อมูลผลลัพธ์ (`outcomes`) ไว้ทางซ้ายของ `~`  และใส่ `.` ทางด้านขวา เพื่อแสดงว่าใช้ทุกตัวแปรในข้อมูลในการทำนาย แล้วใส่ข้อมูลตัวแปรต้นใน `data` จากนั้นใช้คำสั่ง `scale=TRUE` เพื่อเปลี่ยนตัวแปรตามให้เป็นคะแนนมาตรฐาน จะทำให้น้ำหนักของตัวแปรตามแต่ละตัวเท่ากัน ในการแสดงผลจะใช้คำสั่ง `summary` แต่เนื่องจากข้อมูลผลลัพธ์ที่ให้มีจำนวนมาก จึงใส่ `display=c("sp", "bp", "cn")` เพื่อคัดเลือกเฉพาะผลลัพธ์ที่สำคัญ

```{r rda}
library(vegan)
out <- rda(outcomes ~ ., data=predictors, scale=TRUE)
summary(out, display=c("sp", "bp", "cn"))
```

ส่วนต่างๆ ของผลลัพธ์เป็นดังนี้

- `Partitioning of correlations: ` เป็นการแบ่งความแปรปรวนของตัวแปรตาม ว่าทั้งหมดมีสัดส่วน 100% (เทียบเท่ากับตัวแปร 3 ตัว) แบ่งเป็นส่วนที่อธิบายได้ 48.3% (เทียบเท่ากับตัวแปร 1.449 ตัว) และส่วนที่อธิบายไม่ได้ 51.7% (เทียบเท่ากับ 1.551 ตัว)
- `Eigenvalues, and their contribution to the correlations` จะบอกค่า Eigenvalues ของผลรวมแต่ละตัว จำนวนผลรวมที่เป็นไปได้จะเท่ากับจำนวนตัวแปรตาม 3 ตัว ผลรวมทั้ง 3 ตัวมีค่า Eigenvalues คือ 1.404, 0.024, และ 0.021 เทียบเท่ากับสัดส่วนความแปรปรวนรวมของตัวแปรตามทั้งหมดเท่ากับ 46.8%, 0.8% และ 0.7% ตามลำดับ และหากมองเฉพาะความแปรปรวนที่อธิบายได้ ผลรวมแต่ละตัวคิดเป็นสัดส่วน 96.9%, 1.7% และ 1.4% ตามลำดับ
- `Species scores` จะบอกความหมายของผลรวมแต่ละตัว เมื่อเทียบกับตัวแปรตาม เมื่อผลรวมที่ 1 จะแสดงให้เห็นว่าเมื่อผลรวมนี้มีค่าสูงขึ้น ค่าของตัวแปรตามทั้ง 3 ตัวมีค่าน้อยลง เมื่อตรวจดูผลรวมที่ 2 จะพบว่าเมื่อผลรวมมีค่าสูง ค่าของประสิทธิภาพสูงขึ้น แต่ค่าความพยายามพิเศษ และความพึงพอใจจะน้อยลง และเมื่อดูผลรวมที่ 3 จะพบว่าเมื่อผลรวมมีค่าสูง ค่าของประสิทธิภาพ และความพึงพอใจมีค่าสูงขึ้น แต่ค่าความพยายามพิเศษจะมีค่าน้อยลง
- `Biplot scores for constraining variables` จะบ่งบอกค่าสหสัมพันธ์ระหว่างผลรวมทั้งสาม กับตัวแปรทำนายแต่ละตัว (คล้ายกับความหมายน้ำหนักองค์ประกอบ) ค่าสหสัมพันธ์เหล่านี้จะช่วยในการแปลความหมายผลรวมแต่ละตัว ในตัวอย่างนี้ผลรวมที่ 1 จะหมายถึงลักษณะของผู้นำในผู้นำที่ส่งผลทางบวกกับลูกน้อง 

สามารถหาค่าสัมประสิทธิ์การทำนายที่ปรับแล้ว (Adjusted $R^2$) ได้ดังนี้

```{r rsquared}
RsquareAdj(out)
```

การทดสอบนัยสำคัญสามารถตรวจสอบได้ทั้งหมด 3 รูปแบบ แบบแรกคือทดสอบว่าตัวแปรอิสระทั้งหมด สามารถสร้างผลรวมใดผลรวมหนึ่งแล้วทำนายตัวแปรตามได้อย่างมีนัยสำคัญหรือไม่ ซึ่งการทดสอบนี้พบผลนัยสำคัญ

```{r anovarda}
anova.cca(out)
```

ต่อไป สามารถตรวจสอบได้ว่าผลรวมแต่ละตัว สามารถอธิบายตัวแปรตามได้อย่างมีนัยสำคัญหรือไม่ พบว่ามีเพียงผลรวมแรกที่ถึงระดับนัยสำคัญ

```{r anovaeachrda}
anova.cca(out, by="axis")
```

นอกจากนี้ สามารถตรวจสอบได้ว่า ตัวแปรอิสระแต่ละตัว หากลดออกจากโมเดลแล้ว จะทำให้การทำนายตัวแปรตามเปลี่ยนแปลงไปอย่างมีนัยสำคัญหรือไม่ ซึ่งถึงระดับนัยสำคัญ แสดงว่าตัวแปรทำนายดังกล่าว สามารถทำนายตัวแปรตามได้เพิ่มเติมจากตัวแปรอื่น การแปลความหมายคล้ายกับสัมประสิทธิ์ถดถอย

```{r anovaeachpredictor}
anova.cca(out, by="margin")
```

นอกจากนี้ สามารถนำค่า Eigenvalues ของผลรวมต่างๆ ไปสร้างเป็นกราฟเพื่อตรวจสอบการตกลงของค่า Eigenvalues จากผลรวมลำดับต่างๆ เรียกว่า Scree Plot การแปลความหมายคล้ายๆ กับการวิเคราะห์องค์ประกอบเชิงสำรวจ (Exploratory Factor Analysis) ที่ผลรวมท้ายๆ จะเป็นค่าที่ทำนายตัวแปรเดียวโดดๆ แต่ผลรวมแรกๆ จะเป็นการทำนายตัวแปรมากกว่า 1 ตัว เมื่อไรก็ตามที่เห็นการตกลงของค่า Eigenvalues ในกราฟ เหมือนการตกหน้าผา ด้านซ้ายของหน้าผา คือจำนวนผลรวมที่ควรเก็บไว้

```{r screeplot}
screeplot(out)
```

จากข้อมูลนี้ พบว่าผลรวมแรกเพียงอย่างเดียวน่าจะเพียงพอในการอธิบายความสัมพันธ์ระหว่างตัวแปรทำนายและตัวแปรผลลัพธ์ ซึ่งหากอ่านผลรวมกับการทดสอบนัยสำคัญของผลรวมแต่ละตัว พบว่ามีเพียงผลรวมแรกที่ถึงระดับนัยสำคัญ และผลรวมแรกเป็นเพียงตัวเดียวที่มีความสามารถในการทำนายตัวแปรตามเกิน 1% จากเหตุผลเหล่านี้ จึงสนใจเพียงเฉพาะผลรวมแรกเท่านั้น

ในการสร้างค่าที่ทำนายได้ (Predicted Values) โดยใช้คำสั่ง `predict` สามารถสร้างค่าทำนายได้ 2 แบบ แบบแรกคือค่าของตัวแปรตามที่ทำนายได้จากผลรวม 3 ตัว ใส่ `type="response"` เพื่อหาค่าที่ทำนายได้ของตัวแปรตาม

```{r predictresponses3}
predictedresponse <- predict(out, type="response")
head(predictedresponse)
```

อย่างไรก็ตาม ผลรวมที่น่าสนใจมีเพียงผลรวมแรกเพียงอย่างเดียว ดังนั้นจึงอาจเลือกหาค่าที่ทำนายได้ของตัวแปรตาม เมื่อใช้ผลรวมเพียงแค่ตัวแปรเดียว โดยใส่ `rank=1` 

```{r predictresponses1}
predictedresponse2 <- predict(out, type="response", rank=1)
head(predictedresponse2)
```

นอกจากนี้ ยังสามารถหาของผลรวม ซึ่งก็คือค่าดัชนี ในการวัดแบบก่อรูป (Formative Measurement) โดยใส่ `type="lc"` เพื่อหาค่าที่ทำนายได้ของตัวแปรตาม

```{r predictlc}
predictedlc <- predict(out, type="lc")
head(predictedlc)
```

สมมติว่ามีค่าของตัวแปรทำนายใหม่ 2 ชุด ชุดแรก ตัวแปรทำนาย 6 ตัวแรกมีค่าสูงสุด (20) และ 3 ตัวหลังมีค่าต่ำสุด (4) ส่วนชุดที่สอง ตัวแปรทำนาย 6 ตัวแรกมีค่าต่ำสุด (4) และ 3 ตัวหลังมีค่าสูงสุด (20) สร้างชุดข้อมูลในรูปแบบ `data.frame`

 ```{r newdata}
predictvarnames <- c("idealattributed", "idealbehavior", 
                     "individualconsideration", "inspiremotive", 
                     "intelstimulation", "contingentreward", 
                     "managebyexceptactive", "managebyexceptpassive", 
                     "laissezfaire")
predictorvalues <- c(20, 20, 20, 20, 20, 20, 4, 4, 4)
predictorvalues2 <- c(4, 4, 4, 4, 4, 4, 20, 20, 20)
predval <- rbind(predictorvalues, predictorvalues2)
colnames(predval) <- predictvarnames
predval <- data.frame(predval)
```

จากนั้น นำข้อมูลใหม่ที่ได้มา เอาไปสร้างค่าที่ทำนายได้ 

```{r predictnewdata}
predict(out, type="response", newdata=predval)
predict(out, type="response", newdata=predval, rank=1)
predict(out, type="lc", newdata=predval)
```

สุดท้ายนี้ นักวิจัยสามารถที่จะลดจำนวนตัวแปรทำนายได้ เหมือนกับการวิเคราะห์ถดถอย ในที่นี้อาจใช้วิธี Stepwise ในการคัดเลือกตัวแปรทำนาย สามารถทำได้ผ่านคำสั่ง `ordiR2step` โดยค่าแรก ให้ใส่โมเดลเริ่มต้น ซึ่งเป็นโมเดลที่มีแต่จุดตัด ไม่มีตัวแปรต้นเลย (ใส่แค่ 1 ในการทำนาย) ส่วนค่าที่สอง `scope` คือโมเดลสูงที่สุดที่ทำได้ ซึ่งก็คือดึงมาจากโมเดลเดิม ที่ตัวแปรทำนายทั้ง 9 ตัวอยู่ในโมเดล ค่าที่ 3 `direction` คือวิธีการคัดเลือกโมเดล ในที่นี้เลือกวิธี `stepwise` ส่วนค่าที่ 4 `R2scope` เป็นการสร้างเงื่อนไขว่า ค่า Adjusted $R^2$ ของโมเดลสุดท้าย ต้องไม่มากกว่าโมเดลสูงสุดที่เป็นไปได้ หากการใส่ตัวแปรทำให้ Adjusted $R^2$ สูงกว่าโมเดลสุดท้าย ให้ตัดจบ และค่าที่ 5 `trace` คือการแสดงลำดับการคัดเลือก หากผู้ใช้ต้องการลดการแสดงผล ให้เปลี่ยนตัวนี้เป็น `FALSE` ได้

```{r stepwise}
minmodel <- rda(outcomes ~ 1, data = predictors, scale=TRUE)
stepwise <- ordiR2step(minmodel, scope = formula(out), direction = "stepwise", R2scope = TRUE, trace = TRUE) 
stepwise
```

## การคำนวณอัตราการเจอคนที่ใช่ (Hit Rate) ด้วยแนวคิดของ Taylor & Russell (1939)

ส่วนนี้จะแสดง code ที่ใช้ในการคำนวณผลที่อยู่ในตารางของ Taylor & Russell (1939) โดยท่านต้องกำหนดค่า 3 ค่า ดังนี้

- `sr` คือ อัตราส่วนของคนที่จะเลือกต่อคนที่รับสมัคร (Selection Ratio) กล่าวคือ สัดส่วนของตำแหน่งที่ว่างต่อจำนวนคนสมัครทั้งหมด เช่น การคัดเลือกข้าราชการ มีคนมาสมัคร 10,000 คน ต้องการ 200 คน แสดงว่า `sr=0.02` โดย `sr` มีค่าได้ตั้งแต่ 0 ถึง 1
- `br` คือ อัตราส่วนของคนที่ใช่ (Base Rate) ในประชากรทั้งหมด เช่น จากคนไทยทั้งหมด มีโอกาสเจอคนที่เข้าเกณฑ์ 10% หมายความว่าค่า `br = 0.1` โดย `br` มีค่าเป็นไปได้ตั้งแต่ 0 ถึง 1
- `r` คือ ค่าสหสัมพันธ์ (Correlation) ระหว่างคะแนนที่ใช้คัดเลือก และเกณฑ์ที่ใช้ในการตัดสินว่าคนไหนใช่หรือไม่ เช่น บริษัทใช้แบบทดสอบเชาวน์ปัญญาในการคัดเลือกพนักงาน โดยมีความสัมพันธ์กับผลการปฏิบัติงานเท่ากับ .35 แสดงว่า `r=0.35` โดยแม้ `r` จะมีค่าติดลบได้ ให้เอาเครื่องหมายลบออก แล้วเหมาเป็นค่าบวกให้หมด ค่าที่เป็นไปได้คือตั้งแต่ 0 ถึง 1

**ตัวอย่างที่ 1** อัตราส่วนของคนที่จะเลือกต่อคนที่รับสมัครเท่ากับ 20% อัตราส่วนในการเจอคนที่ใช้เท่ากับ 30% และค่าสหสัมพันธ์ของแบบทดสอบเท่ากับ .25 สามารถแทนค่าใน code ได้ดังข้างล่าง พบว่าโอกาสในการเจอคนที่ใช่เท่ากับ 43%

```{r taylor1}
sr <- 0.2
br <- 0.3
r <- 0.25
library(mvtnorm)
zx <- qnorm(1 - sr)
zy <- qnorm(1 - br)
r <- matrix(c(1, r, r, 1), 2, 2)
q1 <- pmvnorm(lower = c(zx, zy), mean=c(0,0), corr=r)
q1q4 <- pmvnorm(lower = c(zx, -Inf), mean=c(0,0), corr=r)
q1/q1q4
```

**ตัวอย่างที่ 2** อัตราส่วนของคนที่จะเลือกต่อคนที่รับสมัครเท่ากับ 60% อัตราส่วนในการเจอคนที่ใช้เท่ากับ 70% และค่าสหสัมพันธ์ของแบบทดสอบเท่ากับ .5 สามารถแทนค่าใน code ได้ดังข้างล่าง พบว่าโอกาสในการเจอคนที่ใช่เท่ากับ 82%

```{r taylor2}
sr <- 0.6
br <- 0.7
r <- 0.5
library(mvtnorm)
zx <- qnorm(1 - sr)
zy <- qnorm(1 - br)
r <- matrix(c(1, r, r, 1), 2, 2)
q1 <- pmvnorm(lower = c(zx, zy), mean=c(0,0), corr=r)
q1q4 <- pmvnorm(lower = c(zx, -Inf), mean=c(0,0), corr=r)
q1/q1q4
```

## การคำนวณค่าเฉลี่ยของเกณฑ์เมื่อใช้แบบวัด ด้วยแนวคิดของ Naylor & Shines (1965)

ส่วนนี้จะแสดง code ที่ใช้ในการคำนวณผลที่อยู่ในตารางของ Naylor & Shines (1965) โดยท่านต้องกำหนดค่า 2 ค่า คือ `sr` และ `r` ซึ่งมีความหมายเหมือนกับการกำหนดตาราง Taylor & Russell

**ตัวอย่างที่ 1** อัตราส่วนของคนที่จะเลือกต่อคนที่รับสมัครเท่ากับ 20% และค่าสหสัมพันธ์ของแบบทดสอบเท่ากับ .25 สามารถแทนค่าใน code ได้ดังข้างล่าง พบว่าค่าเฉลี่ยของคนที่ถูกคัดเลือกมา จะมีค่าเท่ากับคะแนนมาตรฐาน 0.35

```{r naylor1}
sr <- 0.2
r <- 0.25
zx <- qnorm(1 - sr)
lambda <- dnorm(zx)
r * lambda / sr
```

**ตัวอย่างที่ 2** อัตราส่วนของคนที่จะเลือกต่อคนที่รับสมัครเท่ากับ 60% และค่าสหสัมพันธ์ของแบบทดสอบเท่ากับ .5 สามารถแทนค่าใน code ได้ดังข้างล่าง พบว่าค่าเฉลี่ยของคนที่ถูกคัดเลือกมา จะมีค่าเท่ากับคะแนนมาตรฐาน 0.32

```{r naylor2}
sr <- 0.6
r <- 0.5
zx <- qnorm(1 - sr)
lambda <- dnorm(zx)
r * lambda / sr
```