蕭文龍 部落格 Mac Shiau(blog)

把我珍藏的影片更新如下，需要的可以看看，謝謝：
Introduction of Reflective model PLS-SEM:

把我珍藏的影片更新如下，需要的可以看看，謝謝：
20200317 Introduce Mac Shiau:

[元宇宙在中國大爆發]
在國家政策扶持下，各省全力支持元宇宙產業，虛實整合量表協助您加速研究發表期刊論文

【Knowledge】What to do if the AVE is large than the shared variance
請問如果AVE大於shared variance怎麼辦呢？

元宇宙下有那些應用場景的研究？ 虛實整合量表協助您加速研究設計發表期刊論文
浙江工業大學 蕭文龍團隊 整理

資料來源：Aguinis, H., & Cronin, M. A. (2022). It’s the Theory, Stupid. Organizational Psychology Review. https://doi.org/10.1177/20413866221080629
浙江工業大學 蕭文龍團隊 整理

Why a reviewer asks for random? Why we need to report effective size: f square?
[溫馨回答]：

Surprised：
     有資深教授說： PLS已實質含蓋CB-SEM功能。 統計分析上，PLS功能及實用性已大於CB-SEM等......是不精確的知識。

SSpower <- function(modelpop,modelana,n,alpha=.05,anadf=1){
  require(lavaan)

library(lavaan)
popModel <- "

##################################################
################## 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)
 
 

編碼如下：
library(devtools)