헤이즈 모델 다 만들어 보고 싶다. 1번 해보장
#데이터 만들기
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.9756724##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
}
boot1(1,2,3,d,1000)## $moderation_mean_BootSE
## [1] 0.99010252 0.01310065
##
## $moderation_CI
## 0.1% 1% 5% 10% 90% 95% 99% 99.9%
## 0.9544036 0.9606966 0.9685565 0.9726658 1.0066389 1.0120625 1.0207615 1.0248795##다른표기 생각해보니 +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
}
boot1r(1,2,3,d,10)## $moderation_mean_BootSE
## [1] 0.99529671 0.01300555
##
## $moderation_CI
## 0.1% 1% 5% 10% 90% 95% 99% 99.9%
## 0.9765921 0.9771685 0.9797304 0.9829328 1.0128686 1.0160875 1.0186625 1.0192419
##
## $`+1SD_mean_BootSE`
## [1] 4.6345749 0.1970077
##
## $`+1SD_CI`
## 0.1% 1% 5% 10% 90% 95% 99% 99.9%
## 4.316725 4.327053 4.372954 4.430331 4.867672 4.917263 4.956936 4.965862
##
## $`-1SD_mean_BootSE`
## [1] -1.085306 0.200605
##
## $`-1SD_CI`
## 0.1% 1% 5% 10% 90% 95% 99%
## -1.3799588 -1.3694278 -1.3226233 -1.2641177 -0.8033253 -0.7619762 -0.7288970
## 99.9%
## -0.7214541