##################################################

################## 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("&alpha;"),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("&alpha;"),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)