Hayes Model 1

 이전에 Rpub에 정리해 둔 것을 여기 옮겨둔다. 티스토리는 어떻게 사용하는지 모르는 컴맹이라 RMD파일에 있던 것을 복붙 한다. 

---

 사실 나는 Hayes의 부트스트래핑을 활용한 가설검정을 좋아하지 않는다. 애초에 표본을 복원추출해 정규성을 확보하고 CI를 통해 가설검정을 한다는 것 자체가 타당한 방법이라 생각하지 않기 때문이다. 표본분포의 정규성은 이미 표집 된 이상 연구자가 건드리는 것으로 해결 할 문제가 아니며 GLM자체가 정규성을 가정하고 검정해야 한다고 본다. 

 그래도 내 학위논문도 이걸 썼다.... 그때는 패키지를 썼는데 나름대로 함수를 만들어 보았다.

---

#데이터 만들기

library(stats)
#model 1
x<-rnorm(1000)
w<-rnorm(1000)+ x^2
y<-rnorm(1000, 0,1)  + w*x
co1<-rnorm(1000)

d<-data.frame(x,w,y,co1)

#함수1 조절 변수의 CI만 나오게 해보았다. 일단 estimates만 뽑아내는 함수 만들었다. 함수 입력값은 각 변수의 열 숫자를 쓰게 하였다. 복원추출로 하는거라 할때마다 살짝씩 다를 수 있다(당연하다).

boot1_1<-function(xxx,mmm,yyy,d){
    n<-sample(1:nrow(d),nrow(d),replace = T)
    nnk<-d[n,]
    nnk<-as.data.frame(nnk)
    k2<-lm(nnk[,yyy]~ nnk[,xxx]+nnk[,mmm] + nnk[,xxx]*nnk[,mmm], data=nnk)
    s2<-summary(k2)
    coem<-s2$coefficients
    eff<-as.data.frame(coem)
    eff<-eff[nrow(eff),1]
    eff}
boot1_1(1,2,3,d)

이럴경우 조절변수의 회귀계수는 다음과 같이 나온다.

[1] 0.9867756

##boot 반복구문으로 CI 를 구해보자. bootnum은 몇번 반복할지이고 d는 데이터이다. 

boot1<-function(xxx,mmm,yyy,d,bootnum){
  ###estimate a*m
  boot1_1<-function(xxx,mmm,yyy,d){
    n<-sample(1:nrow(d),nrow(d),replace = T)
    nnk<-d[n,]
    nnk<-as.data.frame(nnk)
    k2<-lm(nnk[,yyy]~ nnk[,xxx]+nnk[,mmm] + nnk[,xxx]*nnk[,mmm], data=nnk)
    s2<-summary(k2)
    coem<-s2$coefficients
    eff<-as.data.frame(coem)
    eff<-eff[nrow(eff),1]
    eff}
  k<-1
  l<-matrix(rep(NA,bootnum),ncol = 1)
  l<-as.data.frame(l)
  repeat{
    l[k,]<-boot1_1(xxx,mmm,yyy,d)
    k<-k+1
    if(k>=bootnum+1) break
  }
  estimates<-list(l)
  ci1<-quantile(l[,1],probs = c(.001,0.01,0.05,0.10,0.90,0.95,0.99,.999))
  kmkmkmkm<-list(c(mean(l[,1]),sd(l[,1])),ci1)
  names(kmkmkmkm)<-c("moderation_mean_BootSE", "moderation_CI")
  kmkmkmkm
}

대략 이전 회귀계수 산출 함수를 repeat해서 bootnum까지 구한 뒤 이에 대한 CI를 구하는 것이다. 

결과를 구해보면 다음과 같다. 

boot1(1,2,3,d,1000)

#결과 값
$moderation_mean_BootSE
[1] 0.99227633 0.01413421   ##첫번째는 계수의 평균이고 두번째는 SE다.

$moderation_CI
     0.1%        1%        5%       10%       90%       95%       99%     99.9% 
0.9540202 0.9614998 0.9693394 0.9747261 1.0103719 1.0170332 1.0278129 1.0412244

##+1sd, -1sd로 나눠서 x->y 개수 구하는 경향이 있으니 그 값도 넣어보자

##model 1 +sd,-sd
boot1r<-function(xxx,mmm,yyy,d,bootnum){
  ###estimate a*m
  boot1_1<-function(xxx,mmm,yyy,d){
    n<-sample(1:nrow(d),nrow(d),replace = T)
    nnk<-d[n,]
    nnk<-as.data.frame(nnk)
    nnk2<-nnk
    nnk3<-nnk
    k2<-lm(nnk[,yyy]~ nnk[,xxx]+nnk[,mmm] + nnk[,xxx]*nnk[,mmm], data=nnk)
    s2<-summary(k2)
    coem<-s2$coefficients
    eff<-as.data.frame(coem)
    eff<-eff[nrow(eff),1]
    
    nnk2[,mmm]<-ifelse(nnk2[,mmm]> mean(nnk2[,mmm])+sd(nnk2[,mmm]), nnk2[,mmm], NA)
    nnk2<-na.omit(nnk2)
    k3<-lm(nnk2[,yyy] ~ nnk2[,xxx], data = nnk2)
    s3<-summary(k3)
    coem2<-s3$coefficients
    effh<-as.data.frame(coem2)
    effh<-effh[nrow(effh),1]
    
    nnk3[,mmm]<-ifelse(nnk3[,mmm] < (mean(nnk3[,mmm])-sd(nnk3[,mmm])), nnk3[,mmm], NA)
    nnk3<-na.omit(nnk3)
    k4<-lm(nnk3[,yyy] ~ nnk3[,xxx], data = nnk3)
    s4<-summary(k4)
    coem3<-s4$coefficients
    effl<-as.data.frame(coem3)
    effl<-effl[nrow(effl),1]
    
    effff<-c(eff, effh,effl)
    effff<-matrix(effff,ncol=3)
    effff}
  k<-1
  l<-matrix(rep(NA,bootnum*3),ncol = 3)
  l<-as.data.frame(l)
  repeat{
    l[k,]<-boot1_1(xxx,mmm,yyy,d)
    k<-k+1
    if(k>=bootnum+1) break
  }
  estimates<-list(l)
  ci1<-quantile(l[,1],probs = c(.001,0.01,0.05,0.10,0.90,0.95,0.99,.999))
  ci2<-quantile(l[,2],probs = c(.001,0.01,0.05,0.10,0.90,0.95,0.99,.999))
  ci3<-quantile(l[,3],probs = c(.001,0.01,0.05,0.10,0.90,0.95,0.99,.999))
  kmkmkmkm<-list(c(mean(l[,1]),sd(l[,1])),ci1,
                 c(mean(l[,2]),sd(l[,2])),ci2,
                 c(mean(l[,3]),sd(l[,3])),ci3)
  names(kmkmkmkm)<-c("moderation_mean_BootSE", "moderation_CI","+1SD_mean_BootSE", "+1SD_CI","-1SD_mean_BootSE", "-1SD_CI")
  kmkmkmkm
}

다시 결과를 보자

$moderation_mean_BootSE
[1] 1.00558515 0.01318121

$moderation_CI
     0.1%        1%        5%       10%       90%       95%       99%     99.9% 
0.9924395 0.9924940 0.9927362 0.9930390 1.0210891 1.0246562 1.0275098 1.0281519 

$`+1SD_mean_BootSE`
[1] 4.6999511 0.3937798

$`+1SD_CI`
    0.1%       1%       5%      10%      90%      95%      99%    99.9% 
4.021768 4.051060 4.181247 4.343981 4.989223 5.243296 5.446555 5.492289 

$`-1SD_mean_BootSE`
[1] -1.1135182  0.2269987

$`-1SD_CI`
      0.1%         1%         5%        10%        90%        95%        99%      99.9% 
-1.4088407 -1.4009186 -1.3657093 -1.3216977 -0.8810358 -0.7745794 -0.6894142 -0.6702521

사실 여기에 CI를 포함한 회귀선 그래프를 그려줘야 하지만 이건 다음에 도전해보겠다. 

솔직히 그냥 패키지(속도는 느리지만 lavaan 추천)를 쓰거나 spss를 쓰는걸 추천한다. 이 코드를 쓰려면 회귀식을 알아야 하며 공변을 넣으려면 코드를 수정해야 한다....