{"id":40,"date":"2018-05-18T20:37:19","date_gmt":"2018-05-19T00:37:19","guid":{"rendered":"http:\/\/sites.nd.edu\/miaoyang\/?page_id=40"},"modified":"2018-05-18T20:48:19","modified_gmt":"2018-05-19T00:48:19","slug":"r-program-for-robust-moderation-analysis","status":"publish","type":"page","link":"https:\/\/sites.nd.edu\/miaoyang\/software\/r-program-for-robust-moderation-analysis\/","title":{"rendered":"R program for robust moderation analysis"},"content":{"rendered":"
\n

<\/h1>\n<\/header>\n
\n

This R program is for performing robust moderation\/interaction analysis using a two-level regression\u00a0model as proposed in the paper \u201cRobust methods for moderation analysis with a two-level regression model\u201d.\u00a0Robustness is achieved using Student\u2019s t distribution and Huber-type weights.<\/p>\n

For citation of the paper, please consider<\/p>\n

Yang, M.<\/strong>, & Yuan, K.-H. (2016). Robust methods for moderation analysis with a two-level regression model.\u00a0Multivariate Behavioral Research<\/em>,\u00a051<\/em>(6), 757-771.<\/p>\n

##=======================================##
\n## A simulated data set with 500 cases ##
\n##=======================================##
\nset.seed(1113)
\nDV <- rt(500,5) # the dependent variable
\nIV <- rnorm(500) # the independent variable
\nModerator <- runif(500,1,10) #the moderator
\ndata <- cbind(DV,IV,Moderator)
\nwrite.table(data,\u2019simdata.txt\u2019,row.names=F)<\/p>\n

##==================##
\n## NML estimation ##
\n##==================##
\nrequire(MASS)
\nNML<-function(y,x,z){
\nparam.names <- as.list(match.call()[-1])<\/p>\n

n <- length(y)
\nh <- cbind(x^2,2*x,rep(1,n))
\nC <- cbind(rep(1,n),z,x,x*z)<\/p>\n

LS.MMR<-lm(y~x*z)
\nr00<-LS.MMR$coefficients[1]
\nr01<-LS.MMR$coefficients[3]
\nr10<-LS.MMR$coefficients[2]
\nr11<-LS.MMR$coefficients[4]
\ngamma.start <- as.matrix(c(r00,r01,r10,r11))
\nsigma.start <- c(1,0,1)<\/p>\n

for (iter in 1:500) {<\/p>\n

tau2<-h%*%sigma.start<\/p>\n

A<-diag(as.vector(1\/tau2))
\nB<-diag(as.vector(1\/tau2^2))<\/p>\n

gamma<-ginv(t(C)%*%A%*%C)%*%(t(C)%*%A%*%y)
\ndelta2<-(y-C%*%gamma)^2
\nsigma<-ginv(t(h)%*%B%*%h)%*%(t(h)%*%B%*%delta2)<\/p>\n

diff<-sum(abs(gamma-gamma.start))+sum(abs(sigma-sigma.start))<\/p>\n

gamma.start<-gamma
\nsigma.start<-sigma<\/p>\n

if (diff<1.0E-12) {
\nconv<-1
\nbreak
\n}
\n}<\/p>\n

tau2<-h%*%sigma
\ndelta2<-(y-C%*%gamma)^2
\ndelta<-y-C%*%gamma<\/p>\n

#BIC
\nk <-length(gamma)+length(sigma)
\nBIC <- n*log(2*pi) + sum(log(tau2)) + sum((delta2)\/tau2) + k*log(n)<\/p>\n

#Sandwich SE
\nU11<-t(C)%*%diag(as.vector(1\/tau2))%*%C; U12<-t(C)%*%diag(as.vector(delta\/tau2^2))%*%h
\nU21<-t(U12); U22<-t(h)%*%diag(as.vector((2*delta2-tau2)\/(2*tau2^3)))%*%h<\/p>\n

V11<-t(C)%*%diag(as.vector(delta2\/tau2^2))%*%C;
\nV12<-t(C)%*%diag(as.vector((delta^3-delta*tau2)\/(2*tau2^3)))%*%h;
\nV21<-t(V12); V22<-t(h)%*%diag(as.vector((delta2-tau2)^2\/(4*tau2^4)))%*%h<\/p>\n

U<-rbind(cbind(U11,U12),cbind(U21,U22))
\nV<-rbind(cbind(V11,V12),cbind(V21,V22))
\nCov<-ginv(U)%*%V%*%ginv(U)
\nse.theta<-sqrt(diag(Cov))<\/p>\n

theta<-c(gamma,sigma)
\nz<-theta\/se.theta
\np<-2*pnorm(-abs(z))
\nR.square<-(gamma[4])^2*var(z)\/((gamma[4])^2*var(z)+sigma[1])<\/p>\n

out<-NULL<\/p>\n

Parameter <- cbind(theta,se.theta,z,p)[1:4,]
\nrownames(Parameter) <- c(\u201cIntercept\u201d,param.names[3:2],paste(param.names[2],param.names[3],sep=\u201d*\u201d))
\ncolnames(Parameter) <- c(\u201cEstimate\u201d,\u201dSE_sw\u201d,\u201dz\u201d,\u201dPr(>|z|)\u201d)
\nout <- list(Parameter=Parameter,BIC=BIC,R.square=R.square)
\nreturn(out)
\n}<\/p>\n

##============================================##
\n## Robust estimation based on Student\u2019s distribution\u00a0 \u00a0 \u00a0 ##
\n##============================================##<\/p>\n

TML<-function(y,x,z,m=5){
\nparam.names <- as.list(match.call()[-1])<\/p>\n

n <- length(y)
\nh <- cbind(x^2,2*x,rep(1,n))
\nC <- cbind(rep(1,n),z,x,x*z)<\/p>\n

LS.MMR<-lm(y~x*z)
\nr00<-LS.MMR$coefficients[1]
\nr01<-LS.MMR$coefficients[3]
\nr10<-LS.MMR$coefficients[2]
\nr11<-LS.MMR$coefficients[4]
\ngamma.start <- as.matrix(c(r00,r01,r10,r11))
\nsigma.start <- c(1,0,1)<\/p>\n

for (iter in 1:500) {<\/p>\n

delta<-y-C%*%gamma.start
\ntau2<-h%*%sigma.start*(m-2)\/m<\/p>\n

di2<-delta^2\/tau2
\nwi<-(m+1)\/(m+di2)<\/p>\n

d<-as.matrix(di2)
\nA<-diag(as.vector(wi\/tau2))
\nB<-diag(as.vector(1\/tau2^2))<\/p>\n

gamma<-ginv(t(C)%*%A%*%C)%*%(t(C)%*%A%*%y)
\nsigma<-ginv(t(h)%*%B%*%h)%*%(t(h)%*%A%*%d)*m\/(m-2)<\/p>\n

diff<-sum(abs(gamma-gamma.start))+sum(abs(sigma-sigma.start))<\/p>\n

gamma.start<-gamma
\nsigma.start<-sigma<\/p>\n

if (diff<1.0E-12) {
\nconv<-1
\nbreak
\n}<\/p>\n

}<\/p>\n

delta<-y-C%*%gamma.start;
\ndelta2<-(y-C%*%gamma.start)^2;
\ntau2<-h%*%sigma.start*(m-2)\/m;
\nd<-as.matrix(delta2\/tau2)
\ndi2<-delta^2\/tau2
\nwi<-(m+1)\/(m+di2)<\/p>\n

#BIC
\nk <-length(gamma)+length(sigma)
\nBIC <- 2*n*log(sqrt(m*pi)*gamma(m\/2)\/gamma(m\/2+0.5)) + sum(log(tau2)) + (m+1)*sum(log(1+d\/m)) + k*log(n)<\/p>\n

#Sandwich SE
\nU11.part<-diag(as.vector((m-d)*wi\/((d+m)*tau2))) #d is di^2
\nU12.part<-diag(as.vector((m-2)*wi*delta\/((d+m)*tau2^2)))
\nU21.part<-U12.part
\nU22.part<-diag(as.vector(0.5*(1-2\/m)^2*((d+2*m)\/(d+m)*wi*d-1)\/tau2^2))
\nV11.part<-diag(as.vector(wi^2*d\/tau2))
\nV12.part<-diag(as.vector((m-2)*(wi*d-1)*wi*delta\/(2*m*tau2^2)))
\nV21.part<-V12.part
\nV22.part<-diag(as.vector((m-2)^2*(wi*d-1)^2\/(4*m^2*tau2^2)))
\nU11<-t(C)%*%U11.part%*%C; U12<-t(C)%*%U12.part%*%h; U21<-t(U12); U22<-t(h)%*%U22.part%*%h
\nV11<-t(C)%*%V11.part%*%C; V12<-t(C)%*%V12.part%*%h; V21<-t(V12); V22<-t(h)%*%V22.part%*%h<\/p>\n

U<-rbind(cbind(U11,U12),cbind(U21,U22))
\nV<-rbind(cbind(V11,V12),cbind(V21,V22))
\nCov<-ginv(U)%*%V%*%ginv(U)
\nse.theta<-sqrt(diag(Cov))<\/p>\n

theta<-c(gamma,sigma)
\nz<-theta\/se.theta
\np<-2*pnorm(-abs(z))
\nR.square<-(gamma[4])^2*var(z)\/((gamma[4])^2*var(z)+sigma[1])<\/p>\n

out<-NULL<\/p>\n

Parameter <- cbind(theta,se.theta,z,p)[1:4,]
\nrownames(Parameter) <- c(\u201cIntercept\u201d,param.names[3:2],paste(param.names[2],param.names[3],sep=\u201d*\u201d))
\ncolnames(Parameter) <- c(\u201cEstimate\u201d,\u201dSE_sw\u201d,\u201dz\u201d,\u201dPr(>|z|)\u201d)
\nout <- list(Parameter=Parameter,BIC=BIC,R.square=R.square)
\nreturn(out)
\n}<\/p>\n

##============================================##
\n## Robust estimation based on Huber-type weights\u00a0 \u00a0 \u00a0 \u00a0 \u00a0##
\n##============================================##<\/p>\n

Huber<-function(y,x,z,alpha=0.95){<\/p>\n

param.names <- as.list(match.call()[-1])<\/p>\n

n <- length(y)
\nh <- cbind(x^2,2*x,rep(1,n))
\nC <- cbind(rep(1,n),z,x,x*z)<\/p>\n

LS.MMR<-lm(y~x*z)
\nr00<-LS.MMR$coefficients[1]
\nr01<-LS.MMR$coefficients[3]
\nr10<-LS.MMR$coefficients[2]
\nr11<-LS.MMR$coefficients[4]
\ngamma.start <- as.matrix(c(r00,r01,r10,r11))
\nsigma.start <- c(1,0,1)<\/p>\n

H2<-qchisq(1-alpha,1)
\nH<-sqrt(H2)
\nkappa=pchisq(H2,3)+H2*alpha
\nq=(H^2\/kappa-1)\/2<\/p>\n

for (iter in 1:500) {<\/p>\n

delta<-y-C%*%gamma.start
\ntau2<-h%*%sigma.start
\ntau2[tau2<=0]<-0.0001<\/p>\n

di<-delta\/sqrt(tau2)
\nwi<- ifelse(abs(di)<H,1,H\/abs(di))
\nwi2 <- (wi^2)\/kappa<\/p>\n

d<-as.matrix(di^2)
\nA<-diag(as.vector(wi\/tau2))
\nA2 <- diag(as.vector(wi2\/tau2))
\nB<-diag(as.vector(1\/tau2^2))<\/p>\n

gamma<-ginv(t(C)%*%A%*%C)%*%(t(C)%*%A%*%y)
\nsigma<-ginv(t(h)%*%B%*%h)%*%(t(h)%*%A2%*%d)<\/p>\n

diff<-sum(abs(gamma-gamma.start))+sum(abs(sigma-sigma.start))<\/p>\n

gamma.start<-gamma
\nsigma.start<-sigma<\/p>\n

if (diff<1.0E-12) {
\nconv<-1
\nbreak
\n}<\/p>\n

}<\/p>\n

tau2<-h%*%sigma; delta2<-(y-C%*%gamma)^2; delta<-y-C%*%gamma; di<-delta\/sqrt(tau2)
\nwi<- ifelse(abs(di)<H,1,H\/abs(di))
\nwi2 <- (wi^2)\/kappa<\/p>\n

tau<-sqrt(tau2)
\ndi<-delta\/tau<\/p>\n

#Sandwich SE
\nU11.part<-ifelse(abs(di)<=H, 1\/tau2, 0)
\nU12.part<-ifelse(abs(di)<=H, delta\/tau2^2, H*sign(di)\/(2*tau^3))
\nU21.part<-ifelse(abs(di)<=H, delta\/kappa\/tau2^2, 0)
\nU22.part<-ifelse(abs(di)<=H, delta2\/kappa\/tau2^3-1\/(2*tau2^2), q\/(tau2^2))<\/p>\n

V11.part<-ifelse(abs(di)<=H, delta2\/tau2^2,H^2\/tau2)
\nV12.part<-ifelse(abs(di)<=H, (delta^3\/kappa-delta*tau2)\/(2*tau2^3),q*H*sign(di)\/(tau^3))
\nV22.part<-ifelse(abs(di)<=H, (delta2\/kappa-tau2)^2\/(4*tau2^4),(q^2)\/(tau2^2))<\/p>\n

U11<-t(C)%*%diag(as.vector(U11.part))%*%C;
\nU12<-t(C)%*%diag(as.vector(U12.part))%*%h
\nU21<-t(h)%*%diag(as.vector(U21.part))%*%C;
\nU22<-t(h)%*%diag(as.vector(U22.part))%*%h<\/p>\n

V11<-t(C)%*%diag(as.vector(V11.part))%*%C; V12<-t(C)%*%diag(as.vector(V12.part))%*%h
\nV21<-t(V12); V22<-t(h)%*%diag(as.vector(V22.part))%*%h<\/p>\n

U<-rbind(cbind(U11,U12),cbind(U21,U22))
\nV<-rbind(cbind(V11,V12),cbind(V21,V22))
\nCov<-ginv(U)%*%V%*%ginv(t(U))
\nse.theta<-sqrt(diag(Cov))<\/p>\n

theta<-c(gamma,sigma)
\nz<-theta\/se.theta
\np<-2*pnorm(-abs(z))
\nR.square<-(gamma[4])^2*var(z)\/((gamma[4])^2*var(z)+sigma[1])<\/p>\n

out<-NULL<\/p>\n

Parameter <- cbind(theta,se.theta,z,p)[1:4,]
\nrownames(Parameter) <- c(\u201cIntercept\u201d,param.names[3:2],paste(param.names[2],param.names[3],sep=\u201d*\u201d))
\ncolnames(Parameter) <- c(\u201cEstimate\u201d,\u201dSE_sw\u201d,\u201dz\u201d,\u201dPr(>|z|)\u201d)
\nout <- list(Parameter=Parameter,BIC=NA,R.square=R.square)
\nreturn(out)
\n}<\/p>\n

##========================================##
\n## Code that needs user\u2019s modification ##
\n##========================================##
\nsetwd(\u201cc:\/moderation\u201d) #set working directory
\ndata <- read.table(\u2018simdata.txt\u2019,header=T) #read data
\ny <- data$DV
\nx <- data$IV
\nz <- data$Moderator
\nNML(y,x,z)
\nTML(y,x,z,m=5)
\nHuber(y,x,z,alpha=.05)<\/p>\n<\/div>\n

 <\/p>\n

 <\/p>\n

This R program is for performing robust moderation\/interaction analysis using a two-level regression\u00a0model as proposed in the paper \u201cRobust methods for moderation analysis with a two-level regression model\u201d.\u00a0Robustness is achieved using Student\u2019s t distribution and Huber-type weights. For citation of the paper, please consider Yang, M., & Yuan, K.-H. (2016). Robust methods for moderation analysis […]<\/p>\n","protected":false},"author":2828,"featured_media":0,"parent":16,"menu_order":0,"comment_status":"closed","ping_status":"closed","template":"","meta":{"footnotes":""},"class_list":["post-40","page","type-page","status-publish","hentry"],"_links":{"self":[{"href":"https:\/\/sites.nd.edu\/miaoyang\/wp-json\/wp\/v2\/pages\/40","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/sites.nd.edu\/miaoyang\/wp-json\/wp\/v2\/pages"}],"about":[{"href":"https:\/\/sites.nd.edu\/miaoyang\/wp-json\/wp\/v2\/types\/page"}],"author":[{"embeddable":true,"href":"https:\/\/sites.nd.edu\/miaoyang\/wp-json\/wp\/v2\/users\/2828"}],"replies":[{"embeddable":true,"href":"https:\/\/sites.nd.edu\/miaoyang\/wp-json\/wp\/v2\/comments?post=40"}],"version-history":[{"count":5,"href":"https:\/\/sites.nd.edu\/miaoyang\/wp-json\/wp\/v2\/pages\/40\/revisions"}],"predecessor-version":[{"id":48,"href":"https:\/\/sites.nd.edu\/miaoyang\/wp-json\/wp\/v2\/pages\/40\/revisions\/48"}],"up":[{"embeddable":true,"href":"https:\/\/sites.nd.edu\/miaoyang\/wp-json\/wp\/v2\/pages\/16"}],"wp:attachment":[{"href":"https:\/\/sites.nd.edu\/miaoyang\/wp-json\/wp\/v2\/media?parent=40"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}