||
##################################################
################## power4SEM ####################
########### created by Suzanne Jak ###############
##################################################
# Load packages
if(!require(shiny)){install.packages('shiny')}
if(!require(semTools)){install.packages('semTools')}
if(!require(lavaan)){install.packages('lavaan')}
if(!require(semPlot)){install.packages('semPlot')}
if(!require(shinyAce)){install.packages('shinyAce')}
if(!require(shinyhelper)){install.packages('shinyhelper')}
if(!require(shinycssloaders)){install.packages('shinycssloaders')}
if(!require(shinyBS)){install.packages('shinyBS')}
library(shiny)
library(semTools)
library(lavaan)
library(semPlot)
library(shinyAce)
library(shinyhelper)
library(shinycssloaders)
library(shinyBS)
#source("http://www.suzannejak.nl/SEM_power_functions.R")
#define power function powerSS
SEM_chipower <- function(ncp,df,a = 0.05,plotpower=TRUE){
crit <- qchisq(1-a,df)
power <- 1 - pchisq(crit,df,ncp)
if(plotpower==TRUE){
max <- qchisq(.999,df,ncp)
x <- seq(0,max,.01)
dens0 <- dchisq(x,df)
dens1 <- dchisq(x,df,ncp)
powerplot <- plot(x,dens0, type = "l", col = "red", lty = "dashed", ylab = "density", xlab = "chi-square")
lines(x,dens1,type = "l", col = "blue")
i <- x>=crit & x <= max # select those x's over critical value
polygon(c(crit,x[i],max), c(0,dens1[i],0), col="blue", density = 20)
abline(v=crit, lty = "twodash")
legend("topright", inset=.05, title="Legend",
c("H0 is True","H1 is True", expression(paste("Critical ",chi^2))), lwd=2, lty=c("dashed","solid","twodash"), col=c("red","blue","black"))
title(main = paste("Power =", round(power,3)))
}
return(power)
}
SEM_NforChipower <- function(powd,ncp,df,N,alpha = .05){
Fmin <- ncp/(N-1)
powa <- SEM_chipower(ncp,df,alpha,plotpower=FALSE)
# add 100 to N until actual power is higher than desired power
while (powa < powd){
N <- N + 100
ncp.n <- (N-1) * Fmin
powa <- SEM_chipower(ncp.n,df,alpha,plotpower=FALSE)
}
# half interval, add or substract from N, re-calculate power, until powdif < .000001.
intv <- 100
powdif <- max(powa,powd) - min(powa,powd)
while(powdif > .000001){
if(powa > powd) dir = -1 else dir=1
intv <- intv * .5
N <- N + dir* intv
ncp.n <- (N-1) * Fmin
powa <- SEM_chipower(ncp.n,df,alpha,plotpower=FALSE)
powdif <- max(powa,powd) - min(powa,powd)
}
return(paste0("For a power of ",powd,", the minimum sample size needed is ",ceiling(N)," (NCP = ",round(ncp.n,3),")"))
}
SEM_RMSEApower <- function(RMSEA0,RMSEA1,df,N,alpha = .05){
lambda0 <- (N-1) * df * (RMSEA0)^2
lambda1 <- (N-1) * df * (RMSEA1)^2
max <- qchisq(.999,df,max(lambda0,lambda1))
x <- seq(0,max,.1)
if (RMSEA0 < RMSEA1){
crit <- qchisq(1-alpha,df,lambda0)
RMSEApower <- 1 - pchisq(crit,df,lambda1)
dens0 <- dchisq(x,df,lambda0)
dens1 <- dchisq(x,df,lambda1)
plot <- plot(x,dens0, type = "l", col = "red", lty = "dashed", ylab = "density", xlab = "chi-square", ylim = c(0,max(dens0)))
lines(x,dens1,type = "l", col = "blue")
i <- x>=crit & x <= max # select those x's over critical value
polygon(c(crit,x[i],max), c(0,dens1[i],0), col="blue", density = 20)
}
if (RMSEA1 < RMSEA0){
crit <- qchisq(alpha,df,lambda0)
RMSEApower <- pchisq(crit,df,lambda1)
dens0 <- dchisq(x,df,lambda0)
dens1 <- dchisq(x,df,lambda1)
plot <- plot(x,dens0, type = "l", col = "red", lty = "dashed", ylab = "density", xlab = "chi-square", ylim = c(0,max(dens1)))
lines(x,dens1,type = "l", col = "blue")
i <- x<=crit # select those x's over critical value
polygon(c(0,x[i],crit), c(0,dens1[i],0), col="blue", density = 20)
}
abline(v=crit, lty = "twodash")
legend("topright", inset=.05, title="Legend",
c(paste("RMSEA = ",RMSEA0),paste("RMSEA = ",RMSEA1),expression(paste("Critical ",chi^2))), lwd=2, lty=c("dashed","solid","twodash"), col=c("red","blue","black"))
title(main = paste("Power to reject RMSEA of", RMSEA0, "=", round(RMSEApower,3)))
return(round(RMSEApower,3))
}
# powerSS() is an adaptation of SSpower() from semTools
powerSS <- function (powerModel, n, nparam, popModel, mu, Sigma, fun = "cfa",
alpha = 0.05, ...)
{
if (missing(Sigma)) {
popMoments <- lavaan::fitted(do.call(fun, list(model = popModel)))
if (length(n) > 1L) {
Sigma <- list()
mu <- if (!is.null(popMoments[[1]]$mean))
list()
else NULL
for (g in 1:length(n)) {
Sigma[[g]] <- popMoments$cov
if (!is.null(popMoments[[g]]$mean))
mu[[g]] <- popMoments$mean
}
}
else {
Sigma <- popMoments$cov
mu <- popMoments$mean
}
}
else {
if (is.list(Sigma)) {
nG <- length(Sigma)
if (length(n) < nG)
n <- rep(n, length.out = nG)
if (length(n) > nG)
n <- n[1:nG]
no.mu <- missing(mu)
if (!no.mu)
null.mu <- any(sapply(mu, is.null))
if (no.mu || null.mu) {
mu <- NULL
}
}
else if (is.matrix(Sigma)) {
n <- n[[1]]
if (missing(mu)) {
mu <- NULL
}
else if (!is.numeric(mu) || !!is.null(mu)) {
stop("mu must be a numeric vector, or a list (one vector per group)")
}
}
else stop("Sigma must be a covariance matrix, or a list (one matrix per group)")
}
dots <- list(...)
funArgs <- list(model = powerModel, sample.cov = Sigma, sample.mean = mu,
sample.nobs = n)
useArgs <- c(funArgs, dots[setdiff(names(dots), names(funArgs))])
fit <- do.call(fun, useArgs)
ncp <- lavaan::fitmeasures(fit)["chisq"]
critVal <- qchisq(alpha, df = nparam, lower.tail = FALSE)
list(power = pchisq(critVal, df = nparam, ncp = ncp, lower.tail = FALSE),ncp = round(ncp,3),n = n, popMoments=popMoments)
}
# Define UI ----
ui <- navbarPage(title="Power calculations for SEM",
tabPanel(tags$h4(style = "font-family:Impact; color: green","lavaan input"),
fluidRow(
column(6,
actionButton("fitButton", "Obtain NCP",
style="color: #fff; background-color: green; border-color: green"),
p("Specify the intended sample size, the H1 model and the H0 model below and click the green button to obtain the noncentrality parameter (NCP)"),
shinyhelper::helper(numericInput(inputId = "Nint", "Intended sample size",200),
type = "inline",
buttonLabel = "Close",
title = "Intended sample size",
content = c("Insert the intended sample size")),
p('See',a("lavaan", href="http://lavaan.ugent.be/", target="_blank"),'for instructions about model specification.'),
shinyhelper::helper(p(),
type = "inline",
title = "Model under H1",
buttonLabel = "Close",
content = c("Specify the model under the alternative hypothesis using lavaan syntax.",
"All parameters should be specified as fixed.",
"This is the model that includes ALL parameters."),
size = "s"),
textAreaInput("h1model_lavaan",
label = "Specify the model under H1 (the model with all fixed population values)",
width = "600px",
height = "300px",
value =
"# Regression coefficients
y7 ~ .21*y1
y7 ~ .05*y2
y7 ~ .08*y3
y7 ~ .60*y4
y7 ~ -.13*y5
y7 ~ .12*y6
y8 ~ .80*y7
# Covariances
y1 ~~ .46*y2
y1 ~~ .46*y3
y1 ~~ .51*y4
y1 ~~ .50*y5
y1 ~~ .40*y6
y2 ~~ .28*y3
y2 ~~ .39*y4
y2 ~~ .40*y5
y2 ~~ .27*y6
y3 ~~ .78*y4
y3 ~~ .71*y5
y3 ~~ .54*y6
y4 ~~ .79*y5
y4 ~~ .66*y6
y5 ~~ .78*y6
# Variances
y1 ~~ 1*y1;
y2 ~~ 1*y2;
y1 ~~ 1*y3;
y1 ~~ 1*y4;
y1 ~~ 1*y5;
y1 ~~ 1*y6;
y7 ~~ .36*y7
y8 ~~ .36*y8
# Extra parameters
y8 ~ .10*y1
y8 ~ .10*y2
y8 ~ .10*y3
y8 ~ .10*y4
y8 ~ .10*y5
y8 ~ .10*y6"),
actionButton("View", "View H1 values",
style="color: #fff; background-color: green; border-color: green"),
bsModal("ViewSigma", "H1 model implied covariance matrix and standardized parameter values", "View",size = "large",
"Model implied covariance matrix from model H1.",
verbatimTextOutput("sigmaH1"),
"Variances of the variables are on the diagonal.",br(),
"If all model implied variances are 1, then the population values are in standardized metric already.",br(),
"If not, see below for the standardized population values.",
verbatimTextOutput("stdH1")),
tags$head(tags$style(HTML('
textArea {
background-color: #fbfcfa !important;
font-family:Courier;
border: none;
}')))
),
column(6,
shinyhelper::helper(tags$h4(style = "font-family:Impact; color: green","Result"),
type = "inline",
title = "Obtain NCP",
buttonLabel = "Close",
content = c("If you see an error here, then this is the exact error message as given by lavaan.",
"It may be easiest to find out what went wrong by specifying your models in lavaan directly"),
size = "s"),
wellPanel(textOutput("lavaanNCP"),
tags$style("#lavaanNCP {font-size:14px;
color:green;
display:block; }")),
p("H1 model (all paths should be black because they are fixed)"),
plotOutput("ModelH1plot"))),
fluidRow(
column(6,
shinyhelper::helper(p(),
type = "inline",
title = "Model under H0",
buttonLabel = "Close",
content = c("Specify the model under the null hypothesis using lavaan syntax.",
"This is the model that specifies less parameters than the one above."),
size = "s"),
textAreaInput("h0model_lavaan",
label = "Specify the model under H0",
width = "600px",
height = "300px",
value =
"
# Regression coefficients
y7 ~ y1
y7 ~ y2
y7 ~ y3
y7 ~ y4
y7 ~ y5
y7 ~ y6
y8 ~ y7
# Covariances
y1 ~~ y2
y1 ~~ y3
y1 ~~ y4
y1 ~~ y5
y1 ~~ y6
y2 ~~ y3
y2 ~~ y4
y2 ~~ y5
y2 ~~ y6
y3 ~~ y4
y3 ~~ y5
y3 ~~ y6
y4 ~~ y5
y4 ~~ y6
y5 ~~ y6
# Variances
y1 ~~ y1;
y2 ~~ y2;
y3 ~~ y3;
y4 ~~ y4;
y5 ~~ y5;
y6 ~~ y6;
y7 ~~ y7
y8 ~~ y8")),
column(6,
p("H0 model"),
plotOutput("ModelH0plot")
)
)
),
#### CHI SQUARE PAGE
tabPanel(tags$h4(style = "font-family:Impact; color: blue","Chi-square test"),
fluidRow(
column(4,
wellPanel(
tags$h4(style = "font-family:Impact; color: blue","Input"),
textOutput("calculatedNCP"),
tags$style("#calculatedNCP {font-size:12px;
color:green;
display:block; }"),
shinyhelper::helper(numericInput(inputId = "ncp", "Noncentrality parameter",NA,step = 0.1),
type = "inline",
title = "Noncentrality parameter",
buttonLabel = "Close",
content = c("Insert the noncentrality parameter, as obtained by fitting the H0 model to population data generated from H1.",
"You can use the green tab (lavaan-input) to obtain the noncentrality parameter using lavaan-syntax."),
size = "s"),
shinyhelper::helper(numericInput(inputId = "df", "Degrees of freedom",NA),
type = "inline",
title = "Degrees of Freedom",
buttonLabel = "Close",
content = c("Insert the degrees of freedom associated with the chi-square test.",
"For overall model fit, the df are equal to the number of observed statistics minus the number of freely estimated parameters in H0.",
"For the chi-square difference test, the df is equal to the difference in number of parameters between model H0 and model H1."),
size = "s"),
shinyhelper::helper(numericInput(inputId = "a", HTML("α"),NA,step = 0.01),
type = "inline",
title = "Significance level",
buttonLabel = "Close",
content = c("Choose a significance level for the test."),
size = "s"),
actionButton("hitButton", "Calculate!",
style="color: #fff; background-color: blue; border-color: blue")
)),
column(6,
tags$h4(style = "font-family:Impact; color: blue","Result"),
textOutput("chi_power"),
plotOutput("plotpower")
)),
fluidRow(
column(4,
wellPanel(
tags$h4(style = "font-family:Impact; color: blue","Calculate minimum sample size for desired power"),
shinyhelper::helper(numericInput(inputId = "N", "Sample size used to obtain noncentrality parameter",NA),
type = "inline",
title = "N",
buttonLabel = "Close",
content = c("Insert the sample size that was use to obtain the noncentrality parameter that was inserted above."),
size = "s"),
shinyhelper::helper(numericInput(inputId = "powd", "Desired power",NA, step = 0.05),
type = "inline",
title = "Desired power level",
buttonLabel = "Close",
content = c("Specify the power level for which you want to know the minimal sample size."),
size = "s"),
actionButton("goButton", "Calculate!",
style="color: #fff; background-color: blue; border-color: blue")
)),
column(6,
tags$h4(style = "font-family:Impact; color: blue","Result"),
withSpinner(textOutput("samplesize_needed"))))
),
### RMSEA PAGE
tabPanel(tags$h4(style = "font-family:Impact; color: red","RMSEA"),
fluidRow(
column(4,
wellPanel(
tags$h4(style = "font-family:Impact; color: red","Calculate RMSEA-based power"),
shinyhelper::helper(numericInput(inputId = "RMSEA0", HTML("RMSEA H0"), NA, step = 0.01),
type = "inline",
title = "RMSEA under H0",
buttonLabel = "Close",
content = c("RMSEA value under H0. Following the guidelines of MacCallum et al., this is often .05."),
size = "s"),
shinyhelper::helper(numericInput(inputId = "RMSEA1", "RMSEA H1", NA, step = 0.01),
type = "inline",
title = "Inline Help",
buttonLabel = "Close",
content = c("RMSEA value under H1.",
"Following the guidelines of MacCallum et al., this is .08 for a test of close fit, and .01 for a test of not-close fit."),
size = "s"),
shinyhelper::helper(numericInput(inputId = "df_rmsea", "df", NA),
type = "inline",
buttonLabel = "Close",
title = "Degrees of Freedom",
content = c("The degrees of freedom of the model of interest.",
"Degrees of freedom are equal to the number of observed statistics minus the number of freely estimated model parameters"),
size = "s"),
shinyhelper::helper(numericInput(inputId = "N_rmsea", "N", NA),
type = "inline",
title = "Sample size",
buttonLabel = "Close",
content = c("The sample size for which you want to calculate the power."),
size = "s"),
shinyhelper::helper(numericInput(inputId = "a_rmsea", HTML("α"),NA, step = 0.01),
type = "inline",
title = "Significance level",
buttonLabel = "Close",
content = c("Choose a significance level for the test."),
size = "s"),
actionButton("rmseaButton", "Calculate!",
style="color: #fff; background-color: red; border-color: red")
)),
column(6,
tags$h4(style = "font-family:Impact; color: red","Result"),
plotOutput("plotRMSEApower"))),
fluidRow(
column(4,
wellPanel(
tags$h4(style = "font-family:Impact; color: red","Calculate required sample size for desired power"),
shinyhelper::helper(numericInput(inputId = "dpower", "Desired power", NA,step = 0.05),
type = "inline",
title = "Desired power",
buttonLabel = "Close",
content = c("Specify the desired power level for which you want to know the minimal sample size."),
size = "s"),
actionButton("NrmseaButton", "Calculate!",
style="color: #fff; background-color: red; border-color: red")
)),
column(6,
tags$h4(style = "font-family:Impact; color: red","Result"),
withSpinner(textOutput("RMSEA.samplesize.needed"))))
),
### Documentation page
tabPanel(tags$h4(style = "font-family:Impact; color: Grey","Documentation"),
p("For the tutorial paper about this app see Jak, Jorgensen, Verdam, Oort and Elffers (2020) (not added yet)", a("http://www.suzannejak.nl", href="http://www.suzannejak.nl", target="_blank"))
))
# Define server logic ----
server <- function(input, output) {
# calculate NCP with SSpower() (based on semTools function)
pSS <- eventReactive({input$fitButton}, {powerSS(powerModel=input$h0model_lavaan, popModel=input$h1model_lavaan,n=input$Nint,nparam=1)})
checkH1 <- eventReactive({input$ViewSigma}, {
popModel <- input$h1model_lavaan
fitH1 <- lavaan::cfa(model=popModel)
sigmaH1 <- lavaan::fitted(fitH1)$cov
std <- lavaan::standardizedSolution(fitH1, cov.std = FALSE)
list(StandardizedParameterValues = std, ModelImpliedCovariances = sigmaH1)
})
calcNCP <- reactive({pSS()$ncp})
output$lavaanNCP <- renderText({paste0("The noncentrality parameter obtained by fitting the lavaan-models equals ",calcNCP())})
output$calculatedNCP <- renderText({paste0("(The NCP obtained by fitting the lavaan models in the previous tab is ",calcNCP(),")")})
calc.chi.power <- eventReactive(input$hitButton, {
paste0("With noncentrality parameter = ",input$ncp,", df = ",input$df,", and alpha = ",input$a,", the statistical power is ",round(SEM_chipower(ncp=input$ncp,df=input$df,a=input$a,plotpower=FALSE),3))
})
plot.chi.power <- eventReactive(input$hitButton, {
SEM_chipower(ncp=input$ncp,df=input$df,a=input$a)
})
output$chi_power <- renderText({calc.chi.power()})
output$plotpower <- renderPlot({plot.chi.power()})
sampneed <- eventReactive(input$goButton, {
SEM_NforChipower(powd=input$powd,ncp=input$ncp,df=input$df,N=input$N,alpha=input$a)
})
output$samplesize_needed <- renderText({sampneed()})
calc.rmsea.power <- eventReactive(input$rmseaButton, {
SEM_RMSEApower(input$RMSEA0,input$RMSEA1,input$df_rmsea,input$N_rmsea,input$a_rmsea)
})
output$plotRMSEApower <- renderPlot({calc.rmsea.power()})
calc.n.rmsea <- eventReactive(input$NrmseaButton, {
paste0("For a power of ",input$dpower,", the minimum sample size needed is ",semTools::findRMSEAsamplesize(input$RMSEA0,input$RMSEA1,input$df_rmsea,input$dpower, input$a_rmsea))
})
output$RMSEA.samplesize.needed <- renderText({calc.n.rmsea()})
output$ModelH1plot <- renderPlot({semPlot::semPaths(semPlot::semPlotModel(input$h1model_lavaan, fixed.x=FALSE),
nCharNodes=0, nCharEdges=0,
layout="tree", sizeMan=8,
sizeLat=8, edge.label.cex=1.3,
color="white",freeStyle = c("red",1),
fixedStyle = c("black",1))})
output$ModelH0plot <- renderPlot({semPlot::semPaths(semPlot::semPlotModel(input$h0model_lavaan, fixed.x=FALSE),
nCharNodes=0, nCharEdges=0,
layout="tree", sizeMan=8,
sizeLat=8, edge.label.cex=1.3,
color="white",freeStyle = c("red",1),
fixedStyle = c("black",1))})
output$sigmaH1 <- renderPrint({checkH1()$ModelImpliedCovariances})
output$stdH1 <- renderPrint({checkH1()$StandardizedParameterValues})
shinyhelper::observe_helpers()
}
# Run the app ----
shinyApp(ui = ui, server = server)
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-7-27 22:42
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社