R Markdown 필요 패키지 설치
library("readxl")
## Warning: 패키지 'readxl'는 R 버전 4.2.1에서 작성되었습니다
library("dplyr")
##
## 다음의 패키지를 부착합니다: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library("lme4")
## Warning: 패키지 'lme4'는 R 버전 4.2.1에서 작성되었습니다
## 필요한 패키지를 로딩중입니다: Matrix
library("ordinal")
## Warning: 패키지 'ordinal'는 R 버전 4.2.1에서 작성되었습니다
##
## 다음의 패키지를 부착합니다: 'ordinal'
## The following objects are masked from 'package:lme4':
##
## ranef, VarCorr
## The following object is masked from 'package:dplyr':
##
## slice
library("lmerTest")
## Warning: 패키지 'lmerTest'는 R 버전 4.2.1에서 작성되었습니다
##
## 다음의 패키지를 부착합니다: 'lmerTest'
## The following object is masked from 'package:lme4':
##
## lmer
## The following object is masked from 'package:stats':
##
## step
library("knitr")
## Warning: 패키지 'knitr'는 R 버전 4.2.1에서 작성되었습니다
폴더위치 지정 및 데이터 불러오기
setwd("C:/Users/qustj/Desktop/Startle_test")
dts<-read_xlsx("Startle test data sheet(statistics)_Dr.Steibel.xlsx",skip = 1)
## New names:
## • `` -> `...11`
## • `` -> `...12`
## • `` -> `...13`
head(dts)
## # A tibble: 6 × 13
## `Replicate#` `Pen#` `Pig#` `SR or SV\r\n(S…` `Relaxed-\r\nt…` `Orientation\r…`
## <dbl> <chr> <dbl> <chr> <dbl> <dbl>
## 1 1 ML5 3 SV 3 1
## 2 1 ML6 2 SV 3 3
## 3 1 ML6 3 SR 2 3
## 4 1 ML6 5 SV 2 3
## 5 1 ML6 6 SR 2 1
## 6 1 MR6 7 SR 2 3
## # … with 7 more variables: `Startle\r\nresponse\r\n(Startleresp)` <chr>,
## # `Outcome behavior\r\nafter end freeze\r\n(OAF)` <chr>,
## # `Startle magnitude\r\nscore\r\n(0-4)\r\n(Startlemag)` <dbl>,
## # `Time taken to resume\r\n(sec)\r\n(ttr)` <dbl>, ...11 <lgl>, ...12 <chr>,
## # ...13 <chr>
dts 데이터에서 “unique_pen이란 새로운 열(pen열을 factor로 변환한
것)을 생성, orientation 열을 factor로 변환
dts<-mutate(dts,unique_pen=as.factor(paste(replicate,pen)),orientation=as.factor(orientation))
head(dts)
## # A tibble: 6 × 11
## replicate pen pig stress RTS orientation startleresp OAF startlemag
## <dbl> <chr> <dbl> <chr> <dbl> <fct> <chr> <chr> <dbl>
## 1 1 ML5 3 SV 3 1 SI LD 1
## 2 1 ML6 2 SV 3 3 ST LD 1
## 3 1 ML6 3 SR 2 3 RA FD 3
## 4 1 ML6 5 SV 2 3 RA FD 3
## 5 1 ML6 6 SR 2 1 RA LD 3
## 6 1 MR6 7 SR 2 3 RA NIC 3
## # … with 2 more variables: ttr <dbl>, unique_pen <fct>
dts 데이터에서 pen 열의 숫자만 딴 pen_pos열을 생성
dts$pen_pos<-as.numeric(substr(dts$pen,3,3))
head(dts)
## # A tibble: 6 × 12
## replicate pen pig stress RTS orientation startleresp OAF startlemag
## <dbl> <chr> <dbl> <chr> <dbl> <fct> <chr> <chr> <dbl>
## 1 1 ML5 3 SV 3 1 SI LD 1
## 2 1 ML6 2 SV 3 3 ST LD 1
## 3 1 ML6 3 SR 2 3 RA FD 3
## 4 1 ML6 5 SV 2 3 RA FD 3
## 5 1 ML6 6 SR 2 1 RA LD 3
## 6 1 MR6 7 SR 2 3 RA NIC 3
## # … with 3 more variables: ttr <dbl>, unique_pen <fct>, pen_pos <dbl>
boxplot으로 데이터분포 확인
boxplot(dts$ttr~dts$stress)
함수사용목적 : 행동까지 걸린시간과 스트레스 사이의 유의성을
보기위해
lmer = fit linear mixed-effects model(선형 혼합 효과모델)
lmer 함수로 unique_pen 열을 랜덤효과로, stress group, orientation,
pen position, RTS(이완/긴장상태)를 고정효과로 지정
lmm<-lmer(ttr~stress+orientation+as.factor(RTS)+pen_pos+(1|unique_pen),data=dts)
summary(lmm)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ttr ~ stress + orientation + as.factor(RTS) + pen_pos + (1 |
## unique_pen)
## Data: dts
##
## REML criterion at convergence: 442.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0787 -0.7016 -0.2009 0.6386 1.9318
##
## Random effects:
## Groups Name Variance Std.Dev.
## unique_pen (Intercept) 77.67 8.813
## Residual 680.24 26.081
## Number of obs: 52, groups: unique_pen, 21
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 85.216 19.727 34.271 4.320 0.000127 ***
## stressSV 5.844 8.189 43.753 0.714 0.479271
## orientation2 27.105 11.031 43.314 2.457 0.018079 *
## orientation3 1.887 9.038 41.776 0.209 0.835631
## as.factor(RTS)2 -21.950 14.012 38.211 -1.566 0.125479
## as.factor(RTS)3 -39.039 22.061 43.599 -1.770 0.083783 .
## pen_pos -4.771 3.361 24.025 -1.420 0.168582
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) strsSV orntt2 orntt3 a.(RTS)2 a.(RTS)3
## stressSV -0.288
## orientatin2 -0.211 0.273
## orientatin3 -0.293 0.159 0.436
## as.fc(RTS)2 -0.376 -0.159 -0.250 -0.057
## as.fc(RTS)3 -0.100 -0.312 -0.155 -0.157 0.620
## pen_pos -0.723 0.170 0.157 0.109 -0.241 -0.285
anova(lmm)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## stress 346.4 346.39 1 43.753 0.5092 0.47927
## orientation 4733.1 2366.56 2 43.135 3.4790 0.03974 *
## as.factor(RTS) 2373.5 1186.75 2 43.168 1.7446 0.18681
## pen_pos 1370.8 1370.78 1 24.025 2.0151 0.16858
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 분석결과 스트레스/탄력 상태 측면(RTS)에서 유의미한 차이를 보이지 않아 새로운 요소 추가필요
함수사용목적 : RTS에서 유의미한 것을 발견하기 위해 놀라는 크기를
점수화한 열과 스트레스 사이의 유의성을 보자
lmer 함수로 효과는 동일하지만 가우스변수로 가정한 새로운 식
생성
lmm은 ttr~stress 사이의 유의성확인 / lmm2는 startlemag~stress 사이의
유의성확인
lmm2<-lmer(startlemag~stress+orientation+as.factor(RTS)+pen_pos+(1|unique_pen),data=dts)
summary(lmm2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: startlemag ~ stress + orientation + as.factor(RTS) + pen_pos +
## (1 | unique_pen)
## Data: dts
##
## REML criterion at convergence: 123.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9310 -0.6533 0.1134 0.5695 1.9619
##
## Random effects:
## Groups Name Variance Std.Dev.
## unique_pen (Intercept) 0.1254 0.3542
## Residual 0.5194 0.7207
## Number of obs: 52, groups: unique_pen, 21
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4.81559 0.58636 31.74231 8.213 2.36e-09 ***
## stressSV -0.07038 0.23278 43.03418 -0.302 0.76385
## orientation2 0.13259 0.31963 44.55328 0.415 0.68027
## orientation3 -0.16748 0.26333 42.97034 -0.636 0.52816
## as.factor(RTS)2 -0.33252 0.41112 40.78797 -0.809 0.42332
## as.factor(RTS)3 -1.70314 0.63910 44.43824 -2.665 0.01069 *
## pen_pos -0.29285 0.10198 22.47256 -2.872 0.00875 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) strsSV orntt2 orntt3 a.(RTS)2 a.(RTS)3
## stressSV -0.259
## orientatin2 -0.193 0.258
## orientatin3 -0.281 0.124 0.446
## as.fc(RTS)2 -0.366 -0.186 -0.260 -0.032
## as.fc(RTS)3 -0.111 -0.331 -0.153 -0.131 0.632
## pen_pos -0.732 0.173 0.151 0.088 -0.242 -0.275
anova(lmm2)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## stress 0.0475 0.0475 1 43.034 0.0914 0.763855
## orientation 0.5264 0.2632 2 44.198 0.5068 0.605895
## as.factor(RTS) 4.3526 2.1763 2 44.111 4.1902 0.021572 *
## pen_pos 4.2828 4.2828 1 22.473 8.2460 0.008747 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(lmm2)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## stress 0.0475 0.0475 1 43.034 0.0914 0.763855
## orientation 0.5264 0.2632 2 44.198 0.5068 0.605895
## as.factor(RTS) 4.3526 2.1763 2 44.111 4.1902 0.021572 *
## pen_pos 4.2828 4.2828 1 22.473 8.2460 0.008747 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# RTS에서 유의미한 것을 확인. 하지만 랜덤효과때문에 이러한 유의를 보이는건지 확인필요
함수사용목적 : RTS에서 유의미한 것이 랜덤효과때문인지 아닌지
확인하기 위해서
clm = cumulative link model(누적 링크모델)
clm 함수사용 전에 factor로 변환해야 함수가 사용되는듯
dts$startlemag<-factor(dts$startlemag,ordered = TRUE)
clm 함수로 랜덤효과를 제거한 새로운식 생성
랜덤효과를 제거해도 RTS에서 유의미한 차이 발견. 하지만 여기에서는
랜덤효과가 있는것이 더 낫다라고 결론지음.
lmm3 = clm(startlemag ~ stress+orientation+as.factor(RTS)+pen_pos,data=dts)
summary(lmm3)
## formula: startlemag ~ stress + orientation + as.factor(RTS) + pen_pos
## data: dts
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 52 -53.38 124.77 5(0) 5.85e-09 3.7e+03
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## stressSV -0.1004 0.6162 -0.163 0.87058
## orientation2 1.0193 0.8195 1.244 0.21358
## orientation3 -0.2185 0.6175 -0.354 0.72340
## as.factor(RTS)2 -1.4023 1.0500 -1.335 0.18173
## as.factor(RTS)3 -4.3629 1.7632 -2.474 0.01334 *
## pen_pos -0.7544 0.2555 -2.953 0.00314 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 1|2 -8.078 1.925 -4.196
## 2|3 -6.281 1.793 -3.503
## 3|4 -4.035 1.647 -2.450
anova(lmm3, type="III")
## Type III Analysis of Deviance Table with Wald chi-square tests
##
## Df Chisq Pr(>Chisq)
## stress 1 0.0265 0.870581
## orientation 2 2.4672 0.291239
## as.factor(RTS) 2 6.1596 0.045968 *
## pen_pos 1 8.7219 0.003144 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 랜덤효과 제거하나 안하나 RTS에서 유의미하게 나왔는데 왜 랜덤효과가 있는것이 더 낫다라는 결론이 나왔는지..?
# 랜덤효과가 있는 ordinal regression 식 발견!
lmer 함수와 clmm2 함수의 차이점만 알면 될듯.. 사용된 변수와
고정,랜덤효과가 모두 동일한데 왜 clmm2로 사용하는지?
lmm4 = clmm2(startlemag ~ stress+orientation+as.factor(RTS)+pen_pos,random=unique_pen,data=dts,Hess=TRUE)
## Warning: glm.fit: 적합된 확률값들이 0 또는 1 입니다
lmm4
## Cumulative Link Mixed Model fitted with the Laplace approximation
##
## Call:
## clmm2(location = startlemag ~ stress + orientation + as.factor(RTS) +
## pen_pos, random = unique_pen, data = dts, Hess = TRUE)
##
## Random effects:
## Var Std.Dev
## unique_pen 0.403829 0.6354754
##
## Location coefficients:
## stressSV orientation2 orientation3 as.factor(RTS)2 as.factor(RTS)3
## -0.0999492 0.7377116 -0.4037286 -1.3280641 -4.5642544
## pen_pos
## -0.8185282
##
## No Scale coefficients
##
## Threshold coefficients:
## 1|2 2|3 3|4
## -8.633008 -6.742234 -4.322134
##
## log-likelihood: -53.17019
## AIC: 126.3404
summary(lmm4)
## Cumulative Link Mixed Model fitted with the Laplace approximation
##
## Call:
## clmm2(location = startlemag ~ stress + orientation + as.factor(RTS) +
## pen_pos, random = unique_pen, data = dts, Hess = TRUE)
##
## Random effects:
## Var Std.Dev
## unique_pen 0.403829 0.6354754
##
## Location coefficients:
## Estimate Std. Error z value Pr(>|z|)
## stressSV -0.0999 0.6472 -0.1544 0.87726369
## orientation2 0.7377 0.9968 0.7401 0.45923977
## orientation3 -0.4037 0.7470 -0.5405 0.58886889
## as.factor(RTS)2 -1.3281 1.1658 -1.1392 0.25461906
## as.factor(RTS)3 -4.5643 1.9385 -2.3546 0.01854400
## pen_pos -0.8185 0.3031 -2.7004 0.00692486
##
## No scale coefficients
##
## Threshold coefficients:
## Estimate Std. Error z value
## 1|2 -8.6330 2.2929 -3.7652
## 2|3 -6.7422 2.1020 -3.2075
## 3|4 -4.3221 1.8720 -2.3088
##
## log-likelihood: -53.17019
## AIC: 126.3404
## Condition number of Hessian: 3780.96
OAF(정지 종료후 하는 동작)와 스트레스 사이의 연관성분석은 카이스퀘어
검정으로
다항응답모델이기 때문에 더욱 까다롭다?
indt<-chisq.test(dts$OAF,dts$stress,simulate.p.value = TRUE)
print(indt)
##
## Pearson's Chi-squared test with simulated p-value (based on 2000
## replicates)
##
## data: dts$OAF and dts$stress
## X-squared = 7.2466, df = NA, p-value = 0.2559
OAF와 스트레스의 관측치와 기댓값 비교
indt$observed
## dts$stress
## dts$OAF SR SV
## CP 1 0
## E 1 0
## EX 7 10
## FD 4 8
## L 1 0
## LD 4 5
## NIC 8 3
indt$expected
## dts$stress
## dts$OAF SR SV
## CP 0.5 0.5
## E 0.5 0.5
## EX 8.5 8.5
## FD 6.0 6.0
## L 0.5 0.5
## LD 4.5 4.5
## NIC 5.5 5.5
ttr에 대한 studentized residuals 값을 이용한 plot 그리기
qqnorm(rstudent(lmm))
abline(0,1,col="red",lwd=2)
