# Power analysis for cross-classified model

#### Junes

##### Member
Hi there,

Any help would be greatly appreciated!

We have done a study in which a sample of respondents rated a sample of online profiles on their trustworthiness. There were a total of 189 respondents and 259 profiles. So, each profile was rated on average 7.3 times. We are using a cross-classified linear regression model to account for respondent and profile variability.

We've done some a priori power simulations, and now we want to expand on that with some more post-hoc simulations. We used a program, MLPowSim that can help us with that. This program generates R code based on our answers, which were hopefully correct (see below)

However, the problem is that we want to report a power based on a standardized unit, e.g. "A power of 0.8 to detect a correlation of 0.1". The numbers we entered are in the units of our measuring scales. We are not sure how to convert these and enter them into the R code, so that we can get a power for a given correlation, say, 0.1. Especially given the complex design of our study.

It's a bit complex to explain, hopefully I've made myself somewhat clear. Please feel free to ask any questions.

Code:
#       A programme to obtain the power of parameters in an
#             XC unbalanced model with  Normal response
#  NB: unbalanced with fixed probability of non-response in first level
#                 generated on 19/11/18
#####~~~~~~~~~~~~~~~~~     Required packages  ~~~~~~~~~~~~~~~~~~~~~~~~#####
library(MASS)
library(lme4)
#####~~~~~~~~~~~~~~~~~~~    Initial inputs    ~~~~~~~~~~~~~~~~~~~~~~~~#####

set.seed(1)
siglevel<-0.050
z1score<-abs(qnorm(siglevel))
simus<-100
n1low<-189
n1high<-189
n1step<-0
n2low<-259
n2high<-259
n2step<-0
problev1out<-0.5
npred<-1
randsize<-4
rand1size<-1
rand2size<-1
repliclow<-1
replichigh<-1
replicstep<-0
beta<-c(4.663,0.011)
betasize<-length(beta)
effectbeta<-abs(beta)
sgnbeta<-sign(beta)
meanpred<-c(0,15.700)
varpred<-c(0,0)
varpred1<-c(0,0)
varpred2<-c(0,82.450)
sigma2u<-matrix(c(0.746),rand1size,rand1size)
sigma2v<-matrix(c(0.198),rand2size,rand2size)
sigmae<-sqrt(0.506)
n1range<-seq(n1low,n1high,n1step)
n2range<-seq(n2low,n2high,n2step)
reprange<-seq(repliclow,replichigh,replicstep)
n1size<-length(n1range)
n2size<-length(n2range)
repsize<-length(reprange)
totalsize<-n1size*n2size*repsize
finaloutput<-matrix(0,totalsize,6*betasize)
unbalprob<- 0.961
rowcount<-1
##~~~~~~~~~~~~~~~~~~          Creating grouped data       ~~~~~~~~~~~~~~~~~~##

fixname<-c("x0","x1")
fixform<-"1+x1"
rand1form<-"(1|l1id)"
rand2form<-"(1|l2id)"
expression<-paste(c(fixform,rand1form,rand2form),collapse="+")
modelformula<-formula(paste("y ~",expression))
data<-vector("list",3+length(fixname))
names(data)<-c("l1id","l2id","y",fixname)

#####--------- Initial input for power in two approaches ----------------#####

powaprox<-vector("list",betasize)
names(powaprox)<-c("b0","b1")
powsde<-powaprox

cat("               The programme was executed at", date(),"\n")
cat("--------------------------------------------------------------------\n")
for(n2 in seq(n2low,n2high,n2step)){
n2run<-n2
for(n1 in seq(n1low,n1high,n1step)){
n1run<-n1
for(replication in seq(repliclow,replichigh,replicstep)){

cat("Start of simulation for ",n2run," XC2,",n1run," XC1 and ",replication," 1-level units\n")
powaprox[1:betasize]<-rep(0,betasize)
powsde<-powaprox

sdepower<-matrix(0,betasize,simus)

###-----------------------     Simulation step    -----------------------###

for(iter in 1:simus){
if (iter/10 == floor(iter/10)){
cat(" Iteration remain=",simus-iter,"\n")
}
xcunbal<-rbinom(n1*n2,replication,1-unbalprob)
l1id<-rep(rep(1:n1,each=n2),xcunbal)
l2id<-rep(rep(1:n2,n1),xcunbal)
cumxc<-c(0,cumsum(xcunbal))
length=cumxc[length(cumxc)]
x<-matrix(1,length,betasize)
z1<-matrix(1,length,rand1size)
z2<-matrix(1,length,rand2size)

##~~~~~~~~~~~~~~~~~~~~~       Set up X matrix       ~~~~~~~~~~~~~~~~~~~~~##

micpred<-rnorm(length,meanpred,sqrt(varpred))
macpred1<-rnorm(n1,0,sqrt(varpred1))
macpred2<-rnorm(n2,0,sqrt(varpred2))
x[,2]<-micpred+macpred1[l1id]+macpred2[l2id]

#######----------------------------------------------------########

u<-mvrnorm(n1,rep(0,rand1size),sigma2u)
v<-mvrnorm(n2run,rep(0,rand2size),sigma2v)
fixpart<-x%*%beta
rand1part<-rowSums(z1*u[l1id,])
rand2part<-rowSums(z2*v[l2id,])
e<-rnorm(length,0,sigmae)
y<-fixpart+rand1part+rand2part+e

##-------------------        Inputs for model fitting       ---------------##

data$l1id<-as.factor(l1id) data$l2id<-as.factor(l2id)
data$y<-y data$x0<-x[,1]
data\$x1<-x[,2]
###~~~~~~~~~~      Fitting the model using lmer funtion    ~~~~~~~~~~###

(fitmodel <- lmer(modelformula,data,method="REML"))

#~~~~~~~~~~         To obtain the power of parameter(s)       ~~~~~~~~~~#

estbeta<-fixef(fitmodel)
sdebeta<-sqrt(diag(vcov(fitmodel)))
for(l in 1:betasize)
{
cibeta<-estbeta[l]-sgnbeta[l]*z1score*sdebeta[l]
if(beta[l]*cibeta>0)              powaprox[[l]]<-powaprox[[l]]+1
sdepower[l,iter]<-as.numeric(sdebeta[l])
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

} ##  iteration end here

###---------                  Powers and their CIs             ---------###

for(l in 1:betasize){

meanaprox<-powaprox[[l]]<-unlist(powaprox[[l]]/simus)
Laprox<-meanaprox-z1score*sqrt(meanaprox*(1-meanaprox)/simus)
Uaprox<-meanaprox+z1score*sqrt(meanaprox*(1-meanaprox)/simus)
meansde<-mean(sdepower[l,])
varsde<-var(sdepower[l,])
USDE<-meansde-z1score*sqrt(varsde/simus)
LSDE<-meansde+z1score*sqrt(varsde/simus)
powLSDE<- pnorm(effectbeta[l]/LSDE-z1score)
powUSDE<- pnorm(effectbeta[l]/USDE-z1score)
powsde[[l]]<-pnorm(effectbeta[l]/meansde-z1score)

###        Restrict the CIs within 0 and 1        ###
if(Laprox<0) Laprox<-0
if(Uaprox>1) Uaprox<-1
if(powLSDE<0) powLSDE<-0
if(powUSDE>1) powUSDE<-1

finaloutput[rowcount,(6*l-5):(6*l-3)]<-c(Laprox,meanaprox,Uaprox)
finaloutput[rowcount,(6*l-2):(6*l)]<-c(powLSDE,powsde[[l]],powUSDE)

}

rowcount<-rowcount+1
cat("--------------------------------------------------------------------\n")
} ## end of the loop over the replication
} ## end of the loop over the first XC factor
} ## end of the loop over the second XC factor

###---------         Export output in a text file                 ---------###

finaloutput<-as.data.frame(round(finaloutput,3))
output<-data.frame(cbind(rep(n2range,each=n1size*repsize),rep(n1range,n2size,each=repsize),rep(reprange,n1size*n2size),finaloutput))
names(output)<-c("#XC2","#XC1","#1-level","zLb0","zpb0","zUb0","sLb0","spb0","sUb0","zLb1","zpb1","zUb1","sLb1","spb1","sUb1")
write.table(output,"power.txt",sep="\t",quote=F,eol="\n",dec=".",col.names=T,row.names=F,qmethod="double")