Using simsem to conduct Power Analysis for Structural Equational Modeling
Here is an example of conducting a power analysis for Structural Equation Modeling (SEM) using the simsem program in R. In this example, we estimate the amount of data needed to conduct a SEM with three factors with three indicators, and two path coefficients measuring the direct effect that factor three has on factor one and two.
First, we need to library the simsem program. Then we need to create a matrix for our indicators. Here we are creating a matrix with three unique indicators for three factors meaning that we need a total of nine indicators for the three factors.
For the loadings that we want to assign values to later, we assign them to the value NA for now.
Then we create the starting or initial values in the loading.start matrix. We assign the loading values of .7 to each of the indicators for each of the factors. These starting values can vary depending upon the context of your research.
Next, we need to assign the cross loading among loadings with different factors. We could assign them to zero and assume that there is no correlation between loadings of different factors; however, below, we have chosen to model the correlation between the loading as essentially random with a uniform distribution.
Finally, we combine the loading matrix, the starting values, and the cross-loading correlations matrices into one loading matrix titled LY
library(simsem)
loading <- matrix(0, 9, 3)
loading[1:3, 1] <- NA
loading[4:6, 2] <- NA
loading[7:9, 3] <- NA; loading## [,1] [,2] [,3]
## [1,] NA 0 0
## [2,] NA 0 0
## [3,] NA 0 0
## [4,] 0 NA 0
## [5,] 0 NA 0
## [6,] 0 NA 0
## [7,] 0 0 NA
## [8,] 0 0 NA
## [9,] 0 0 NAloading.start <- matrix("", 9, 3)
loading.start[1:3, 1] <- 0.7
loading.start[4:6, 2] <- 0.7
loading.start[7:9, 3] <- 0.7
#Cross loading or loadings among indicators of different factors are random instead of zero
loading.trivial <- matrix("runif(1, -0.2, 0.2)", 9, 3)
loading.trivial[is.na(loading)] <- 0
LY <- bind(loading, loading.start, misspec=loading.trivial); LY## [,1] [,2]
## [1,] "NA:0.7+0" "0+runif(1, -0.2, 0.2)"
## [2,] "NA:0.7+0" "0+runif(1, -0.2, 0.2)"
## [3,] "NA:0.7+0" "0+runif(1, -0.2, 0.2)"
## [4,] "0+runif(1, -0.2, 0.2)" "NA:0.7+0"
## [5,] "0+runif(1, -0.2, 0.2)" "NA:0.7+0"
## [6,] "0+runif(1, -0.2, 0.2)" "NA:0.7+0"
## [7,] "0+runif(1, -0.2, 0.2)" "0+runif(1, -0.2, 0.2)"
## [8,] "0+runif(1, -0.2, 0.2)" "0+runif(1, -0.2, 0.2)"
## [9,] "0+runif(1, -0.2, 0.2)" "0+runif(1, -0.2, 0.2)"
## [,3]
## [1,] "0+runif(1, -0.2, 0.2)"
## [2,] "0+runif(1, -0.2, 0.2)"
## [3,] "0+runif(1, -0.2, 0.2)"
## [4,] "0+runif(1, -0.2, 0.2)"
## [5,] "0+runif(1, -0.2, 0.2)"
## [6,] "0+runif(1, -0.2, 0.2)"
## [7,] "NA:0.7+0"
## [8,] "NA:0.7+0"
## [9,] "NA:0.7+0"Next, we have to create the error variables for each of the indicators. We simply assign a random error using and random normal distribution for each of the error indicators and assign them a variance of one. We then combine the error variables for the indicator variables with the variances for the error variable into an error indicator matrix called RTE.
error.cor.trivial <- matrix("rnorm(1, 0, 0.1)", 9, 9)
diag(error.cor.trivial) <- 1
RTE <- binds(diag(9), misspec=error.cor.trivial); RTE## [,1] [,2] [,3]
## [1,] "1+1" "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)"
## [2,] "0+rnorm(1, 0, 0.1)" "1+1" "0+rnorm(1, 0, 0.1)"
## [3,] "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)" "1+1"
## [4,] "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)"
## [5,] "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)"
## [6,] "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)"
## [7,] "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)"
## [8,] "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)"
## [9,] "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)"
## [,4] [,5] [,6]
## [1,] "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)"
## [2,] "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)"
## [3,] "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)"
## [4,] "1+1" "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)"
## [5,] "0+rnorm(1, 0, 0.1)" "1+1" "0+rnorm(1, 0, 0.1)"
## [6,] "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)" "1+1"
## [7,] "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)"
## [8,] "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)"
## [9,] "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)"
## [,7] [,8] [,9]
## [1,] "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)"
## [2,] "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)"
## [3,] "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)"
## [4,] "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)"
## [5,] "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)"
## [6,] "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)"
## [7,] "1+1" "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)"
## [8,] "0+rnorm(1, 0, 0.1)" "1+1" "0+rnorm(1, 0, 0.1)"
## [9,] "0+rnorm(1, 0, 0.1)" "0+rnorm(1, 0, 0.1)" "1+1"Next we need to assign the correlation between the factors. For simplicity, we are assuming that all the correlation between the factors are zero.
factor.cor <- diag(3)
RPS <- binds(factor.cor, 0.0); RPS## [,1] [,2] [,3]
## [1,] "1" "0" "0"
## [2,] "0" "1" "0"
## [3,] "0" "0" "1"Next, we need to create the path coefficients matrix. Below we have two paths that are not zero. These are paths from factor three to factor and two, which we initially set to NA. We have assigned these paths to be from a normal and uniform distribution with initial parameters shown below. Then we assign those initial values to the path coefficient matrix by combining them into one path coefficient matrix called BE.
path <- matrix(0, 3, 3)
path[3, 1:2] <- NA
path.start <- matrix(0, 3, 3)
path.start[3, 1] <- "rnorm(1, 0.6, 0.05)"
path.start[3, 2] <- "runif(1, 0.3, 0.5)"
BE <- bind(path, path.start); BE## [,1] [,2] [,3]
## [1,] "0" "0" "0"
## [2,] "0" "0" "0"
## [3,] "NA:rnorm(1, 0.6, 0.05)" "NA:runif(1, 0.3, 0.5)" "0"Finally, we can create the sem model that we will simulate to find the appropriate power. We combine can each of the pieces that we developed earlier in the example, into one model titled SEM.model.
Then we simulate this model with 50 to 500 participants. We then plot the output.
We can evaluate the simulation by reviewing graphs. As we can see as n gets higher, the values for the sample size sensitive statistics increase. For example, the chi-square, the statistic becomes larger with larger samples sizes. We can see for statistics like the RMSEA, that this value hovers around .1, which relatively low.
Finally, we can review what we want, which is the sample size requirements to conduct this SEM project. We can look at the values for the two path coefficients regressed upon factor 1 and factor 2. We can see that a sample size of 93 is required for this model, because that is the largest sample requirement.
SEM.model <- model(BE=BE, LY=LY, RPS=RPS, RTE=RTE, modelType="SEM")
Output <- sim(NULL, n=50:500, SEM.model) ## Progress: 1 / 451
## Progress: 2 / 451
## Progress: 3 / 451
## Progress: 4 / 451
## Progress: 5 / 451
## Progress: 6 / 451
## Progress: 7 / 451
## Progress: 8 / 451
## Progress: 9 / 451
## Progress: 10 / 451
## Progress: 11 / 451
## Progress: 12 / 451
## Progress: 13 / 451
## Progress: 14 / 451
## Progress: 15 / 451
## Progress: 16 / 451
## Progress: 17 / 451
## Progress: 18 / 451
## Progress: 19 / 451
## Progress: 20 / 451
## Progress: 21 / 451
## Progress: 22 / 451
## Progress: 23 / 451
## Progress: 24 / 451
## Progress: 25 / 451
## Progress: 26 / 451
## Progress: 27 / 451
## Progress: 28 / 451
## Progress: 29 / 451
## Progress: 30 / 451
## Progress: 31 / 451
## Progress: 32 / 451
## Progress: 33 / 451
## Progress: 34 / 451
## Progress: 35 / 451
## Progress: 36 / 451
## Progress: 37 / 451
## Progress: 38 / 451
## Progress: 39 / 451
## Progress: 40 / 451
## Progress: 41 / 451
## Progress: 42 / 451
## Progress: 43 / 451
## Progress: 44 / 451
## Progress: 45 / 451
## Progress: 46 / 451
## Progress: 47 / 451
## Progress: 48 / 451
## Progress: 49 / 451
## Progress: 50 / 451
## Progress: 51 / 451
## Progress: 52 / 451
## Progress: 53 / 451
## Progress: 54 / 451
## Progress: 55 / 451
## Progress: 56 / 451
## Progress: 57 / 451
## Progress: 58 / 451
## Progress: 59 / 451
## Progress: 60 / 451
## Progress: 61 / 451
## Progress: 62 / 451
## Progress: 63 / 451
## Progress: 64 / 451
## Progress: 65 / 451
## Progress: 66 / 451
## Progress: 67 / 451
## Progress: 68 / 451
## Progress: 69 / 451
## Progress: 70 / 451
## Progress: 71 / 451
## Progress: 72 / 451
## Progress: 73 / 451
## Progress: 74 / 451
## Progress: 75 / 451
## Progress: 76 / 451
## Progress: 77 / 451
## Progress: 78 / 451
## Progress: 79 / 451
## Progress: 80 / 451
## Progress: 81 / 451
## Progress: 82 / 451
## Progress: 83 / 451
## Progress: 84 / 451
## Progress: 85 / 451
## Progress: 86 / 451
## Progress: 87 / 451
## Progress: 88 / 451
## Progress: 89 / 451
## Progress: 90 / 451
## Progress: 91 / 451
## Progress: 92 / 451
## Progress: 93 / 451
## Progress: 94 / 451
## Progress: 95 / 451
## Progress: 96 / 451
## Progress: 97 / 451
## Progress: 98 / 451
## Progress: 99 / 451
## Progress: 100 / 451
## Progress: 101 / 451
## Progress: 102 / 451
## Progress: 103 / 451
## Progress: 104 / 451
## Progress: 105 / 451
## Progress: 106 / 451
## Progress: 107 / 451
## Progress: 108 / 451
## Progress: 109 / 451
## Progress: 110 / 451
## Progress: 111 / 451
## Progress: 112 / 451
## Progress: 113 / 451
## Progress: 114 / 451
## Progress: 115 / 451
## Progress: 116 / 451
## Progress: 117 / 451
## Progress: 118 / 451
## Progress: 119 / 451
## Progress: 120 / 451
## Progress: 121 / 451
## Progress: 122 / 451
## Progress: 123 / 451
## Progress: 124 / 451
## Progress: 125 / 451
## Progress: 126 / 451
## Progress: 127 / 451
## Progress: 128 / 451
## Progress: 129 / 451
## Progress: 130 / 451
## Progress: 131 / 451
## Progress: 132 / 451
## Progress: 133 / 451
## Progress: 134 / 451
## Progress: 135 / 451
## Progress: 136 / 451
## Progress: 137 / 451
## Progress: 138 / 451
## Progress: 139 / 451
## Progress: 140 / 451
## Progress: 141 / 451
## Progress: 142 / 451
## Progress: 143 / 451
## Progress: 144 / 451
## Progress: 145 / 451
## Progress: 146 / 451
## Progress: 147 / 451
## Progress: 148 / 451
## Progress: 149 / 451
## Progress: 150 / 451
## Progress: 151 / 451
## Progress: 152 / 451
## Progress: 153 / 451
## Progress: 154 / 451
## Progress: 155 / 451
## Progress: 156 / 451
## Progress: 157 / 451
## Progress: 158 / 451
## Progress: 159 / 451
## Progress: 160 / 451
## Progress: 161 / 451
## Progress: 162 / 451
## Progress: 163 / 451
## Progress: 164 / 451
## Progress: 165 / 451
## Progress: 166 / 451
## Progress: 167 / 451
## Progress: 168 / 451
## Progress: 169 / 451
## Progress: 170 / 451
## Progress: 171 / 451
## Progress: 172 / 451
## Progress: 173 / 451
## Progress: 174 / 451
## Progress: 175 / 451
## Progress: 176 / 451
## Progress: 177 / 451
## Progress: 178 / 451
## Progress: 179 / 451
## Progress: 180 / 451
## Progress: 181 / 451
## Progress: 182 / 451
## Progress: 183 / 451
## Progress: 184 / 451
## Progress: 185 / 451
## Progress: 186 / 451
## Progress: 187 / 451
## Progress: 188 / 451
## Progress: 189 / 451
## Progress: 190 / 451
## Progress: 191 / 451
## Progress: 192 / 451
## Progress: 193 / 451
## Progress: 194 / 451
## Progress: 195 / 451
## Progress: 196 / 451
## Progress: 197 / 451
## Progress: 198 / 451
## Progress: 199 / 451
## Progress: 200 / 451
## Progress: 201 / 451
## Progress: 202 / 451
## Progress: 203 / 451
## Progress: 204 / 451
## Progress: 205 / 451
## Progress: 206 / 451
## Progress: 207 / 451
## Progress: 208 / 451
## Progress: 209 / 451
## Progress: 210 / 451
## Progress: 211 / 451
## Progress: 212 / 451
## Progress: 213 / 451
## Progress: 214 / 451
## Progress: 215 / 451
## Progress: 216 / 451
## Progress: 217 / 451
## Progress: 218 / 451
## Progress: 219 / 451
## Progress: 220 / 451
## Progress: 221 / 451
## Progress: 222 / 451
## Progress: 223 / 451
## Progress: 224 / 451
## Progress: 225 / 451
## Progress: 226 / 451
## Progress: 227 / 451
## Progress: 228 / 451
## Progress: 229 / 451
## Progress: 230 / 451
## Progress: 231 / 451
## Progress: 232 / 451
## Progress: 233 / 451
## Progress: 234 / 451
## Progress: 235 / 451
## Progress: 236 / 451
## Progress: 237 / 451
## Progress: 238 / 451
## Progress: 239 / 451
## Progress: 240 / 451
## Progress: 241 / 451
## Progress: 242 / 451
## Progress: 243 / 451
## Progress: 244 / 451
## Progress: 245 / 451
## Progress: 246 / 451
## Progress: 247 / 451
## Progress: 248 / 451
## Progress: 249 / 451
## Progress: 250 / 451
## Progress: 251 / 451
## Progress: 252 / 451
## Progress: 253 / 451
## Progress: 254 / 451
## Progress: 255 / 451
## Progress: 256 / 451
## Progress: 257 / 451
## Progress: 258 / 451
## Progress: 259 / 451
## Progress: 260 / 451
## Progress: 261 / 451
## Progress: 262 / 451
## Progress: 263 / 451
## Progress: 264 / 451
## Progress: 265 / 451
## Progress: 266 / 451
## Progress: 267 / 451
## Progress: 268 / 451
## Progress: 269 / 451
## Progress: 270 / 451
## Progress: 271 / 451
## Progress: 272 / 451
## Progress: 273 / 451
## Progress: 274 / 451
## Progress: 275 / 451
## Progress: 276 / 451
## Progress: 277 / 451
## Progress: 278 / 451
## Progress: 279 / 451
## Progress: 280 / 451
## Progress: 281 / 451
## Progress: 282 / 451
## Progress: 283 / 451
## Progress: 284 / 451
## Progress: 285 / 451
## Progress: 286 / 451
## Progress: 287 / 451
## Progress: 288 / 451
## Progress: 289 / 451
## Progress: 290 / 451
## Progress: 291 / 451
## Progress: 292 / 451
## Progress: 293 / 451
## Progress: 294 / 451
## Progress: 295 / 451
## Progress: 296 / 451
## Progress: 297 / 451
## Progress: 298 / 451
## Progress: 299 / 451
## Progress: 300 / 451
## Progress: 301 / 451
## Progress: 302 / 451
## Progress: 303 / 451
## Progress: 304 / 451
## Progress: 305 / 451
## Progress: 306 / 451
## Progress: 307 / 451
## Progress: 308 / 451
## Progress: 309 / 451
## Progress: 310 / 451
## Progress: 311 / 451
## Progress: 312 / 451
## Progress: 313 / 451
## Progress: 314 / 451
## Progress: 315 / 451
## Progress: 316 / 451
## Progress: 317 / 451
## Progress: 318 / 451
## Progress: 319 / 451
## Progress: 320 / 451
## Progress: 321 / 451
## Progress: 322 / 451
## Progress: 323 / 451
## Progress: 324 / 451
## Progress: 325 / 451
## Progress: 326 / 451
## Progress: 327 / 451
## Progress: 328 / 451
## Progress: 329 / 451
## Progress: 330 / 451
## Progress: 331 / 451
## Progress: 332 / 451
## Progress: 333 / 451
## Progress: 334 / 451
## Progress: 335 / 451
## Progress: 336 / 451
## Progress: 337 / 451
## Progress: 338 / 451
## Progress: 339 / 451
## Progress: 340 / 451
## Progress: 341 / 451
## Progress: 342 / 451
## Progress: 343 / 451
## Progress: 344 / 451
## Progress: 345 / 451
## Progress: 346 / 451
## Progress: 347 / 451
## Progress: 348 / 451
## Progress: 349 / 451
## Progress: 350 / 451
## Progress: 351 / 451
## Progress: 352 / 451
## Progress: 353 / 451
## Progress: 354 / 451
## Progress: 355 / 451
## Progress: 356 / 451
## Progress: 357 / 451
## Progress: 358 / 451
## Progress: 359 / 451
## Progress: 360 / 451
## Progress: 361 / 451
## Progress: 362 / 451
## Progress: 363 / 451
## Progress: 364 / 451
## Progress: 365 / 451
## Progress: 366 / 451
## Progress: 367 / 451
## Progress: 368 / 451
## Progress: 369 / 451
## Progress: 370 / 451
## Progress: 371 / 451
## Progress: 372 / 451
## Progress: 373 / 451
## Progress: 374 / 451
## Progress: 375 / 451
## Progress: 376 / 451
## Progress: 377 / 451
## Progress: 378 / 451
## Progress: 379 / 451
## Progress: 380 / 451
## Progress: 381 / 451
## Progress: 382 / 451
## Progress: 383 / 451
## Progress: 384 / 451
## Progress: 385 / 451
## Progress: 386 / 451
## Progress: 387 / 451
## Progress: 388 / 451
## Progress: 389 / 451
## Progress: 390 / 451
## Progress: 391 / 451
## Progress: 392 / 451
## Progress: 393 / 451
## Progress: 394 / 451
## Progress: 395 / 451
## Progress: 396 / 451
## Progress: 397 / 451
## Progress: 398 / 451
## Progress: 399 / 451
## Progress: 400 / 451
## Progress: 401 / 451
## Progress: 402 / 451
## Progress: 403 / 451
## Progress: 404 / 451
## Progress: 405 / 451
## Progress: 406 / 451
## Progress: 407 / 451
## Progress: 408 / 451
## Progress: 409 / 451
## Progress: 410 / 451
## Progress: 411 / 451
## Progress: 412 / 451
## Progress: 413 / 451
## Progress: 414 / 451
## Progress: 415 / 451
## Progress: 416 / 451
## Progress: 417 / 451
## Progress: 418 / 451
## Progress: 419 / 451
## Progress: 420 / 451
## Progress: 421 / 451
## Progress: 422 / 451
## Progress: 423 / 451
## Progress: 424 / 451
## Progress: 425 / 451
## Progress: 426 / 451
## Progress: 427 / 451
## Progress: 428 / 451
## Progress: 429 / 451
## Progress: 430 / 451
## Progress: 431 / 451
## Progress: 432 / 451
## Progress: 433 / 451
## Progress: 434 / 451
## Progress: 435 / 451
## Progress: 436 / 451
## Progress: 437 / 451
## Progress: 438 / 451
## Progress: 439 / 451
## Progress: 440 / 451
## Progress: 441 / 451
## Progress: 442 / 451
## Progress: 443 / 451
## Progress: 444 / 451
## Progress: 445 / 451
## Progress: 446 / 451
## Progress: 447 / 451
## Progress: 448 / 451
## Progress: 449 / 451
## Progress: 450 / 451
## Progress: 451 / 451plotCutoff(Output, 0.05)
Cpow <- getPower(Output)
findPower(Cpow, "N", 0.80)## f1=~y1 f1=~y2 f1=~y3 f2=~y4 f2=~y5 f2=~y6 f3=~y7 f3=~y8 f3=~y9 y1~~y1
## Inf Inf 52 Inf Inf Inf 67 59 59 50
## y2~~y2 y3~~y3 y4~~y4 y5~~y5 y6~~y6 y7~~y7 y8~~y8 y9~~y9 f3~f1 f3~f2
## 50 50 50 50 50 50 50 50 71 93
## y1~1 y2~1 y3~1 y4~1 y5~1 y6~1 y7~1 y8~1 y9~1
## NA NA NA NA NA NA NA NA NA