Skip to content

Commit

Permalink
version 0.9.2
Browse files Browse the repository at this point in the history
  • Loading branch information
@johnnyzhz @cran-robot
johnnyzhz authored and cran-robot committed May 10, 2023
1 parent 8096ed5 commit 654077f
Show file tree
Hide file tree
Showing 12 changed files with 247 additions and 250 deletions.
8 changes: 4 additions & 4 deletions DESCRIPTION
View file
@@ -1,7 +1,7 @@
Package: WebPower
Title: Basic and Advanced Statistical Power Analysis
Version: 0.9.0
Date: 2023-04-16
Version: 0.9.2
Date: 2023-05-09
Authors@R: c(person("Zhiyong", "Zhang", role = c("aut", "cre"),
email = "johnnyzhz@gmail.com"),
person("Yujiao", "Mai", role = "aut"),
Expand All @@ -25,6 +25,6 @@ Encoding: UTF-8
LazyLoad: yes
LazyData: yes
NeedsCompilation: no
Packaged: 2023-04-20 17:05:43 UTC; zzhang4
Packaged: 2023-05-10 01:01:20 UTC; zzhang4
Repository: CRAN
Date/Publication: 2023-04-20 18:02:37 UTC
Date/Publication: 2023-05-10 02:30:02 UTC
22 changes: 11 additions & 11 deletions MD5
View file
@@ -1,10 +1,10 @@
c789f5404bdf5b6774c5fde70fe240b2 *DESCRIPTION
268d8f29a27a853530f1423bdad590fd *DESCRIPTION
7d64335ba38d220709a921137814ee14 *NAMESPACE
a0bcec493c022a6e90c045aff311c3a6 *R/model14.R
e27c5aff1b059e0363768f2d08fa7e1e *R/model15.R
5208bbf4ba29a6d91130aa89f2a4060c *R/model58.R
bf344a714d712897d11ca09860ab631a *R/model7.R
df3ac23988eb7f24ce567e8c9f2ea0c6 *R/model8.R
950a9b89b49627f0f4d37baf4b901a1b *R/model14.R
11d33652db6615837b9b682f2ab7a2da *R/model15.R
ee4d00e15593b6ec2b9f583171095aef *R/model58.R
cf095a79c2a8266843b3d147df174096 *R/model7.R
19cd026ef07667091f1bd61ffc17b15a *R/model8.R
7e04ca484da8d6ceb462867d83b1ada8 *R/webpower.R
83b0027813eefbbbb56cdfae00a9adcf *data/webpower.RData
f4c53917e7a6ce6720539abe8c6f18a1 *inst/CRT2.txt
Expand Down Expand Up @@ -44,11 +44,11 @@ f93dcccad006d428281f7e0a359bacb0 *man/wp.mc.sem.power.curve.Rd
e3aa54355131249a67530e681459f10a *man/wp.mc.t.Rd
60375052e10b52d9268e0edb729742e2 *man/wp.mediation.Rd
b82ac428a511dcf6b45e9af35422e00b *man/wp.mmrm.Rd
dda9ad6387dcd7282f099fc61ff2e8f9 *man/wp.modmed.m14.Rd
997b1223aa006f8f095ca550cef30a6d *man/wp.modmed.m15.Rd
7bb3072286b579188a32716c571c1617 *man/wp.modmed.m58.Rd
0a91119994e601cb5b8dde2c92595d4f *man/wp.modmed.m7.Rd
71cb09b42c462ecb0c5ddde411530cc1 *man/wp.modmed.m8.Rd
c97d1a521efad8bfe8b5a2f8a4c9b598 *man/wp.modmed.m14.Rd
6a70dca0f9a5479a08ea3053943961cc *man/wp.modmed.m15.Rd
6014bbfc4c4a41fcc2f617a0681cb093 *man/wp.modmed.m58.Rd
241b0ad20c9c7ba24f9ad71e8ca0880a *man/wp.modmed.m7.Rd
cb4817b654955c55b6c5aba5a4c72886 *man/wp.modmed.m8.Rd
25eb6c557b4a4accd999ef15bc8cac40 *man/wp.mrt2arm.Rd
9015883be37879a88e9585d072e03de6 *man/wp.mrt3arm.Rd
cf30c06a5737b7fa17b8f3f65fd7a6f0 *man/wp.poisson.Rd
Expand Down
45 changes: 22 additions & 23 deletions R/model14.R
View file
Expand Up @@ -5,10 +5,10 @@
#' power analysis of model 14 in Introduction to Mediation, Moderation, and Conditional Process Analysis
#'
#' @param a1 regression coefficient of mediator (m) on predictor (x)
#' @param cc regression coefficient of outcome (y) on predictor (x)
#' @param cp regression coefficient of outcome (y) on predictor (x)
#' @param b1 regression coefficient of outcome (y) on mediator (m)
#' @param c1 regression coefficient of outcome (y) on moderator (w)
#' @param c2 regression coefficient of outcome (y) on the product (mw)
#' @param d1 regression coefficient of outcome (y) on moderator (w)
#' @param b2 regression coefficient of outcome (y) on the product (mw)
#' @param sigx2 variance of predictor (x)
#' @param sigw2 variance of moderator (w)
#' @param sige12 variance of error in the first regression equation
Expand All @@ -28,13 +28,13 @@
#' @return power of indirect effect, direct effect, and moderation
#' @export
#' @examples
#' test = wp.modmed.m14(a1 = 0.2, cc = 0.2, b1 = 0.5, c1 = 0.5, c2 = 0.2, sigx2 = 1,
#' test = wp.modmed.m14(a1 = 0.2, cp = 0.2, b1 = 0.5, d1 = 0.5, b2 = 0.2, sigx2 = 1,
#' sigw2 = 1, sige12 = 1, sige22 = 1, sigx_w = 0.5, n = 50,
#' nrep = 100, alpha = 0.05, b = 1000, ncore = 1)
#' print(test)
wp.modmed.m14 <- function (a1 = 0.2, cc = 0.2, b1 = 0.5, c1 = 0.5, c2 = 0.2,
wp.modmed.m14 <- function (a1 = 0.2, cp = 0.2, b1 = 0.5, d1 = 0.5, b2 = 0.2,
sigx2 = 1, sigw2 = 1, sige12 = 1, sige22 = 1,
sigx_w = 0.5, n = 100, nrep = 1000, alpha = 0.05,
sigx_w = 0.5, n = 100, nrep = 100, alpha = 0.05,
b = 1000, nb = n, w_value = 0, method = "value",
ncore = 1, pop.cov = NULL, mu = NULL,
varnames = c('y', 'x', 'w', 'm', 'mw'))
Expand All @@ -43,15 +43,15 @@ wp.modmed.m14 <- function (a1 = 0.2, cc = 0.2, b1 = 0.5, c1 = 0.5, c2 = 0.2,
if (is.null(pop.cov) || is.null(mu)) {
sigm2 = a1^2*sigx2 + sige12
sigxw2 = sigx2*sigw2 + sigx_w^2
sigy2 = (cc + a1*b1)^2*sigx2 + c1^2*sigw2 + c2^2*sige12 *
sigw2 + (a1*c2)^2*sigxw2 + b1^2*sige12 + sige22 + 2*(cc + a1*b1)*c1*sigx_w
sigy2 = (cp + a1*b1)^2*sigx2 + d1^2*sigw2 + b2^2*sige12 *
sigw2 + (a1*b2)^2*sigxw2 + b1^2*sige12 + sige22 + 2*(cp + a1*b1)*d1*sigx_w
sigmw2 = a1^2*sigxw2 + sige12*sigw2
sigx_m = a1*sigx2
sigx_y = (cc + a1*b1)*sigx2 + c1*sigx_w
sigm_y = a1*(cc + a1*b1)*sigx2 + a1*c1*sigx_w + b1*sige12
sigx_y = (cp + a1*b1)*sigx2 + d1*sigx_w
sigm_y = a1*(cp + a1*b1)*sigx2 + a1*d1*sigx_w + b1*sige12
sigm_w = a1*sigx_w
sigy_w = (cc + a1*b1)*sigx_w + c1*sigw2
sigy_mw = a1^2*c2*sigxw2 + c2*sige12*sigw2
sigy_w = (cp + a1*b1)*sigx_w + d1*sigw2
sigy_mw = a1^2*b2*sigxw2 + b2*sige12*sigw2

# y, x, w, m, mw
pop.cov = array(
Expand All @@ -64,7 +64,7 @@ wp.modmed.m14 <- function (a1 = 0.2, cc = 0.2, b1 = 0.5, c1 = 0.5, c2 = 0.2,
)

u_mw = a1*sigx_w
u_y = a1*c2*sigx_w
u_y = a1*b2*sigx_w
colnames(pop.cov) = rownames(pop.cov) = c('y', 'x', 'w', 'm', 'mw')
mu = c(u_y, 0, 0, 0, u_mw)
}else{
Expand All @@ -86,21 +86,21 @@ wp.modmed.m14 <- function (a1 = 0.2, cc = 0.2, b1 = 0.5, c1 = 0.5, c2 = 0.2,
test_boot2 = lm(y ~ x + m + w + mw, data = boot_data)
boot_CI = test_boot1$coefficients[2]*(test_boot2$coefficients[3] + test_boot2$coefficients[5]*w_value)
boot_CD = test_boot2$coefficients[2]
boot_c2 = as.numeric(test_boot2$coefficients[5])
boot_b2 = as.numeric(test_boot2$coefficients[5])
boot_CI1 = test_boot1$coefficients[2]
boot_CI2 = (test_boot2$coefficients[3] + test_boot2$coefficients[5]*w_value)
return(list(boot_CI, boot_CD, boot_c2, boot_CI1, boot_CI2))
return(list(boot_CI, boot_CD, boot_b2, boot_CI1, boot_CI2))
}
boot_effect = lapply(1:b, bootstrap)
boot_CI = matrix(0, ncol = 1, nrow = b)
boot_CD = matrix(0, ncol = 1, nrow = b)
boot_c2 = matrix(0, ncol = 1, nrow = b)
boot_b2 = matrix(0, ncol = 1, nrow = b)
boot_CI1 = matrix(0, ncol = 1, nrow = b)
boot_CI2 = matrix(0, ncol = 1, nrow = b)

boot_CI = t(sapply(1:b, function(i) unlist(boot_effect[[i]][1])))
boot_CD = t(sapply(1:b, function(i) unlist(boot_effect[[i]][2])))
boot_c2 = t(sapply(1:b, function(i) unlist(boot_effect[[i]][3])))
boot_b2 = t(sapply(1:b, function(i) unlist(boot_effect[[i]][3])))
boot_CI1 = t(sapply(1:b, function(i) unlist(boot_effect[[i]][4])))
boot_CI2 = t(sapply(1:b, function(i) unlist(boot_effect[[i]][5])))

Expand All @@ -109,7 +109,7 @@ wp.modmed.m14 <- function (a1 = 0.2, cc = 0.2, b1 = 0.5, c1 = 0.5, c2 = 0.2,
interval_CI1 = matrix(0, ncol = 1, nrow = 2)
interval_CI2 = matrix(0, ncol = 1, nrow = 2)
interval_CD = matrix(0, ncol = 1, nrow = 2)
interval_c2 = matrix(0, ncol = 1, nrow = 2)
interval_b2 = matrix(0, ncol = 1, nrow = 2)

interval_CI[, 1] = quantile(boot_CI,
probs = c(alpha / 2, 1 - alpha / 2),
Expand All @@ -123,18 +123,18 @@ wp.modmed.m14 <- function (a1 = 0.2, cc = 0.2, b1 = 0.5, c1 = 0.5, c2 = 0.2,
interval_CD[, 1] = quantile(boot_CD,
probs = c(alpha / 2, 1 - alpha / 2),
names = T)
interval_c2[, 1] = quantile(boot_c2,
interval_b2[, 1] = quantile(boot_b2,
probs = c(alpha / 2, 1 - alpha / 2),
names = T)


r_CI = as.numeric(!sapply(1, function(i) dplyr::between(0, interval_CI[1, i], interval_CI[2, i])))
r_CD = as.numeric(!sapply(1, function(i) dplyr::between(0, interval_CD[1, i], interval_CD[2, i])))
r_c2 = as.numeric(!sapply(1, function(i) dplyr::between(0, interval_c2[1, i], interval_c2[2, i])))
r_b2 = as.numeric(!sapply(1, function(i) dplyr::between(0, interval_b2[1, i], interval_b2[2, i])))
if (method == "joint") {
r_CI = as.numeric(!dplyr::between(0, interval_CI1[1, 1], interval_CI1[2, 1]))*as.numeric(!dplyr::between(0, interval_CI2[1, 1], interval_CI2[2, 1]))
}
power = c(r_CI, r_CD, r_c2)
power = c(r_CI, r_CD, r_b2)
return(power)
}

Expand All @@ -143,7 +143,7 @@ wp.modmed.m14 <- function (a1 = 0.2, cc = 0.2, b1 = 0.5, c1 = 0.5, c2 = 0.2,

parallel::clusterExport(
CL1,
c('a1', 'cc', 'b1', 'c1', 'c2', 'sigx2', 'sigw2',
c('a1', 'cp', 'b1', 'd1', 'b2', 'sigx2', 'sigw2',
'sige12', 'sige22', 'sigx_w', 'n', 'nrep',
'alpha', 'b','nb', 'pop.cov', 'mu', 'method'
),
Expand Down Expand Up @@ -182,4 +182,3 @@ power3 is the power of moderation on the path m to y."
}



47 changes: 24 additions & 23 deletions R/model15.R
View file
Expand Up @@ -4,11 +4,11 @@
#' power analysis of model 15 in Introduction to Mediation, Moderation, and Conditional Process Analysis
#'
#' @param a1 regression coefficient of mediator (m) on predictor (x)
#' @param c regression coefficient of outcome (y) on predictor (x)
#' @param cp regression coefficient of outcome (y) on predictor (x)
#' @param b1 regression coefficient of outcome (y) on mediator (m)
#' @param b2 regression coefficient of outcome (y) on the product (mw)
#' @param c1 regression coefficient of outcome (y) on moderator (w)
#' @param c2 regression coefficient of outcome (y) on the product (xw)
#' @param d1 regression coefficient of outcome (y) on moderator (w)
#' @param d2 regression coefficient of outcome (y) on the product (xw)
#' @param sigx2 variance of predictor (x)
#' @param sigw2 variance of moderator (w)
#' @param sige12 variance of error in the first regression equation
Expand All @@ -28,32 +28,32 @@
#' @return power of indirect effect, direct effect, and moderation
#' @export
#' @examples
#' test = wp.modmed.m15(a1 = 0.6, c = 0.2, b1 = 0.3, b2 = 0.2, c1 = 0.2, c2 = 0.1,
#' test = wp.modmed.m15(a1 = 0.6, cp = 0.2, b1 = 0.3, b2 = 0.2, d1 = 0.2, d2 = 0.1,
#' sigx2 = 1, sigw2 = 1, sige12 = 1, sige22 = 1, sigx_w = 0.4,
#' n = 50, nrep = 100, alpha = 0.05, b = 1000, ncore = 1)
#'print(test)
wp.modmed.m15 <- function(a1 = 0.6, c = 0.2, b1 = 0.3, b2 = 0.5, c1 = 0.2, c2 = 0.1,
wp.modmed.m15 <- function(a1 = 0.6, cp = 0.2, b1 = 0.3, b2 = 0.5, d1 = 0.2, d2 = 0.1,
sigx2 = 1, sigw2 = 1, sige12 = 1, sige22 = 1, sigx_w = 0.4,
n = 500, nrep = 1000, alpha = 0.05, b = 1000, nb = n,
n = 100, nrep = 100, alpha = 0.05, b = 1000, nb = n,
w_value = 0, method = "value", ncore = 1,
pop.cov = NULL, mu = NULL,
varnames = c('y','x','w','m','xw','mw'))
{
if (is.null(pop.cov) || is.null(mu)){
sigxw2 = sigx2*sigw2 + sigx_w^2
sigm_xw = 0
sigy_xw = (c2 + a1*b2)*sigxw2
sigy_xw = (d2 + a1*b2)*sigxw2
sigm_x = a1*sigx2
sigm_w = a1*sigx_w
sigy_x = c1*sigx_w + (c + a1*b1)*sigx2
sigy_w = c1*sigw2 + (c + a1*b1)*sigx_w
sigy_x = d1*sigx_w + (cp + a1*b1)*sigx2
sigy_w = d1*sigw2 + (cp + a1*b1)*sigx_w
sige1w2 = sige12*sigw2
sigy2 = (c + a1*b1)^2*sigx2 + (c2 + a1*b2)^2*sigxw2 + b2 ^
2*sige1w2 + c1^2*sigw2 + 2*c1*(c + a1*b1)*sigx_w + b1^2 *
sigy2 = (cp + a1*b1)^2*sigx2 + (d2 + a1*b2)^2*sigxw2 + b2 ^
2*sige1w2 + d1^2*sigw2 + 2*d1*(cp + a1*b1)*sigx_w + b1^2 *
sige12 + sige22
sigm2 = a1^2*sigx2 + sige12
sigy_m = a1*(c + a1*b1)*sigx2 + a1*c1*sigx_w + b1*sige12
sigy_mw = a1*(c2 + a1*b2)*sigxw2 + b2*sige1w2
sigy_m = a1*(cp + a1*b1)*sigx2 + a1*d1*sigx_w + b1*sige12
sigy_mw = a1*(d2 + a1*b2)*sigxw2 + b2*sige1w2
sigxw_mw = a1*sigxw2
sigmw2 = a1^2*sigxw2 + sige1w2

Expand All @@ -71,7 +71,7 @@ wp.modmed.m15 <- function(a1 = 0.6, c = 0.2, b1 = 0.3, b2 = 0.5, c1 = 0.2, c2 =
rownames(pop.cov) = colnames(pop.cov) = c('y', 'x', 'w', 'm', 'xw', 'mw')
u_xw = sigx_w
u_mw = sigm_w
u_y = c2*u_xw + b2*u_mw
u_y = d2*u_xw + b2*u_mw
mu = c(u_y, 0, 0, 0, u_xw, u_mw)
}else{
pop.cov = pop.cov
Expand All @@ -93,23 +93,23 @@ wp.modmed.m15 <- function(a1 = 0.6, c = 0.2, b1 = 0.3, b2 = 0.5, c1 = 0.2, c2 =
test_boot2 = lm(y ~ x + m + w + xw + mw, data = boot_data)
boot_CI = (test_boot2$coefficients[3] + test_boot2$coefficients[6]*w_value)*test_boot1$coefficients[2]
boot_CD = test_boot2$coefficients[2] + test_boot2$coefficients[5]*w_value
boot_c2 = test_boot2$coefficients[5]
boot_d2 = test_boot2$coefficients[5]
boot_b2 = test_boot2$coefficients[6]
boot_CI1 = test_boot1$coefficients[2]
boot_CI2 = (test_boot2$coefficients[3] + test_boot2$coefficients[6]*w_value)
return(list(boot_CI, boot_CD, boot_c2, boot_b2, boot_CI1, boot_CI2))
return(list(boot_CI, boot_CD, boot_d2, boot_b2, boot_CI1, boot_CI2))
}
boot_effect = lapply(1:b, bootstrap)
boot_CI = matrix(0, ncol = 1, nrow = b)
boot_CI1 = matrix(0, ncol = 1, nrow = b)
boot_CI2 = matrix(0, ncol = 1, nrow = b)
boot_CD = matrix(0, ncol = 1, nrow = b)
boot_c2 = matrix(0, ncol = 1, nrow = b)
boot_d2 = matrix(0, ncol = 1, nrow = b)
boot_b2 = matrix(0, ncol = 1, nrow = b)

boot_CI=t(sapply(1:b,function(i) unlist(boot_effect[[i]][1])))
boot_CD=t(sapply(1:b,function(i) unlist(boot_effect[[i]][2])))
boot_c2=t(sapply(1:b,function(i) unlist(boot_effect[[i]][3])))
boot_d2=t(sapply(1:b,function(i) unlist(boot_effect[[i]][3])))
boot_b2=t(sapply(1:b,function(i) unlist(boot_effect[[i]][4])))
boot_CI1=t(sapply(1:b,function(i) unlist(boot_effect[[i]][5])))
boot_CI2=t(sapply(1:b,function(i) unlist(boot_effect[[i]][6])))
Expand All @@ -118,7 +118,7 @@ wp.modmed.m15 <- function(a1 = 0.6, c = 0.2, b1 = 0.3, b2 = 0.5, c1 = 0.2, c2 =
interval_CI1=matrix(0,ncol=1,nrow=2)
interval_CI2=matrix(0,ncol=1,nrow=2)
interval_CD=matrix(0,ncol=1,nrow=2)
interval_c2=matrix(0,ncol=1,nrow=2)
interval_d2=matrix(0,ncol=1,nrow=2)
interval_b2=matrix(0,ncol=1,nrow=2)

interval_CI[, 1] = quantile(boot_CI,
Expand All @@ -133,7 +133,7 @@ wp.modmed.m15 <- function(a1 = 0.6, c = 0.2, b1 = 0.3, b2 = 0.5, c1 = 0.2, c2 =
interval_CD[, 1] = quantile(boot_CD,
probs = c(alpha / 2, 1 - alpha / 2),
names = T)
interval_c2[, 1] = quantile(boot_c2,
interval_d2[, 1] = quantile(boot_d2,
probs = c(alpha / 2, 1 - alpha / 2),
names = T)
interval_b2[, 1] = quantile(boot_b2,
Expand All @@ -143,17 +143,17 @@ wp.modmed.m15 <- function(a1 = 0.6, c = 0.2, b1 = 0.3, b2 = 0.5, c1 = 0.2, c2 =

r_CI=as.numeric(!sapply(1,function(i) dplyr::between(0,interval_CI[1,i],interval_CI[2,i])))
r_DI=as.numeric(!sapply(1,function(i) dplyr::between(0,interval_CD[1,i],interval_CD[2,i])))
r_c2=as.numeric(!sapply(1:1,function(i) dplyr::between(0,interval_c2[1,i],interval_c2[2,i])))
r_d2=as.numeric(!sapply(1:1,function(i) dplyr::between(0,interval_d2[1,i],interval_d2[2,i])))
r_b2=as.numeric(!sapply(1:1,function(i) dplyr::between(0,interval_b2[1,i],interval_b2[2,i])))
if (method =="joint"){
r_CI=as.numeric(!dplyr::between(0,interval_CI1[1,1],interval_CI1[2,1]))*as.numeric(!dplyr::between(0,interval_CI2[1,1],interval_CI2[2,1]))
}
power=c(r_CI,r_DI, r_c2, r_b2)
power=c(r_CI,r_DI, r_d2, r_b2)
return(power)
}
if (ncore > 1){
CL1=parallel::makeCluster(ncore)
parallel::clusterExport(CL1,c('a1','c','b1','b2','c1','c2',
parallel::clusterExport(CL1,c('a1','c','b1','b2','d1','d2',
'sigx2','sigw2','sige12','sige22','sigx_w',
'n','nrep','alpha','b','nb','pop.cov','u_y','u_xw',
'u_mw', 'mu', 'method', 'w_value'),envir = environment())
Expand Down Expand Up @@ -187,3 +187,4 @@ power4 is the power of moderation on the path m to y."), class = "webpower")
}



0 comments on commit 654077f

Please sign in to comment.