Linear Regression_bw_LME_tutorial 1 & 2 실습

Data 1: sex, pitch, my.df

1. Create dummy data

pitch <- c(233,204,242,130,112,142) #categorial var.
sex <- c(rep("female", 3), rep("male",3))
my.df <- data.frame(sex, pitch); my.df

##      sex pitch
## 1 female   233
## 2 female   204
## 3 female   242
## 4   male   130
## 5   male   112
## 6   male   142

2. Change vector name

data.frame(gender = sex, pitch)

##   gender pitch
## 1 female   233
## 2 female   204
## 3 female   242
## 4   male   130
## 5   male   112
## 6   male   142

• sex를 gender로 이름 변경

3. Make linear model

xmdl <- lm(pitch ~ sex, my.df)
# 모델
summary(xmdl)

## 
## Call:
## lm(formula = pitch ~ sex, data = my.df)
## 
## Residuals:
##       1       2       3       4       5       6 
##   6.667 -22.333  15.667   2.000 -16.000  14.000 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   226.33      10.18  22.224 2.43e-05 ***
## sexmale       -98.33      14.40  -6.827  0.00241 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17.64 on 4 degrees of freedom
## Multiple R-squared:  0.921,  Adjusted R-squared:  0.9012 
## F-statistic: 46.61 on 1 and 4 DF,  p-value: 0.002407

• my.df에서 sex에 따른 pitch 변화를 보는 linear model
• dependent var: pitch, fixed effect(predictor): sex
• p-value가 유의미함.
• fixed effect가 1개이므로,
  overall model의 p-vlaue(0.002407)와 coefficient p-value(0.00241)가 차이가 거의 없음.
• coeffiecient table은 (만약 fixed effect가 여러개라면) 각 fixed effect 하나씩 봄.
• Intercept(상수항) 값: 226.33
  
• sexmale 값: -98.33
  y = -98.3x +b

4. Calculate mean in terms of sex (sex지정하여)

with(my.df, mean(pitch[sex=='female']))

## [1] 226.3333

mean(my.df[my.df$sex=="female",]$pitch)

## [1] 226.3333

• 두 가지 방법이 있음.
• my.df에서 sex가 female인 사람의 pitch 평균

5. dplyr 활용하여 sex별 pitch 평균 한꺼번에 보기

library(dplyr)

## Warning: 패키지 'dplyr'는 R 버전 4.2.2에서 작성되었습니다

## 
## 다음의 패키지를 부착합니다: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

tapply(my.df$pitch, sex, mean)

##   female     male 
## 226.3333 128.0000

my.df %>% 
  group_by(sex) %>% 
  summarise(mean=mean(pitch))

## # A tibble: 2 × 2
##   sex     mean
##   <chr>  <dbl>
## 1 female  226.
## 2 male    128

• tapply활용하거나
• pipe이용

Data 2: age, pitch2, my.df2

1. age에 따른 pitch 변화보기_make data

age <- c(14,23,35,48,52,67) #continuous var.
pitch2 <- c(252,244,240,233,212,204)
my.df2 <- data.frame(age, pitch2,stringAsFactor = T)
xmdl2 <- lm(pitch2 ~ age, my.df2)
summary(xmdl2)

## 
## Call:
## lm(formula = pitch2 ~ age, data = my.df2)
## 
## Residuals:
##      1      2      3      4      5      6 
## -2.338 -2.149  4.769  9.597 -7.763 -2.115 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 267.0765     6.8522   38.98 2.59e-06 ***
## age          -0.9099     0.1569   -5.80  0.00439 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.886 on 4 degrees of freedom
## Multiple R-squared:  0.8937, Adjusted R-squared:  0.8672 
## F-statistic: 33.64 on 1 and 4 DF,  p-value: 0.004395

• my.df2에서 age에 따른 pitch2 변화를 보는 linear model
• dependent var: pitch, fixed effect(predictor): age
• Intercept 값: 267.0765
  시작점(0살ㅇ)이 높은 이유-> predictor가 age여서 0(살)에 맞추었을 때 값
• age 값: -0.9099
  한살씩 높아지면 intercept 값에서 이만큼 낮아짐. y = -0.9099x +b

#plot
plot(age, pitch2, data= my.df2)

## Warning in plot.window(...): "data"는 그래픽 매개변수가 아닙니다

## Warning in plot.xy(xy, type, ...): "data"는 그래픽 매개변수가 아닙니다

## Warning in axis(side = side, at = at, labels = labels, ...): "data"는 그래픽 매
## 개변수가 아닙니다

## Warning in axis(side = side, at = at, labels = labels, ...): "data"는 그래픽 매
## 개변수가 아닙니다

## Warning in box(...): "data"는 그래픽 매개변수가 아닙니다

## Warning in title(...): "data"는 그래픽 매개변수가 아닙니다

• 1개 이상의 독립변수(predictor)를 사용하면 multiple linear regression(=다중선형회귀)이 됨
  
  • 코드: pitch ~ sex + age + e

linear model이 되기 위한 조건

1. Linearity

• 잔차 그래프(residual plot)가 선형이 되어야 함.

#plot of residuals(rotated)
#잔차(residuals): 표본집단에서의 (기댓값-관측값) 차이=fitted values

plot(fitted(xmdl2), residuals(xmdl2)) #  (rotated!) residual plot: 잔차(fitted vlaues)가 y=0인 x절편위에 있음.
abline(a=0,b=0, col=5, lty="dashed") #좌표에 직선을 그음 0콤마0

2. Absence of collinearity=공산성이 없어야 함

• About 다중 선형 회귀
• 회귀 모형에 포함될 설명변수들 사이에 서로 밀접한 상관관계가 있어서 다중 선형 회귀 모델에서 이들 각각의 개별 효과를 파악하기 힘들게 되는 현상

ex) "how average talking speed affects intelligence ratings?"에 관해 알고자 한다.
아래와 같은 회귀식 만들 수 있음.
intelligence rating ~ talking speed

talking speed를 알아보기 위해 초당 발화하는 syllables/words/sentences를 측정
  --> 말을 빨리하면 초당 발화하는 syllables/words/sentences이 많아지는 것이 당연하다.
    --> 그래서 초당 syllables/words/sentences를 predictors로 사용하면 collinearity 문제가 생김!

• 위의 예시처럼 공산성에 문제가 있으면 모델이 불안정해짐.
• 즉, 만약 multiple predictors가 서로 너무 비슷하다면 어떤 요소가 significant 한지 알지 어려워짐.
• 공산성이 없도록 해야 한다.

3. Homoskedasticity(등분산성) OR ‘absence of heteroskedasticity(이분산성)’

• 분석하는 집단들의 분산이 같아야 함.
• ideal shape: 동그라미 모양
• 좋은 잔차 그래프 예시

plot(rnorm(100),rnorm(100))

4. Normality of residuals

• 잔차는 정규성을 보여야 함.

hist(residuals(xmdl2)) #히스토그램으로 확인 #ideal shape: bell 모양

qqnorm(residuals(xmdl2)) #qqplot으로 확인 #ideal shape: 정규분포처럼 1자

5. Absence of influential data points

• 영향줄만한 데이터가 없어야 한다.
• 만약 특정 데이터가 배제되었다면, coefficients(계수값)가 조정되어야 한다.

dfbeta(xmdl) #DFbeta 값

##     (Intercept)   sexmale
## 1  3.333333e+00 -3.333333
## 2 -1.116667e+01 11.166667
## 3  7.833333e+00 -7.833333
## 4 -1.359740e-16  1.000000
## 5  1.087792e-15 -8.000000
## 6 -9.518180e-16  7.000000

6. Indepencence !

• 통계 테스트에서 가장 중요한 요소..
• 동전의 앞/뒷면 처럼 공존할 수 없음
• 한 참여자로부터 여러개의 응답을 받았다면, 해당 응답들은 서로 독립적이라고 할 수 없음.
• 각각의 참여자는 하나의 data point에 기여해야 한다.

Data 3: Politness(subject, gender, scenario, attitude, freq)

1. Import Politeness Data

politness <- read.csv("C:\\Users\\csjja\\Desktop\\ChosunSL-main\\politeness_data.csv", header=T)
summary(politness)

##    subject             gender             scenario   attitude        
##  Length:84          Length:84          Min.   :1   Length:84         
##  Class :character   Class :character   1st Qu.:2   Class :character  
##  Mode  :character   Mode  :character   Median :4   Mode  :character  
##                                        Mean   :4                     
##                                        3rd Qu.:6                     
##                                        Max.   :7                     
##                                                                      
##    frequency    
##  Min.   : 82.2  
##  1st Qu.:131.6  
##  Median :203.9  
##  Mean   :193.6  
##  3rd Qu.:248.6  
##  Max.   :306.8  
##  NA's   :1

2. scenario를 factor로 다시 지정

• 번외로 해본 practice

politness$scenario<- as.factor(politness$scenario)
class(politness$scenario)

## [1] "factor"

3. missing value 확인

which(is.na(politness$frequency)) #방법1

## [1] 39

which(!complete.cases(politness)) #방법2

## [1] 39

4. 공손성(태도)와 피치 간의 관계를 평균을 통해 boxplot으로 보기

boxplot(frequency ~ attitude*gender, col=c("white","lightgray"),politness)

5. 모델 만들기(ver1)

• dependent var: freq, predictor: attitude, random effect: 참여자와 7가지 시나리오

library(lme4)

## Warning: 패키지 'lme4'는 R 버전 4.2.2에서 작성되었습니다

## 필요한 패키지를 로딩중입니다: Matrix

politeness.model <- lmer(frequency ~ attitude + (1|subject) + (1|scenario), data=politness)
summary(politeness.model)

## Linear mixed model fit by REML ['lmerMod']
## Formula: frequency ~ attitude + (1 | subject) + (1 | scenario)
##    Data: politness
## 
## REML criterion at convergence: 793.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2006 -0.5817 -0.0639  0.5625  3.4385 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  scenario (Intercept)  219     14.80   
##  subject  (Intercept) 4015     63.36   
##  Residual              646     25.42   
## Number of obs: 83, groups:  scenario, 7; subject, 6
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  202.588     26.754   7.572
## attitudepol  -19.695      5.585  -3.527
## 
## Correlation of Fixed Effects:
##             (Intr)
## attitudepol -0.103

• subject의 variation: 219
• intercept: attitude가 informal
• intercept값: 202.588

6. 모델 만들기(ver2)

• gender을 fixed effect로 추가

politeness.model = lmer(frequency ~ attitude + gender + (1|subject) + (1|scenario), data=politness)
summary(politeness.model)

## Linear mixed model fit by REML ['lmerMod']
## Formula: frequency ~ attitude + gender + (1 | subject) + (1 | scenario)
##    Data: politness
## 
## REML criterion at convergence: 775.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2591 -0.6236 -0.0772  0.5388  3.4795 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  scenario (Intercept) 219.5    14.81   
##  subject  (Intercept) 615.6    24.81   
##  Residual             645.9    25.41   
## Number of obs: 83, groups:  scenario, 7; subject, 6
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  256.846     16.116  15.938
## attitudepol  -19.721      5.584  -3.532
## genderM     -108.516     21.013  -5.164
## 
## Correlation of Fixed Effects:
##             (Intr) atttdp
## attitudepol -0.173       
## genderM     -0.652  0.004

• subject의 variation: 205.2(감소).
• intercept:attitude가 informal 이고 gender을 Female
• intercept값: 256.847(증가)

<추가 정보>

• lmetest 사용하면 p-vlaue볼 수 있음.
• general summary statistics
• AIC(Akaike’s Information Criterion)이 적을수록 model fit 더 좋음.