R/pcurve2.R
In MathiasHarrer/dmetar: Companion R Package for the Guide 'Doing Meta-Analysis in R'
Defines functions
pcurve
Documented in
pcurve
#' Perform a \emph{p}-curve analysis
#'
#' This function performs a \eqn{p}-curve analysis using a \code{meta} object or calculated effect size data.
#'
#' @usage pcurve(x, effect.estimation = FALSE, N, dmin = 0, dmax = 1)
#'
#' @param x Either an object of class \code{meta}, generated by the \code{metagen}, \code{metacont},
#' \code{metacor}, \code{metainc}, or \code{metabin} function, or a dataframe containing the calculated effect size
#' (named \code{TE}, log-transformed if based on a ratio), standard error (named \code{seTE}) and study label (named \code{studlab})
#' for each study.
#' @param effect.estimation Logical. Should the true effect size underlying the \emph{p}-curve be estimated?
#' If set to \code{TRUE}, a vector containing the total sample size for each study must be provided for
#' \code{N}. \code{FALSE} by default.
#' @param N A numeric vector of same length as the number of effect sizes included in \code{x} specifying the
#' total sample size \eqn{N} corresponding to each effect. Only needed if \code{effect.estimation = TRUE}.
#' @param dmin If \code{effect.estimation = TRUE}: lower limit for the effect size (\eqn{d}) space in which
#' the true effect size should be searched. Must be greater or equal to 0. Default is 0.
#' @param dmax If \code{effect.estimation = TRUE}: upper limit for the effect size (\eqn{d}) space in which
#' the true effect size should be searched. Must be greater than 0. Default is 1.
#'
#' @details
#' \strong{P-curve Analysis}
#'
#' \eqn{P}-curve analysis (Simonsohn, Simmons & Nelson, 2014, 2015) has been proposed as a method
#' to detect \eqn{p}-hacking and publication bias in meta-analyses.
#'
#' \eqn{P}-Curve assumes that publication bias
#' is not only generated because researchers do not publish non-significant results,
#' but also because analysts “play” around with their data ("\eqn{p}-hacking"; e.g., selectively removing outliers,
#' choosing different outcomes, controlling for different variables) until a non-significant
#' finding becomes significant (i.e., \eqn{p<0.05}).
#'
#' The method assumes that for a specific research
#' question, \eqn{p}-values smaller 0.05 of included studies should follow a right-skewed distribution
#' if a true effect exists, even when the power in single studies was (relatively) low. Conversely,
#' a left-skewed \eqn{p}-value distribution indicates the presence of \eqn{p}-hacking and absence of
#' a true underlying effect. To control for "ambitious" \eqn{p}-hacking, \eqn{P}-curve also incorporates a
#' "half-curve" test (Simonsohn, Simmons & Nelson, 2014, 2015).
#'
#' Simonsohn et al. (2014)
#' stress that \eqn{p}-curve analysis should only be used for test statistics which were actually of interest
#' in the context of the included study, and that a detailed table documenting the reported results
#' used in for the \eqn{p}-curve analysis should be created before communicating
#' results (\href{http://www.p-curve.com/Supplement/}{link}).
#'
#' \strong{Implementation in the function}
#'
#' To generate the \eqn{p}-curve and conduct the analysis, this function reuses parts of the \emph{R} code underlying
#' the \href{http://p-curve.com/app4/pcurve_app4.052.r}{P-curve App 4.052} (Simonsohn, 2017). The effect sizes
#' included in the \code{meta} object or \code{data.frame} provided for \code{x} are transformed
#' into \eqn{z}-values internally, which are then used to calculate {p}-values and conduct the
#' Stouffer and Binomial test used for the \eqn{p}-curve analysis. Interpretations of the function
#' concerning the presence or absence/inadequateness of evidential value are made according to the
#' guidelines described by Simonsohn, Simmons and Nelson (2015):
#'
#' \itemize{
#' \item \strong{Evidential value present}: The right-skewness test is significant for the half curve with
#' \eqn{p<0.05} \strong{or} the \eqn{p}-value of the right-skewness test is \eqn{<0.1} for both the half and full curve.
#' \item \strong{Evidential value absent or inadequate}: The flatness test is \eqn{p<0.05} for the full curve
#' \strong{or} the flatness test for the half curve and the binomial test are \eqn{p<0.1}.
#'}
#'
#' For effect size estimation, the \code{pcurve} function implements parts of the loss function
#' presented in Simonsohn, Simmons and Nelson (2014b).
#' The function generates a loss function for candidate effect sizes \eqn{\hat{d}}, using \eqn{D}-values in
#' a Kolmogorov-Smirnov test as the metric of fit, and the value of \eqn{\hat{d}} which minimizes \eqn{D}
#' as the estimated true effect.
#'
#' It is of note that a lack of robustness of \eqn{p}-curve analysis results
#' has been noted for meta-analyses with substantial heterogeneity (van Aert, Wicherts, & van Assen, 2016).
#' Following van Aert et al., adjusted effect size estimates should only be
#' reported and interpreted for analyses with \eqn{I^2} values below 50 percent.
#' A warning message is therefore printed by
#' the \code{pcurve} function when \code{x} is of class \code{meta} and the between-study heterogeneity
#' of the meta-analysis is substantial (i.e., \eqn{I^2} greater than 50 percent).
#'
#'
#' @references Harrer, M., Cuijpers, P., Furukawa, T.A, & Ebert, D. D. (2019).
#' \emph{Doing Meta-Analysis in R: A Hands-on Guide}. DOI: 10.5281/zenodo.2551803.
#' \href{https://bookdown.org/MathiasHarrer/Doing_Meta_Analysis_in_R/pcurve.html}{Chapter 9.2}.
#'
#' Simonsohn, U., Nelson, L. D., & Simmons, J. P. (2014a). P-curve: a Key to the File-drawer.
#' \emph{Journal of Experimental Psychology, 143}(2), 534.
#'
#' Simonsohn, U., Nelson, L. D. & Simmons, J. P. (2014b). P-Curve and Effect Size:
#' Correcting for Publication Bias Using Only Significant Results.
#' \emph{Perspectives on Psychological Science 9}(6), 666–81.
#'
#' Simonsohn, U., Nelson, L. D. & Simmons, J. P. (2015). Better P-Curves: Making P-Curve
#' Analysis More Robust to Errors, Fraud, and Ambitious P-Hacking, a Reply to Ulrich and Miller (2015).
#' \emph{Journal of Experimental Psychology, 144}(6), 1146-1152.
#'
#' Simonsohn, U. (2017). R code for the P-Curve App 4.052. http://p-curve.com/app4/pcurve_app4.052.r (Accessed 2019-08-16).
#'
#' Van Aert, R. C., Wicherts, J. M., & van Assen, M. A. (2016).
#' Conducting meta-analyses based on p values: Reservations and recommendations for applying
#' \emph{p}-uniform and \emph{p}-curve. \emph{Perspectives on Psychological Science, 11}(5), 713-729.
#'
#' @author Mathias Harrer & David Daniel Ebert
#'
#' @return Returns a plot and main results of the pcurve analysis:
#' \itemize{
#' \item \strong{P-curve plot}: A plot displaying the observed \eqn{p}-curve and significance results
#' for the right-skewness and flatness test.
#' \item \strong{Number of studies}: The number of studies provided for the analysis, the number
#' of significant \eqn{p}-values included in the analysis, and the number of studies with \eqn{p<0.025}
#' used for the half-curve tests.
#' \item \strong{Test results}: The results for the right-skewness and flatness test, including the
#' \eqn{p_{binomial}} value, as well as the \eqn{z} and \eqn{p} value for the full and half-curve test.
#' \item \strong{Power Estimate}: The power estimate and 95\% confidence interval.
#' \item \strong{Evidential value}: Two lines displaying if evidential value is present and/or absent/inadequate based
#' on the results (using the guidelines by Simonsohn et al., 2015, see Details).
#' \item \strong{True effect estimate}: If \code{effect.estimation} is set to \code{TRUE}, the estimated true effect
#' \eqn{\hat{d}} is returned additionally.
#'}
#'
#'
#' If results are saved to a variable, a list of class \code{pcurve} containing the following objects is returned:
#' \itemize{
#' \item \code{pcurveResults}: A data frame containing the results for the right-skewness and flatness test, including the
#' \eqn{p_{binomial}} value, as well as the \eqn{z} and \eqn{p} value for the full and half-curve test.
#' \item \code{Power}: The power estimate and 95\% confidence interval.
#' \item \code{PlotData}: A data frame with the data used in the \eqn{p}-curve plot.
#' \item \code{Input}: A data frame containing the provided effect sizes, calculated \eqn{p}-values and individual results for each included (significant) effect.
#' \item \code{EvidencePresent}, \code{EvidenceAbsent}, \code{kInput}, \code{kAnalyzed}, \code{kp0.25}: Further results of the \eqn{p}-curve analysis, including the presence/absence of evidence interpretation,
#' and number of provided/significant/\eqn{p<0.025} studies.
#' \item \code{I2}: \eqn{I^2}-Heterogeneity of the studies provided as input (only when \code{x} is of class \code{meta}).
#' \item \code{class.meta.object}: \code{class} of the original object provided in \code{x}.
#'}
#'
#'
#' @import stringr poibin
#' @importFrom graphics abline axis lines mtext par plot points rect segments text
#' @importFrom stats as.formula hat influence ks.test optimize pbinom pchisq pf pnorm pt punif qchisq qf qnorm qt reformulate reorder setNames uniroot
#'
#' @export pcurve
#'
#' @seealso
#' \code{\link{eggers.test}}
#'
#' @examples
#' # Example 1: Use metagen object, do not estimate d
#' suppressPackageStartupMessages(library(meta))
#' data("ThirdWave")
#' meta1 = metagen(TE,seTE, studlab=ThirdWave$Author, data=ThirdWave)
#' pcurve(meta1)
#'
#' # Example 2: Provide Ns, calculate d estimate
#' N = c(105, 161, 60, 37, 141, 82, 97, 61, 200, 79, 124, 25, 166, 59, 201, 95, 166, 144)
#' pcurve(meta1, effect.estimation = TRUE, N = N)
#'
#' # Example 3: Use metacont object, calculate d estimate
#' data("amlodipine")
#' meta2 <- metacont(n.amlo, mean.amlo, sqrt(var.amlo),
#' n.plac, mean.plac, sqrt(var.plac),
#' data=amlodipine, studlab=study, sm="SMD")
#' N = amlodipine$n.amlo + amlodipine$n.plac
#' pcurve(meta2, effect.estimation = TRUE, N = N, dmin = 0, dmax = 1)
#'
#' # Example 4: Construct x object from scratch
#' sim = data.frame("studlab" = c(paste("Study_", 1:18, sep = "")),
#' "TE" = c(0.561, 0.296, 0.648, 0.362, 0.770, 0.214, 0.476,
#' 0.459, 0.343, 0.804, 0.357, 0.476, 0.638, 0.396, 0.497,
#' 0.384, 0.568, 0.415),
#' "seTE" = c(0.338, 0.297, 0.264, 0.258, 0.279, 0.347, 0.271, 0.319,
#' 0.232, 0.237, 0.385, 0.398, 0.342, 0.351, 0.296, 0.325,
#' 0.322, 0.225))
#' pcurve(sim)
pcurve = function(x, effect.estimation = FALSE, N, dmin = 0, dmax = 1){
# Rename x to metaobject, remove x
metaobject = x
rm(x)
# Stop if metaobject is not meta or does not contain TE or seTE column
if (!(class(metaobject)[1] %in% c("metagen", "metabin", "metacont", "metacor", "metainc", "meta", "metaprop"))){
for (i in 1:length(colnames(metaobject))){
te.exists = FALSE
if (colnames(metaobject)[i]=="TE"){
te.exists = TRUE
break
} else {
}
}
for (i in 1:length(colnames(metaobject))){
sete.exists = FALSE
if (colnames(metaobject)[i]=="seTE"){
sete.exists = TRUE
break
} else {
}
}
for (i in 1:length(colnames(metaobject))){
studlab.exists = FALSE
if (colnames(metaobject)[i]=="studlab"){
studlab.exists = TRUE
break
} else {
}
}
if(te.exists == FALSE | sete.exists ==FALSE | studlab.exists ==FALSE){
stop("x must be a meta-analysis object generated by meta functions or a data.frame with columns labeled studlab, TE, and seTE.")
}
}
#Disable scientific notation
options(scipen=999)
# Calculate Z
zvalues.input = abs(metaobject$TE/metaobject$seTE)
##############################################
# 1. Functions ###############################
##############################################
getncp.f =function(df1,df2, power) {
error = function(ncp_est, power, x, df1,df2) pf(x, df1 = df1, df2=df2, ncp = ncp_est) - (1-power)
xc=qf(p=.95, df1=df1,df2=df2)
return(uniroot(error, c(0, 1000), x = xc, df1 = df1,df2=df2, power=power)$root) }
getncp.c =function(df, power) {
xc=qchisq(p=.95, df=df)
error = function(ncp_est, power, x, df) pchisq(x, df = df, ncp = ncp_est) - (1-power)
return(uniroot(error, c(0, 1000), x = xc, df = df, power=power)$root) }
getncp=function(family,df1,df2,power) {
if (family=="f") ncp=getncp.f(df1=df1,df2=df2,power=power)
if (family=="c") ncp=getncp.c(df=df1,power=power)
return(ncp) }
percent <- function(x, digits = 0, format = "f", ...) {
paste(formatC(100 * x, format = format, digits = digits, ...), "%", sep = "")
}
pbound=function(p) pmin(pmax(p,2.2e-16),1-2.2e-16)
prop33=function(pc)
{
prop=ifelse(family=="f" & p<.05,1-pf(qf(1-pc,df1=df1, df2=df2),df1=df1, df2=df2, ncp=ncp33),NA)
prop=ifelse(family=="c" & p<.05,1-pchisq(qchisq(1-pc,df=df1), df=df1, ncp=ncp33),prop)
prop
}
stouffer=function(pp) sum(qnorm(pp),na.rm=TRUE)/sqrt(sum(!is.na(pp)))
###############################################################################
# 2. Process data ############################################################
###############################################################################
# Note: due to reliance on the pcurve-app function, z-scores are pasted into characters first
# and then screened to generate variables necessary for further computation
zvalues.input = paste("z=", zvalues.input, sep="")
filek = "input"
raw = zvalues.input
raw=tolower(raw)
ktot=length(raw)
k=seq(from=1,to=length(raw))
stat=substring(raw,1,1)
test=ifelse(stat=="r","t",stat)
# Create family
family=test
family=ifelse(test=="t","f",family)
family=ifelse(test=="z","c",family)
#family: f,c converting t-->f and z-->c
# Find comma,parentheses,equal sign
par1 =str_locate(raw,"\\(")[,1]
par2 =str_locate(raw,"\\)")[,1]
comma=str_locate(raw,",")[,1]
eq =str_locate(raw,"=")[,1]
# DF for t-tests
df=as.numeric(ifelse(test=="t",substring(raw,par1+1,par2 -1),NA))
# DF1
df1=as.numeric(ifelse(test=="f",substring(raw,par1+1,comma-1),NA))
df1=as.numeric(ifelse(test=="z",1,df1))
df1=as.numeric(ifelse(test=="t",1,df1))
df1=as.numeric(ifelse(test=="c",substring(raw,par1+1,par2 -1),df1))
# DF2
df2=as.numeric(ifelse(test=="f",substring(raw,comma+1,par2-1),NA))
df2=as.numeric(ifelse(test=="t",df,df2))
equal=abs(as.numeric(substring(raw,eq+1)))
value=ifelse((stat=="f" | stat=="c"),equal,NA)
value=ifelse(stat=="r", (equal/(sqrt((1-equal**2)/df2)))**2,value)
value=ifelse(stat=="t", equal**2 ,value)
value=ifelse(stat=="z", equal**2 ,value)
p=ifelse(family=="f",1-pf(value,df1=df1,df2=df2),NA)
p=ifelse(family=="c",1-pchisq(value,df=df1),p)
p=pbound(p) #Bound it to level of precision, see function 3 above
ksig= sum(p<.05,na.rm=TRUE) #significant studies
khalf=sum(p<.025,na.rm=TRUE) #half p-curve studies
if (ksig <= 2){
stop("Two or less significant (p<0.05) effect sizes were detected, so p-curve analysis cannot be conducted.")
}
##############################################################################
# 3. PP-values ###############################################################
##############################################################################
# Right Skew, Full p-curve
ppr=as.numeric(ifelse(p<.05,20*p,NA))
ppr=pbound(ppr)
# Right Skew, half p-curve
ppr.half=as.numeric(ifelse(p<.025,40*p,NA))
ppr.half=pbound(ppr.half)
# Power of 33%
ncp33=mapply(getncp,df1=df1,df2=df2,power=1/3,family=family)
# Full-p-curve
pp33=ifelse(family=="f" & p<.05,3*(pf(value, df1=df1, df2=df2, ncp=ncp33)-2/3),NA)
pp33=ifelse(family=="c" & p<.05,3*(pchisq(value, df=df1, ncp=ncp33)-2/3),pp33)
pp33=pbound(pp33)
# half p-curve
prop25=3*prop33(.025)
prop25.sig=prop25[p<.05]
#Compute pp-values for the half
pp33.half=ifelse(family=="f" & p<.025, (1/prop25)*(pf(value,df1=df1,df2=df2,ncp=ncp33)-(1-prop25)),NA)
pp33.half=ifelse(family=="c" & p<.025, (1/prop25)*(pchisq(value,df=df1, ncp=ncp33)-(1-prop25)),pp33.half)
pp33.half=pbound(pp33.half)
##############################################################################
# 4. Stouffer & Binomial test ################################################
##############################################################################
# Convert pp-values to Z scores, using Stouffer function above
Zppr = stouffer(ppr)
Zpp33 = stouffer(pp33)
Zppr.half = stouffer(ppr.half)
Zpp33.half = stouffer(pp33.half)
# Overall p-values from Stouffer test
p.Zppr = pnorm(Zppr)
p.Zpp33 = pnorm(Zpp33)
p.Zppr.half = pnorm(Zppr.half)
p.Zpp33.half = pnorm(Zpp33.half)
# Save results to file
main.results=as.numeric(c(ktot, ksig, khalf, Zppr,
p.Zppr, Zpp33, p.Zpp33, Zppr.half,
p.Zppr.half, Zpp33.half, p.Zpp33.half))
# BINOMIAL
# Observed share of p<.025
prop25.obs=sum(p<.025)/sum(p<.05)
# Flat null
binom.r=1-pbinom(q=prop25.obs*ksig- 1, prob=.5, size=ksig)
# Power of 33% null
binom.33=ppoibin(kk=prop25.obs*ksig,pp=prop25[p<.05])
# Save binomial results
binomial=c(mean(prop25.sig), prop25.obs, binom.r, binom.33)
# Beautifyier Function
cleanp=function(p)
{
p.clean=round(p,4) #Round it
p.clean=substr(p.clean,2,6) #Drop the 0
p.clean=paste0("= ",p.clean)
if (p < .0001) p.clean= " < .0001"
if (p > .9999) p.clean= " > .9999"
return(p.clean)
}
#If there are zero p<.025, change Stouffer values for half-p-curve tests for "N/A" messages
if (khalf==0) {
Zppr.half ="N/A"
p.Zppr.half ="=N/A"
Zpp33.half ="N/A"
p.Zpp33.half ="=N/A"
}
#If there are more than 1 p<.025, round the Z and beutify the p-values
if (khalf>0) {
Zppr.half =round(Zppr.half,2)
Zpp33.half =round(Zpp33.half,2)
p.Zppr.half=cleanp(p.Zppr.half)
p.Zpp33.half=cleanp(p.Zpp33.half)
}
#Clean results for full test
Zppr=round(Zppr,2)
Zpp33=round(Zpp33,2)
p.Zppr=cleanp(p.Zppr)
p.Zpp33=cleanp(p.Zpp33)
binom.r=cleanp(binom.r)
binom.33=cleanp(binom.33)
################################################
# 5. Power ####################################
################################################
powerfit=function(power_est)
{
ncp_est=mapply(getncp,df1=df1,df2=df2,power=power_est,family=family)
pp_est=ifelse(family=="f" & p<.05,(pf(value,df1=df1,df2=df2,ncp=ncp_est)-(1-power_est))/power_est,NA)
pp_est=ifelse(family=="c" & p<.05,(pchisq(value,df=df1,ncp=ncp_est)-(1-power_est))/power_est,pp_est)
pp_est=pbound(pp_est)
return(stouffer(pp_est))
}
fit=c()
fit=abs(powerfit(.051))
for (i in 6:99) fit=c(fit,abs(powerfit(i/100)))
mini=match(min(fit,na.rm=TRUE),fit)
hat=(mini+4)/100
x.power=seq(from=5,to=99)/100
get.power_pct =function(pct) {
#Function that finds power that gives p-value=pct for the Stouffer test
#for example, get.power_pct(.5) returns the level of power that leads to p=.5 for the stouffer test.
#half the time we would see p-curves more right skewed than the one we see, and half the time
#less right-skewed, if the true power were that get.power_pct(.5). So it is the median estimate of power
#similarliy, get.power_pct(.1) gives the 10th percentile estimate of power...
#Obtain the normalized equivalent of pct, e.g., for 5% it is -1.64, for 95% it is 1.64
z=qnorm(pct) #convert to z because powerfit() outputs a z-score.
#Quantify gap between computed p-value and desired pct
error = function(power_est, z) powerfit(power_est) - z
#Find the value of power that makes that gap zero, (root)
return(uniroot(error, c(.0501, .99),z)$root) }
# Boundary conditions
p.power.05=pnorm(powerfit(.051)) #Proability p-curve would be at least at right-skewed if power=.051
p.power.99=pnorm(powerfit(.99)) #Proability p-curve would be at least at right-skewed if power=.99
# Lower end of ci
if (p.power.05<=.95) power.ci.lb=.05
if (p.power.99>=.95) power.ci.lb=.99
if (p.power.05>.95 && p.power.99<.95) power.ci.lb=get.power_pct(.95)
# Higher end of CI
if (p.power.05<=.05) power.ci.ub=.05
if (p.power.99>=.05) power.ci.ub=.99
if (p.power.05>.05 && p.power.99<.05) power.ci.ub=get.power_pct(.05)
# Save power fit
power_results=c(power.ci.lb,hat,power.ci.ub)
##############################################################################
# 6. Plot ###################################################################
##############################################################################
# Green line (Expected p-curve for 33% power)
gcdf1=prop33(.01)
gcdf2=prop33(.02)
gcdf3=prop33(.03)
gcdf4=prop33(.04)
green1=mean(gcdf1,na.rm=TRUE)*3
green2=mean(gcdf2-gcdf1,na.rm=TRUE)*3
green3=mean(gcdf3-gcdf2,na.rm=TRUE)*3
green4=mean(gcdf4-gcdf3,na.rm=TRUE)*3
green5=mean(1/3-gcdf4,na.rm=TRUE)*3
green=100*c(green1,green2,green3,green4,green5)
# Blue line (observed p-curve)
ps=ceiling(p[p<.05]*100)/100
blue=c()
for (i in c(.01,.02,.03,.04,.05)) blue=c(blue,sum(ps==i,na.rm=TRUE)/ksig*100)
# Red line
red=c(20,20,20,20,20)
# Make the graph
x = c(.01,.02,.03,.04,.05)
par(mar=c(6,5.5,1.5,3))
moveup=max(max(blue[2:5])-66,0)
ylim=c(0,105+moveup)
legend.top=100+moveup
plot(x,blue, type='l', col='dodgerblue2', main="",
lwd=2, xlab="", ylab="", xaxt="n",yaxt="n", xlim=c(0.01,0.051),
ylim=ylim, bty='L', las=1,axes=F)
x_=c(".01",".02",".03",".04",".05")
axis(1,at=x,labels=x_)
y_=c("0%","25%","50%","75%","100%")
y=c(0,25,50,75,100)
axis(2,at=y,labels=y_,las=1,cex.axis=1.2)
mtext("Percentage of test results",font=2,side=2,line=3.85,cex=1.25)
mtext("p ",font=4,side=1,line=2.3,cex=1.25)
mtext(" -value", font=2,side=1,line=2.3,cex=1.25)
points(x,blue,type="p",pch=20,bg="dodgerblue2",col="dodgerblue2")
text(x+.00075,blue+3.5,percent(round(blue)/100),col='black', cex=.75)
lines(x,red, type='l', col='firebrick2', lwd=1.5, lty=3)
lines(x,green, type='l', col='springgreen4', lwd=1.5, lty=5)
tab1=.017 #Labels for line at p=.023 in x-axis
tab2=tab1+.0015 #Test results and power esimates at tab1+.0015
gap1=9 #between labels
gap2=4 #between lable and respective test (e.g., "OBserved p-curve" and "power estimate")
font.col='gray44'
text.blue=paste0("Power estimate: ",percent(hat),", CI(",
percent(power.ci.lb),",",
percent(power.ci.ub),")")
text(tab1,legend.top, adj=0,cex=.85,bquote("Observed "*italic(p)*"-curve"))
text(tab2,legend.top-gap2,adj=0,cex=.68,text.blue,col=font.col)
text.red=bquote("Tests for right-skewness: "*italic(p)*""[Full]~.(p.Zppr)*", "*italic(p)*""[Half]~.(p.Zppr.half))
#note: .() within bquote prints the value rather than the variable name
text(tab1,legend.top-gap1, adj=0,cex=.85, "Null of no effect" )
text(tab2,legend.top-gap1-gap2, adj=0,cex=.68, text.red, col=font.col )
text.green=bquote("Tests for flatness: "*italic(p)*""[Full]~.(p.Zpp33)*", "*italic(p)*""[half]~.(p.Zpp33.half)*", "*italic(p)*""[Binomial]~.(binom.33))
text(tab1,legend.top-2*gap1, adj=0,cex=.85,"Null of 33% power")
text(tab2,legend.top-2*gap1-gap2, adj=0,cex=.68,text.green,col=font.col)
segments(x0=tab1-.005,x1=tab1-.001,y0=legend.top,y1=legend.top, col='dodgerblue2',lty=1,lwd=1.5)
segments(x0=tab1-.005,x1=tab1-.001,y0=legend.top-gap1, y1=legend.top-gap1,col='firebrick2',lty=3,lwd=1.5)
segments(x0=tab1-.005,x1=tab1-.001,y0=legend.top-2*gap1,y1=legend.top-2*gap1,col='springgreen4',lty=2,lwd=1.5)
rect(tab1-.0065,legend.top-2*gap1-gap2-3,tab1+.032,legend.top+3,border='gray85')
msgx=bquote("Note: The observed "*italic(p)*"-curve includes "*.(ksig)*
" statistically significant ("*italic(p)*" < .05) results, of which "*.(khalf)*
" are "*italic(p)*" < .025.")
mtext(msgx,side=1,line=4,cex=.65,adj=0)
kns=ktot-ksig
if (kns==0) ns_msg="There were no non-significant results entered."
if (kns==1) ns_msg=bquote("There was one additional result entered but excluded from "*italic(p)*"-curve because it was "*italic(p)*" > .05.")
if (kns>1) ns_msg=bquote("There were "*.(kns)*" additional results entered but excluded from "*italic(p)*"-curve because they were "*italic(p)*" > .05.")
mtext(ns_msg,side=1,line=4.75,cex=.65,adj=0)
##############################################################################
# 7 Save Calculations #######################################################
##############################################################################
# table_calc
table_calc=data.frame(raw, p, ppr, ppr.half, pp33, pp33.half,
qnorm(ppr), qnorm(ppr.half), qnorm(pp33), qnorm(pp33.half))
headers1=c("Entered statistic","p-value", "ppr", "ppr half", "pp33%","pp33 half",
"Z-R","Z-R half","Z-33","z-33 half")
table_calc=setNames(table_calc,headers1)
# table_figure
headers2=c("p-value","Observed (blue)","Power 33% (Green)", "Flat (Red)")
table_figure=setNames(data.frame(x,blue,green,red),headers2)
################################################
# 8. Cumulative p-curves (Deprecated) ##########
################################################
#7.1 FUNCTION THAT RECOMPUTES OVERALL STOUFFER TEST WITHOUT (K) MOST EXTREME VALUES, ADJUSTING THE UNIFORM TO ONLY INCLUDE RANGE THAT REMAINS
dropk=function(pp,k,droplow)
{
#Syntax:
#pp: set of pp-values to analyze sensitivity to most extremes
#k: # of most extreme values to exclude
#dropsmall: 1 to drop smallest, 0 to drop largest
pp=pp[!is.na(pp)] #Drop missing values
n=length(pp) #See how many studies are left
pp=sort(pp) #Sort the pp-value from small to large
if (k==0) ppk=pp #If k=0 do nothing for nothing is being dropped
#If we are dropping low values
if (droplow==1 & k>0)
{
#Eliminate lowest k from the vector of pp-values
ppk=(pp[(1+k):n])
ppmin=min(pp[k],k/(n+1)) #Boundary used to define possible range of values after exclusion
ppk=(ppk-ppmin)/(1-ppmin) #Take the k+1 smallest pp-value up to the highest, subtract from each the boundary value, divide by the range, ~U(0,1) under the null
#This is explained in Supplement 1 of Simonsohn, Simmons Nelson, JEPG 2016 "Better p-curves" paper. See https://osf.io/mbw5g/
}
#If we are dropping high values
if (droplow==0 & k>0)
{
#Eliminate lowest k from the vector of pp-values
ppk=pp[1:(n-k)]
ppmax=max(pp[n-k+1],(n-k)/(n+1)) #Find new boundary of range
ppk=ppk/ppmax #Redefine range to make U(0,1)
}
#In case of a tie with two identical values we would have the ppk be 0 or 1, let's replace that with almost 0 and almost 1
ppk=pmax(ppk,.00001) #Adds small constant to the smallest redefined p-value, avoids problem if dropped p-value is "equal" to next highest, then that pp-value becomes 0
ppk=pmin(ppk,.99999) #Subtract small constant to the largest redefined pp-value, same reason
Z=sum(qnorm(ppk))/sqrt(n-k)
return(pnorm(Z))
} #End function dropk
#7.2 Apply function, in loop with increasing number of exclusions, to full p-curve
#Empty vectors for results
droplow.r=droplow.33=drophigh.r=drophigh.33=c()
#Loop over full p-curves
for (i in 0:(round(ksig/2)-1))
{
#Drop the lowest k studies in terms of respective overall test
#Right skew
droplow.r= c(droplow.r, dropk(pp=ppr,k=i,droplow=1))
drophigh.r=c(drophigh.r, dropk(pp=ppr,k=i,droplow=0))
#Power of 33%
droplow.33=c(droplow.33, dropk(pp=pp33,k=i,droplow=1))
drophigh.33=c(drophigh.33, dropk(pp=pp33,k=i,droplow=0))
}
#Half p-curves
if (khalf>0)
{
droplow.halfr=drophigh.halfr=c()
for (i in 0:(round(khalf/2)-1))
{
#Drop the lowest k studies in terms of respective overall test
droplow.halfr= c(droplow.halfr, dropk(pp=ppr.half,k=i,droplow=1))
drophigh.halfr=c(drophigh.halfr, dropk(pp=ppr.half,k=i,droplow=0))
} #End loop
}#End if that runs calculations only if khalf>0
#7.3 FUNCTION THAT DOES THE PLOT OF RESULTS
plotdrop=function(var,col)
{
k=length(var)
#Plot the dots
plot(0:(k-1),var,xlab="",ylab="",type="b",yaxt="n",xaxt="n",main="",
cex.main=1.15,ylim=c(0,1),col=col)
#Add marker in results with 0 drops
points(0,var[1],pch=19,cex=1.6)
#Red line at p=.05
abline(h=.05,col="red")
#Y-axis value labels
axis(2,c(.05,2:9/10),labels=c('.05','.2','.3','.4','.5','6','7','.8','.9'),las=1,cex.axis=1.5)
axis(1,c(0:(k-1)),las=1,cex.axis=1.4)
}
######################################################################################
# 9. Effect Estimation ###############################################################
######################################################################################
if (effect.estimation == TRUE){
# Define ci.to.t function
ci.to.t = function(TE, lower, upper, n){
z.to.d = function(z, n){
d = (2*z)/sqrt(n)
return(abs(d))
}
ci.to.p = function(est, lower, upper){
SE = (upper-lower)/(2*1.96)
z = abs(est/SE)
p = exp(-0.717*z - 0.416*z^2)
return(p)
}
d.to.t = function(d, n){
df = n-2
t = (d*sqrt(df))/2
return(t)
}
p = ci.to.p(TE, lower, upper)
z = abs(qnorm(p/2))
d = z.to.d(z, n)
t = d.to.t(d, n)
return(t)
}
#Function 13 - loss function
loss=function(t_obs,df_obs,d_est) {
#1.Convert all ts to the same sign (for justification see Supplement 5)
t_obs=abs(t_obs)
#2 Compute p-values
p_obs=2*(1-pt(t_obs,df=df_obs))
#3 Keep significant t-values and corresponding df.
t.sig=subset(t_obs,p_obs<.05)
df.sig=subset(df_obs,p_obs<.05)
#4.Compute non-centrality parameter implied by d_est and df_obs
#df+2 is total N.
#Becuase the noncentrality parameter for the student distribution is ncp=sqrt(n/2)*d,
#we add 2 to d.f. to get N, divide by 2 to get n, and by 2 again for ncp, so -->df+2/4
ncp_est=sqrt((df.sig+2)/4)*d_est
#5.Find critical t-value for p=.05 (two-sided)
#this is used below to compute power, it is a vector as different tests have different dfs
#and hence different critical values
tc=qt(.975,df.sig)
#4.Find power for ncp given tc, again, this is a vector of implied power, for ncp_est, for each test
power_est=1-pt(tc,df.sig,ncp_est)
#5.Compute pp-values
#5.1 First get the overall probability of a t>tobs, given ncp
p_larger=pt(t.sig,df=df.sig,ncp=ncp_est)
#5.2 Now, condition on p<.05
ppr=(p_larger-(1-power_est))/power_est #this is the pp-value for right-skew
#6. Compute the gap between the distribution of observed pp-values and a uniform distribution 0,1
KSD=ks.test(ppr,punif)$statistic #this is the D statistic outputted by the KS test against uniform
return(KSD)
}
if(missing(N)){
stop("If 'effect.estimation=TRUE', argument 'N' must be provided.")
}
if (length(N) != length(metaobject$TE)){
stop("N must be of same length as the number of studies contained in x.")
}
lower = metaobject$TE - (metaobject$seTE*1.96)
upper = metaobject$TE + (metaobject$seTE*1.96)
t_obs = ci.to.t(metaobject$TE, lower, upper, N)
df_obs = N-2
#Results will be stored in these vectors, create them first
loss.all=c()
di=c()
#Compute loss for effect sizes between d=c(dmin,dmax) in steps of .01
for (i in 0:((dmax-dmin)*100))
{
d=dmin+i/100 #effect size being considered
di=c(di,d) #add it to the vector (kind of silly, but kept for symmetry)
options(warn=-1) #turn off warning becuase R does not like its own pt() function!
loss.all=c(loss.all,loss(df_obs=df_obs,t_obs=t_obs,d_est=d))
#apply loss function so that effect size, store result
options(warn=0) #turn warnings back on
}
#find the effect leading to smallest loss in that set, that becomes the starting point in the optimize command
imin=match(min(loss.all),loss.all) #which i tested effect size lead to the overall minimum?
dstart=dmin+imin/100 #convert that i into a d.
#optimize around the global minimum
dhat=optimize(loss,c(dstart-.1,dstart+.1), df_obs=df_obs,t_obs=t_obs)
options(warn=-0)
#Plot results
plot(di,loss.all,xlab="Effect size\nCohen-d", ylab="Loss (D stat in KS test)",ylim=c(0,1), main="How well does each effect size fit? (lower is better)")
points(dhat$minimum,dhat$objective,pch=19,col="red",cex=2)
text(dhat$minimum,dhat$objective-.08,paste0("p-curve's estimate of effect size:\nd=",round(dhat$minimum,3)),col="red")
}
######################################################################################
# 10. Prepare Results for Return #####################################################
######################################################################################
# Get results
main.results = round(main.results, 3)
ktotal = round(main.results[1]) # Get the total number of inserted TEs
k.sign = round(main.results[2]) # Get the total number of significant TEs
k.025 = round(main.results[3]) # Get the number of p<0.25 TEs
skew.full.z = main.results[4] # Get the Z-score for the full curve skewness test
skew.full.p = main.results[5] # Get the p-value for the full curve skewness test
flat.full.z = main.results[6] # Get the Z-score for the full curve flatness test
flat.full.p = main.results[7] # Get the p-value for the full curve flatness test
skew.half.z = main.results[8] # Get the Z-score for the half curve skewness test
skew.half.p = main.results[9] # Get the p-value for the half curve skewness test
flat.half.z = main.results[10] # Get the Z-score for the half curve flatness test
flat.half.p = main.results[11] # Get the p-value for the half curve flatness test
skew.binomial.p = round(binomial[3], 3) # Get the skewness binomial p-value
flat.binomial.p = round(binomial[4], 3) # Get the flatness binomial p-value
# Make data.frame
skewness = c(skew.binomial.p, skew.full.z, skew.full.p, skew.half.z, skew.half.p)
flatness = c(flat.binomial.p, flat.full.z, flat.full.p, flat.half.z, flat.half.p)
colnames.df = c("pBinomial", "zFull", "pFull", "zHalf", "pHalf")
rownames.df = c("Right-skewness test", "Flatness test")
pcurveResults = rbind(skewness, flatness)
colnames(pcurveResults) = colnames.df
rownames(pcurveResults) = rownames.df
# Power results
power_results = round(power_results, 3)
powerEstimate = power_results[2]
powerLower = power_results[1]
powerUpper = power_results[3]
Power = as.data.frame(cbind(powerEstimate, powerLower, powerUpper))
rownames(Power) = ""
# Presence and absence of evidential value
# - If the half p-curve test is right-skewed with p<.05 or both the half and full test
# are right-skewed with p<.1, then p-curve analysis indicates the presence of evidential value
# - Evidential value is inadequate or absent if the 33% power test is p<.05 for the full p-curve
# or both the half p-curve and binomial 33% power test are p<.1
if (skew.half.p < 0.05 | (skew.half.p < 0.1 & skew.full.p < 0.1)){
presence.ev = "yes"
} else {
presence.ev = "no"
}
if (flat.full.p < 0.05 | (flat.half.p < 0.1 & flat.binomial.p < 0.1)){
absence.ev = "yes"
} else {
absence.ev = "no"
}
# Plot Data
PlotData = round(table_figure, 3)
# Input Data
table_calc[,1] = NULL
colnames(table_calc) = c("p", "ppSkewFull", "ppSkewHalf", "ppFlatFull", "ppFlatHalf", "zSkewFull", "zSkewHalf",
"zFlatFull", "zFlatHalf")
Input = cbind(metaobject$TE, round(table_calc,3))
rownames(Input) = paste(1:length(metaobject$TE), metaobject$studlab)
colnames(Input)[1] = "TE"
if (effect.estimation==TRUE){
dEstimate = round(dhat$minimum, 3)
return.list = list("pcurveResults" = pcurveResults,
"Power" = Power,
"PlotData" = PlotData,
"Input" = Input,
"EvidencePresent" = presence.ev,
"EvidenceAbsent" = absence.ev,
"kInput" = ktot,
"kAnalyzed" = k.sign,
"kp0.25" = k.025,
"dEstimate" = dEstimate,
"I2" = metaobject$I2,
"class.meta.object" = class(metaobject)[1])
class(return.list) = c("pcurve", "effect.estimation")
} else {
return.list = list("pcurveResults" = pcurveResults,
"Power" = Power,
"PlotData" = PlotData,
"Input" = Input,
"EvidencePresent" = presence.ev,
"EvidenceAbsent" = absence.ev,
"kInput" = ktot,
"kAnalyzed" = k.sign,
"kp0.25" = k.025,
"I2" = metaobject$I2,
"class.meta.object" = class(metaobject)[1])
class(return.list) = c("pcurve", "no.effect.estimation")
}
cat(" ", "\n")
invisible(return.list)
return.list
}
MathiasHarrer/dmetar documentation built on April 4, 2024, 6:57 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions pcurve
Documented in pcurve
#' Perform a \emph{p}-curve analysis
#'
#' This function performs a \eqn{p}-curve analysis using a \code{meta} object or calculated effect size data.
#'
#' @usage pcurve(x, effect.estimation = FALSE, N, dmin = 0, dmax = 1)
#'
#' @param x Either an object of class \code{meta}, generated by the \code{metagen}, \code{metacont},
#' \code{metacor}, \code{metainc}, or \code{metabin} function, or a dataframe containing the calculated effect size
#' (named \code{TE}, log-transformed if based on a ratio), standard error (named \code{seTE}) and study label (named \code{studlab})
#' for each study.
#' @param effect.estimation Logical. Should the true effect size underlying the \emph{p}-curve be estimated?
#' If set to \code{TRUE}, a vector containing the total sample size for each study must be provided for
#' \code{N}. \code{FALSE} by default.
#' @param N A numeric vector of same length as the number of effect sizes included in \code{x} specifying the
#' total sample size \eqn{N} corresponding to each effect. Only needed if \code{effect.estimation = TRUE}.
#' @param dmin If \code{effect.estimation = TRUE}: lower limit for the effect size (\eqn{d}) space in which
#' the true effect size should be searched. Must be greater or equal to 0. Default is 0.
#' @param dmax If \code{effect.estimation = TRUE}: upper limit for the effect size (\eqn{d}) space in which
#' the true effect size should be searched. Must be greater than 0. Default is 1.
#'
#' @details
#' \strong{P-curve Analysis}
#'
#' \eqn{P}-curve analysis (Simonsohn, Simmons & Nelson, 2014, 2015) has been proposed as a method
#' to detect \eqn{p}-hacking and publication bias in meta-analyses.
#'
#' \eqn{P}-Curve assumes that publication bias
#' is not only generated because researchers do not publish non-significant results,
#' but also because analysts “play” around with their data ("\eqn{p}-hacking"; e.g., selectively removing outliers,
#' choosing different outcomes, controlling for different variables) until a non-significant
#' finding becomes significant (i.e., \eqn{p<0.05}).
#'
#' The method assumes that for a specific research
#' question, \eqn{p}-values smaller 0.05 of included studies should follow a right-skewed distribution
#' if a true effect exists, even when the power in single studies was (relatively) low. Conversely,
#' a left-skewed \eqn{p}-value distribution indicates the presence of \eqn{p}-hacking and absence of
#' a true underlying effect. To control for "ambitious" \eqn{p}-hacking, \eqn{P}-curve also incorporates a
#' "half-curve" test (Simonsohn, Simmons & Nelson, 2014, 2015).
#'
#' Simonsohn et al. (2014)
#' stress that \eqn{p}-curve analysis should only be used for test statistics which were actually of interest
#' in the context of the included study, and that a detailed table documenting the reported results
#' used in for the \eqn{p}-curve analysis should be created before communicating
#' results (\href{http://www.p-curve.com/Supplement/}{link}).
#'
#' \strong{Implementation in the function}
#'
#' To generate the \eqn{p}-curve and conduct the analysis, this function reuses parts of the \emph{R} code underlying
#' the \href{http://p-curve.com/app4/pcurve_app4.052.r}{P-curve App 4.052} (Simonsohn, 2017). The effect sizes
#' included in the \code{meta} object or \code{data.frame} provided for \code{x} are transformed
#' into \eqn{z}-values internally, which are then used to calculate {p}-values and conduct the
#' Stouffer and Binomial test used for the \eqn{p}-curve analysis. Interpretations of the function
#' concerning the presence or absence/inadequateness of evidential value are made according to the
#' guidelines described by Simonsohn, Simmons and Nelson (2015):
#'
#' \itemize{
#' \item \strong{Evidential value present}: The right-skewness test is significant for the half curve with
#' \eqn{p<0.05} \strong{or} the \eqn{p}-value of the right-skewness test is \eqn{<0.1} for both the half and full curve.
#' \item \strong{Evidential value absent or inadequate}: The flatness test is \eqn{p<0.05} for the full curve
#' \strong{or} the flatness test for the half curve and the binomial test are \eqn{p<0.1}.
#'}
#'
#' For effect size estimation, the \code{pcurve} function implements parts of the loss function
#' presented in Simonsohn, Simmons and Nelson (2014b).
#' The function generates a loss function for candidate effect sizes \eqn{\hat{d}}, using \eqn{D}-values in
#' a Kolmogorov-Smirnov test as the metric of fit, and the value of \eqn{\hat{d}} which minimizes \eqn{D}
#' as the estimated true effect.
#'
#' It is of note that a lack of robustness of \eqn{p}-curve analysis results
#' has been noted for meta-analyses with substantial heterogeneity (van Aert, Wicherts, & van Assen, 2016).
#' Following van Aert et al., adjusted effect size estimates should only be
#' reported and interpreted for analyses with \eqn{I^2} values below 50 percent.
#' A warning message is therefore printed by
#' the \code{pcurve} function when \code{x} is of class \code{meta} and the between-study heterogeneity
#' of the meta-analysis is substantial (i.e., \eqn{I^2} greater than 50 percent).
#'
#'
#' @references Harrer, M., Cuijpers, P., Furukawa, T.A, & Ebert, D. D. (2019).
#' \emph{Doing Meta-Analysis in R: A Hands-on Guide}. DOI: 10.5281/zenodo.2551803.
#' \href{https://bookdown.org/MathiasHarrer/Doing_Meta_Analysis_in_R/pcurve.html}{Chapter 9.2}.
#'
#' Simonsohn, U., Nelson, L. D., & Simmons, J. P. (2014a). P-curve: a Key to the File-drawer.
#' \emph{Journal of Experimental Psychology, 143}(2), 534.
#'
#' Simonsohn, U., Nelson, L. D. & Simmons, J. P. (2014b). P-Curve and Effect Size:
#' Correcting for Publication Bias Using Only Significant Results.
#' \emph{Perspectives on Psychological Science 9}(6), 666–81.
#'
#' Simonsohn, U., Nelson, L. D. & Simmons, J. P. (2015). Better P-Curves: Making P-Curve
#' Analysis More Robust to Errors, Fraud, and Ambitious P-Hacking, a Reply to Ulrich and Miller (2015).
#' \emph{Journal of Experimental Psychology, 144}(6), 1146-1152.
#'
#' Simonsohn, U. (2017). R code for the P-Curve App 4.052. http://p-curve.com/app4/pcurve_app4.052.r (Accessed 2019-08-16).
#'
#' Van Aert, R. C., Wicherts, J. M., & van Assen, M. A. (2016).
#' Conducting meta-analyses based on p values: Reservations and recommendations for applying
#' \emph{p}-uniform and \emph{p}-curve. \emph{Perspectives on Psychological Science, 11}(5), 713-729.
#'
#' @author Mathias Harrer & David Daniel Ebert
#'
#' @return Returns a plot and main results of the pcurve analysis:
#' \itemize{
#' \item \strong{P-curve plot}: A plot displaying the observed \eqn{p}-curve and significance results
#' for the right-skewness and flatness test.
#' \item \strong{Number of studies}: The number of studies provided for the analysis, the number
#' of significant \eqn{p}-values included in the analysis, and the number of studies with \eqn{p<0.025}
#' used for the half-curve tests.
#' \item \strong{Test results}: The results for the right-skewness and flatness test, including the
#' \eqn{p_{binomial}} value, as well as the \eqn{z} and \eqn{p} value for the full and half-curve test.
#' \item \strong{Power Estimate}: The power estimate and 95\% confidence interval.
#' \item \strong{Evidential value}: Two lines displaying if evidential value is present and/or absent/inadequate based
#' on the results (using the guidelines by Simonsohn et al., 2015, see Details).
#' \item \strong{True effect estimate}: If \code{effect.estimation} is set to \code{TRUE}, the estimated true effect
#' \eqn{\hat{d}} is returned additionally.
#'}
#'
#'
#' If results are saved to a variable, a list of class \code{pcurve} containing the following objects is returned:
#' \itemize{
#' \item \code{pcurveResults}: A data frame containing the results for the right-skewness and flatness test, including the
#' \eqn{p_{binomial}} value, as well as the \eqn{z} and \eqn{p} value for the full and half-curve test.
#' \item \code{Power}: The power estimate and 95\% confidence interval.
#' \item \code{PlotData}: A data frame with the data used in the \eqn{p}-curve plot.
#' \item \code{Input}: A data frame containing the provided effect sizes, calculated \eqn{p}-values and individual results for each included (significant) effect.
#' \item \code{EvidencePresent}, \code{EvidenceAbsent}, \code{kInput}, \code{kAnalyzed}, \code{kp0.25}: Further results of the \eqn{p}-curve analysis, including the presence/absence of evidence interpretation,
#' and number of provided/significant/\eqn{p<0.025} studies.
#' \item \code{I2}: \eqn{I^2}-Heterogeneity of the studies provided as input (only when \code{x} is of class \code{meta}).
#' \item \code{class.meta.object}: \code{class} of the original object provided in \code{x}.
#'}
#'
#'
#' @import stringr poibin
#' @importFrom graphics abline axis lines mtext par plot points rect segments text
#' @importFrom stats as.formula hat influence ks.test optimize pbinom pchisq pf pnorm pt punif qchisq qf qnorm qt reformulate reorder setNames uniroot
#'
#' @export pcurve
#'
#' @seealso
#' \code{\link{eggers.test}}
#'
#' @examples
#' # Example 1: Use metagen object, do not estimate d
#' suppressPackageStartupMessages(library(meta))
#' data("ThirdWave")
#' meta1 = metagen(TE,seTE, studlab=ThirdWave$Author, data=ThirdWave)
#' pcurve(meta1)
#'
#' # Example 2: Provide Ns, calculate d estimate
#' N = c(105, 161, 60, 37, 141, 82, 97, 61, 200, 79, 124, 25, 166, 59, 201, 95, 166, 144)
#' pcurve(meta1, effect.estimation = TRUE, N = N)
#'
#' # Example 3: Use metacont object, calculate d estimate
#' data("amlodipine")
#' meta2 <- metacont(n.amlo, mean.amlo, sqrt(var.amlo),
#' n.plac, mean.plac, sqrt(var.plac),
#' data=amlodipine, studlab=study, sm="SMD")
#' N = amlodipine$n.amlo + amlodipine$n.plac
#' pcurve(meta2, effect.estimation = TRUE, N = N, dmin = 0, dmax = 1)
#'
#' # Example 4: Construct x object from scratch
#' sim = data.frame("studlab" = c(paste("Study_", 1:18, sep = "")),
#' "TE" = c(0.561, 0.296, 0.648, 0.362, 0.770, 0.214, 0.476,
#' 0.459, 0.343, 0.804, 0.357, 0.476, 0.638, 0.396, 0.497,
#' 0.384, 0.568, 0.415),
#' "seTE" = c(0.338, 0.297, 0.264, 0.258, 0.279, 0.347, 0.271, 0.319,
#' 0.232, 0.237, 0.385, 0.398, 0.342, 0.351, 0.296, 0.325,
#' 0.322, 0.225))
#' pcurve(sim)
pcurve = function(x, effect.estimation = FALSE, N, dmin = 0, dmax = 1){
# Rename x to metaobject, remove x
metaobject = x
rm(x)
# Stop if metaobject is not meta or does not contain TE or seTE column
if (!(class(metaobject)[1] %in% c("metagen", "metabin", "metacont", "metacor", "metainc", "meta", "metaprop"))){
for (i in 1:length(colnames(metaobject))){
te.exists = FALSE
if (colnames(metaobject)[i]=="TE"){
te.exists = TRUE
break
} else {
}
}
for (i in 1:length(colnames(metaobject))){
sete.exists = FALSE
if (colnames(metaobject)[i]=="seTE"){
sete.exists = TRUE
break
} else {
}
}
for (i in 1:length(colnames(metaobject))){
studlab.exists = FALSE
if (colnames(metaobject)[i]=="studlab"){
studlab.exists = TRUE
break
} else {
}
}
if(te.exists == FALSE | sete.exists ==FALSE | studlab.exists ==FALSE){
stop("x must be a meta-analysis object generated by meta functions or a data.frame with columns labeled studlab, TE, and seTE.")
}
}
#Disable scientific notation
options(scipen=999)
# Calculate Z
zvalues.input = abs(metaobject$TE/metaobject$seTE)
##############################################
# 1. Functions ###############################
##############################################
getncp.f =function(df1,df2, power) {
error = function(ncp_est, power, x, df1,df2) pf(x, df1 = df1, df2=df2, ncp = ncp_est) - (1-power)
xc=qf(p=.95, df1=df1,df2=df2)
return(uniroot(error, c(0, 1000), x = xc, df1 = df1,df2=df2, power=power)$root) }
getncp.c =function(df, power) {
xc=qchisq(p=.95, df=df)
error = function(ncp_est, power, x, df) pchisq(x, df = df, ncp = ncp_est) - (1-power)
return(uniroot(error, c(0, 1000), x = xc, df = df, power=power)$root) }
getncp=function(family,df1,df2,power) {
if (family=="f") ncp=getncp.f(df1=df1,df2=df2,power=power)
if (family=="c") ncp=getncp.c(df=df1,power=power)
return(ncp) }
percent <- function(x, digits = 0, format = "f", ...) {
paste(formatC(100 * x, format = format, digits = digits, ...), "%", sep = "")
}
pbound=function(p) pmin(pmax(p,2.2e-16),1-2.2e-16)
prop33=function(pc)
{
prop=ifelse(family=="f" & p<.05,1-pf(qf(1-pc,df1=df1, df2=df2),df1=df1, df2=df2, ncp=ncp33),NA)
prop=ifelse(family=="c" & p<.05,1-pchisq(qchisq(1-pc,df=df1), df=df1, ncp=ncp33),prop)
prop
}
stouffer=function(pp) sum(qnorm(pp),na.rm=TRUE)/sqrt(sum(!is.na(pp)))
###############################################################################
# 2. Process data ############################################################
###############################################################################
# Note: due to reliance on the pcurve-app function, z-scores are pasted into characters first
# and then screened to generate variables necessary for further computation
zvalues.input = paste("z=", zvalues.input, sep="")
filek = "input"
raw = zvalues.input
raw=tolower(raw)
ktot=length(raw)
k=seq(from=1,to=length(raw))
stat=substring(raw,1,1)
test=ifelse(stat=="r","t",stat)
# Create family
family=test
family=ifelse(test=="t","f",family)
family=ifelse(test=="z","c",family)
#family: f,c converting t-->f and z-->c
# Find comma,parentheses,equal sign
par1 =str_locate(raw,"\\(")[,1]
par2 =str_locate(raw,"\\)")[,1]
comma=str_locate(raw,",")[,1]
eq =str_locate(raw,"=")[,1]
# DF for t-tests
df=as.numeric(ifelse(test=="t",substring(raw,par1+1,par2 -1),NA))
# DF1
df1=as.numeric(ifelse(test=="f",substring(raw,par1+1,comma-1),NA))
df1=as.numeric(ifelse(test=="z",1,df1))
df1=as.numeric(ifelse(test=="t",1,df1))
df1=as.numeric(ifelse(test=="c",substring(raw,par1+1,par2 -1),df1))
# DF2
df2=as.numeric(ifelse(test=="f",substring(raw,comma+1,par2-1),NA))
df2=as.numeric(ifelse(test=="t",df,df2))
equal=abs(as.numeric(substring(raw,eq+1)))
value=ifelse((stat=="f" | stat=="c"),equal,NA)
value=ifelse(stat=="r", (equal/(sqrt((1-equal**2)/df2)))**2,value)
value=ifelse(stat=="t", equal**2 ,value)
value=ifelse(stat=="z", equal**2 ,value)
p=ifelse(family=="f",1-pf(value,df1=df1,df2=df2),NA)
p=ifelse(family=="c",1-pchisq(value,df=df1),p)
p=pbound(p) #Bound it to level of precision, see function 3 above
ksig= sum(p<.05,na.rm=TRUE) #significant studies
khalf=sum(p<.025,na.rm=TRUE) #half p-curve studies
if (ksig <= 2){
stop("Two or less significant (p<0.05) effect sizes were detected, so p-curve analysis cannot be conducted.")
}
##############################################################################
# 3. PP-values ###############################################################
##############################################################################
# Right Skew, Full p-curve
ppr=as.numeric(ifelse(p<.05,20*p,NA))
ppr=pbound(ppr)
# Right Skew, half p-curve
ppr.half=as.numeric(ifelse(p<.025,40*p,NA))
ppr.half=pbound(ppr.half)
# Power of 33%
ncp33=mapply(getncp,df1=df1,df2=df2,power=1/3,family=family)
# Full-p-curve
pp33=ifelse(family=="f" & p<.05,3*(pf(value, df1=df1, df2=df2, ncp=ncp33)-2/3),NA)
pp33=ifelse(family=="c" & p<.05,3*(pchisq(value, df=df1, ncp=ncp33)-2/3),pp33)
pp33=pbound(pp33)
# half p-curve
prop25=3*prop33(.025)
prop25.sig=prop25[p<.05]
#Compute pp-values for the half
pp33.half=ifelse(family=="f" & p<.025, (1/prop25)*(pf(value,df1=df1,df2=df2,ncp=ncp33)-(1-prop25)),NA)
pp33.half=ifelse(family=="c" & p<.025, (1/prop25)*(pchisq(value,df=df1, ncp=ncp33)-(1-prop25)),pp33.half)
pp33.half=pbound(pp33.half)
##############################################################################
# 4. Stouffer & Binomial test ################################################
##############################################################################
# Convert pp-values to Z scores, using Stouffer function above
Zppr = stouffer(ppr)
Zpp33 = stouffer(pp33)
Zppr.half = stouffer(ppr.half)
Zpp33.half = stouffer(pp33.half)
# Overall p-values from Stouffer test
p.Zppr = pnorm(Zppr)
p.Zpp33 = pnorm(Zpp33)
p.Zppr.half = pnorm(Zppr.half)
p.Zpp33.half = pnorm(Zpp33.half)
# Save results to file
main.results=as.numeric(c(ktot, ksig, khalf, Zppr,
p.Zppr, Zpp33, p.Zpp33, Zppr.half,
p.Zppr.half, Zpp33.half, p.Zpp33.half))
# BINOMIAL
# Observed share of p<.025
prop25.obs=sum(p<.025)/sum(p<.05)
# Flat null
binom.r=1-pbinom(q=prop25.obs*ksig- 1, prob=.5, size=ksig)
# Power of 33% null
binom.33=ppoibin(kk=prop25.obs*ksig,pp=prop25[p<.05])
# Save binomial results
binomial=c(mean(prop25.sig), prop25.obs, binom.r, binom.33)
# Beautifyier Function
cleanp=function(p)
{
p.clean=round(p,4) #Round it
p.clean=substr(p.clean,2,6) #Drop the 0
p.clean=paste0("= ",p.clean)
if (p < .0001) p.clean= " < .0001"
if (p > .9999) p.clean= " > .9999"
return(p.clean)
}
#If there are zero p<.025, change Stouffer values for half-p-curve tests for "N/A" messages
if (khalf==0) {
Zppr.half ="N/A"
p.Zppr.half ="=N/A"
Zpp33.half ="N/A"
p.Zpp33.half ="=N/A"
}
#If there are more than 1 p<.025, round the Z and beutify the p-values
if (khalf>0) {
Zppr.half =round(Zppr.half,2)
Zpp33.half =round(Zpp33.half,2)
p.Zppr.half=cleanp(p.Zppr.half)
p.Zpp33.half=cleanp(p.Zpp33.half)
}
#Clean results for full test
Zppr=round(Zppr,2)
Zpp33=round(Zpp33,2)
p.Zppr=cleanp(p.Zppr)
p.Zpp33=cleanp(p.Zpp33)
binom.r=cleanp(binom.r)
binom.33=cleanp(binom.33)
################################################
# 5. Power ####################################
################################################
powerfit=function(power_est)
{
ncp_est=mapply(getncp,df1=df1,df2=df2,power=power_est,family=family)
pp_est=ifelse(family=="f" & p<.05,(pf(value,df1=df1,df2=df2,ncp=ncp_est)-(1-power_est))/power_est,NA)
pp_est=ifelse(family=="c" & p<.05,(pchisq(value,df=df1,ncp=ncp_est)-(1-power_est))/power_est,pp_est)
pp_est=pbound(pp_est)
return(stouffer(pp_est))
}
fit=c()
fit=abs(powerfit(.051))
for (i in 6:99) fit=c(fit,abs(powerfit(i/100)))
mini=match(min(fit,na.rm=TRUE),fit)
hat=(mini+4)/100
x.power=seq(from=5,to=99)/100
get.power_pct =function(pct) {
#Function that finds power that gives p-value=pct for the Stouffer test
#for example, get.power_pct(.5) returns the level of power that leads to p=.5 for the stouffer test.
#half the time we would see p-curves more right skewed than the one we see, and half the time
#less right-skewed, if the true power were that get.power_pct(.5). So it is the median estimate of power
#similarliy, get.power_pct(.1) gives the 10th percentile estimate of power...
#Obtain the normalized equivalent of pct, e.g., for 5% it is -1.64, for 95% it is 1.64
z=qnorm(pct) #convert to z because powerfit() outputs a z-score.
#Quantify gap between computed p-value and desired pct
error = function(power_est, z) powerfit(power_est) - z
#Find the value of power that makes that gap zero, (root)
return(uniroot(error, c(.0501, .99),z)$root) }
# Boundary conditions
p.power.05=pnorm(powerfit(.051)) #Proability p-curve would be at least at right-skewed if power=.051
p.power.99=pnorm(powerfit(.99)) #Proability p-curve would be at least at right-skewed if power=.99
# Lower end of ci
if (p.power.05<=.95) power.ci.lb=.05
if (p.power.99>=.95) power.ci.lb=.99
if (p.power.05>.95 && p.power.99<.95) power.ci.lb=get.power_pct(.95)
# Higher end of CI
if (p.power.05<=.05) power.ci.ub=.05
if (p.power.99>=.05) power.ci.ub=.99
if (p.power.05>.05 && p.power.99<.05) power.ci.ub=get.power_pct(.05)
# Save power fit
power_results=c(power.ci.lb,hat,power.ci.ub)
##############################################################################
# 6. Plot ###################################################################
##############################################################################
# Green line (Expected p-curve for 33% power)
gcdf1=prop33(.01)
gcdf2=prop33(.02)
gcdf3=prop33(.03)
gcdf4=prop33(.04)
green1=mean(gcdf1,na.rm=TRUE)*3
green2=mean(gcdf2-gcdf1,na.rm=TRUE)*3
green3=mean(gcdf3-gcdf2,na.rm=TRUE)*3
green4=mean(gcdf4-gcdf3,na.rm=TRUE)*3
green5=mean(1/3-gcdf4,na.rm=TRUE)*3
green=100*c(green1,green2,green3,green4,green5)
# Blue line (observed p-curve)
ps=ceiling(p[p<.05]*100)/100
blue=c()
for (i in c(.01,.02,.03,.04,.05)) blue=c(blue,sum(ps==i,na.rm=TRUE)/ksig*100)
# Red line
red=c(20,20,20,20,20)
# Make the graph
x = c(.01,.02,.03,.04,.05)
par(mar=c(6,5.5,1.5,3))
moveup=max(max(blue[2:5])-66,0)
ylim=c(0,105+moveup)
legend.top=100+moveup
plot(x,blue, type='l', col='dodgerblue2', main="",
lwd=2, xlab="", ylab="", xaxt="n",yaxt="n", xlim=c(0.01,0.051),
ylim=ylim, bty='L', las=1,axes=F)
x_=c(".01",".02",".03",".04",".05")
axis(1,at=x,labels=x_)
y_=c("0%","25%","50%","75%","100%")
y=c(0,25,50,75,100)
axis(2,at=y,labels=y_,las=1,cex.axis=1.2)
mtext("Percentage of test results",font=2,side=2,line=3.85,cex=1.25)
mtext("p ",font=4,side=1,line=2.3,cex=1.25)
mtext(" -value", font=2,side=1,line=2.3,cex=1.25)
points(x,blue,type="p",pch=20,bg="dodgerblue2",col="dodgerblue2")
text(x+.00075,blue+3.5,percent(round(blue)/100),col='black', cex=.75)
lines(x,red, type='l', col='firebrick2', lwd=1.5, lty=3)
lines(x,green, type='l', col='springgreen4', lwd=1.5, lty=5)
tab1=.017 #Labels for line at p=.023 in x-axis
tab2=tab1+.0015 #Test results and power esimates at tab1+.0015
gap1=9 #between labels
gap2=4 #between lable and respective test (e.g., "OBserved p-curve" and "power estimate")
font.col='gray44'
text.blue=paste0("Power estimate: ",percent(hat),", CI(",
percent(power.ci.lb),",",
percent(power.ci.ub),")")
text(tab1,legend.top, adj=0,cex=.85,bquote("Observed "*italic(p)*"-curve"))
text(tab2,legend.top-gap2,adj=0,cex=.68,text.blue,col=font.col)
text.red=bquote("Tests for right-skewness: "*italic(p)*""[Full]~.(p.Zppr)*", "*italic(p)*""[Half]~.(p.Zppr.half))
#note: .() within bquote prints the value rather than the variable name
text(tab1,legend.top-gap1, adj=0,cex=.85, "Null of no effect" )
text(tab2,legend.top-gap1-gap2, adj=0,cex=.68, text.red, col=font.col )
text.green=bquote("Tests for flatness: "*italic(p)*""[Full]~.(p.Zpp33)*", "*italic(p)*""[half]~.(p.Zpp33.half)*", "*italic(p)*""[Binomial]~.(binom.33))
text(tab1,legend.top-2*gap1, adj=0,cex=.85,"Null of 33% power")
text(tab2,legend.top-2*gap1-gap2, adj=0,cex=.68,text.green,col=font.col)
segments(x0=tab1-.005,x1=tab1-.001,y0=legend.top,y1=legend.top, col='dodgerblue2',lty=1,lwd=1.5)
segments(x0=tab1-.005,x1=tab1-.001,y0=legend.top-gap1, y1=legend.top-gap1,col='firebrick2',lty=3,lwd=1.5)
segments(x0=tab1-.005,x1=tab1-.001,y0=legend.top-2*gap1,y1=legend.top-2*gap1,col='springgreen4',lty=2,lwd=1.5)
rect(tab1-.0065,legend.top-2*gap1-gap2-3,tab1+.032,legend.top+3,border='gray85')
msgx=bquote("Note: The observed "*italic(p)*"-curve includes "*.(ksig)*
" statistically significant ("*italic(p)*" < .05) results, of which "*.(khalf)*
" are "*italic(p)*" < .025.")
mtext(msgx,side=1,line=4,cex=.65,adj=0)
kns=ktot-ksig
if (kns==0) ns_msg="There were no non-significant results entered."
if (kns==1) ns_msg=bquote("There was one additional result entered but excluded from "*italic(p)*"-curve because it was "*italic(p)*" > .05.")
if (kns>1) ns_msg=bquote("There were "*.(kns)*" additional results entered but excluded from "*italic(p)*"-curve because they were "*italic(p)*" > .05.")
mtext(ns_msg,side=1,line=4.75,cex=.65,adj=0)
##############################################################################
# 7 Save Calculations #######################################################
##############################################################################
# table_calc
table_calc=data.frame(raw, p, ppr, ppr.half, pp33, pp33.half,
qnorm(ppr), qnorm(ppr.half), qnorm(pp33), qnorm(pp33.half))
headers1=c("Entered statistic","p-value", "ppr", "ppr half", "pp33%","pp33 half",
"Z-R","Z-R half","Z-33","z-33 half")
table_calc=setNames(table_calc,headers1)
# table_figure
headers2=c("p-value","Observed (blue)","Power 33% (Green)", "Flat (Red)")
table_figure=setNames(data.frame(x,blue,green,red),headers2)
################################################
# 8. Cumulative p-curves (Deprecated) ##########
################################################
#7.1 FUNCTION THAT RECOMPUTES OVERALL STOUFFER TEST WITHOUT (K) MOST EXTREME VALUES, ADJUSTING THE UNIFORM TO ONLY INCLUDE RANGE THAT REMAINS
dropk=function(pp,k,droplow)
{
#Syntax:
#pp: set of pp-values to analyze sensitivity to most extremes
#k: # of most extreme values to exclude
#dropsmall: 1 to drop smallest, 0 to drop largest
pp=pp[!is.na(pp)] #Drop missing values
n=length(pp) #See how many studies are left
pp=sort(pp) #Sort the pp-value from small to large
if (k==0) ppk=pp #If k=0 do nothing for nothing is being dropped
#If we are dropping low values
if (droplow==1 & k>0)
{
#Eliminate lowest k from the vector of pp-values
ppk=(pp[(1+k):n])
ppmin=min(pp[k],k/(n+1)) #Boundary used to define possible range of values after exclusion
ppk=(ppk-ppmin)/(1-ppmin) #Take the k+1 smallest pp-value up to the highest, subtract from each the boundary value, divide by the range, ~U(0,1) under the null
#This is explained in Supplement 1 of Simonsohn, Simmons Nelson, JEPG 2016 "Better p-curves" paper. See https://osf.io/mbw5g/
}
#If we are dropping high values
if (droplow==0 & k>0)
{
#Eliminate lowest k from the vector of pp-values
ppk=pp[1:(n-k)]
ppmax=max(pp[n-k+1],(n-k)/(n+1)) #Find new boundary of range
ppk=ppk/ppmax #Redefine range to make U(0,1)
}
#In case of a tie with two identical values we would have the ppk be 0 or 1, let's replace that with almost 0 and almost 1
ppk=pmax(ppk,.00001) #Adds small constant to the smallest redefined p-value, avoids problem if dropped p-value is "equal" to next highest, then that pp-value becomes 0
ppk=pmin(ppk,.99999) #Subtract small constant to the largest redefined pp-value, same reason
Z=sum(qnorm(ppk))/sqrt(n-k)
return(pnorm(Z))
} #End function dropk
#7.2 Apply function, in loop with increasing number of exclusions, to full p-curve
#Empty vectors for results
droplow.r=droplow.33=drophigh.r=drophigh.33=c()
#Loop over full p-curves
for (i in 0:(round(ksig/2)-1))
{
#Drop the lowest k studies in terms of respective overall test
#Right skew
droplow.r= c(droplow.r, dropk(pp=ppr,k=i,droplow=1))
drophigh.r=c(drophigh.r, dropk(pp=ppr,k=i,droplow=0))
#Power of 33%
droplow.33=c(droplow.33, dropk(pp=pp33,k=i,droplow=1))
drophigh.33=c(drophigh.33, dropk(pp=pp33,k=i,droplow=0))
}
#Half p-curves
if (khalf>0)
{
droplow.halfr=drophigh.halfr=c()
for (i in 0:(round(khalf/2)-1))
{
#Drop the lowest k studies in terms of respective overall test
droplow.halfr= c(droplow.halfr, dropk(pp=ppr.half,k=i,droplow=1))
drophigh.halfr=c(drophigh.halfr, dropk(pp=ppr.half,k=i,droplow=0))
} #End loop
}#End if that runs calculations only if khalf>0
#7.3 FUNCTION THAT DOES THE PLOT OF RESULTS
plotdrop=function(var,col)
{
k=length(var)
#Plot the dots
plot(0:(k-1),var,xlab="",ylab="",type="b",yaxt="n",xaxt="n",main="",
cex.main=1.15,ylim=c(0,1),col=col)
#Add marker in results with 0 drops
points(0,var[1],pch=19,cex=1.6)
#Red line at p=.05
abline(h=.05,col="red")
#Y-axis value labels
axis(2,c(.05,2:9/10),labels=c('.05','.2','.3','.4','.5','6','7','.8','.9'),las=1,cex.axis=1.5)
axis(1,c(0:(k-1)),las=1,cex.axis=1.4)
}
######################################################################################
# 9. Effect Estimation ###############################################################
######################################################################################
if (effect.estimation == TRUE){
# Define ci.to.t function
ci.to.t = function(TE, lower, upper, n){
z.to.d = function(z, n){
d = (2*z)/sqrt(n)
return(abs(d))
}
ci.to.p = function(est, lower, upper){
SE = (upper-lower)/(2*1.96)
z = abs(est/SE)
p = exp(-0.717*z - 0.416*z^2)
return(p)
}
d.to.t = function(d, n){
df = n-2
t = (d*sqrt(df))/2
return(t)
}
p = ci.to.p(TE, lower, upper)
z = abs(qnorm(p/2))
d = z.to.d(z, n)
t = d.to.t(d, n)
return(t)
}
#Function 13 - loss function
loss=function(t_obs,df_obs,d_est) {
#1.Convert all ts to the same sign (for justification see Supplement 5)
t_obs=abs(t_obs)
#2 Compute p-values
p_obs=2*(1-pt(t_obs,df=df_obs))
#3 Keep significant t-values and corresponding df.
t.sig=subset(t_obs,p_obs<.05)
df.sig=subset(df_obs,p_obs<.05)
#4.Compute non-centrality parameter implied by d_est and df_obs
#df+2 is total N.
#Becuase the noncentrality parameter for the student distribution is ncp=sqrt(n/2)*d,
#we add 2 to d.f. to get N, divide by 2 to get n, and by 2 again for ncp, so -->df+2/4
ncp_est=sqrt((df.sig+2)/4)*d_est
#5.Find critical t-value for p=.05 (two-sided)
#this is used below to compute power, it is a vector as different tests have different dfs
#and hence different critical values
tc=qt(.975,df.sig)
#4.Find power for ncp given tc, again, this is a vector of implied power, for ncp_est, for each test
power_est=1-pt(tc,df.sig,ncp_est)
#5.Compute pp-values
#5.1 First get the overall probability of a t>tobs, given ncp
p_larger=pt(t.sig,df=df.sig,ncp=ncp_est)
#5.2 Now, condition on p<.05
ppr=(p_larger-(1-power_est))/power_est #this is the pp-value for right-skew
#6. Compute the gap between the distribution of observed pp-values and a uniform distribution 0,1
KSD=ks.test(ppr,punif)$statistic #this is the D statistic outputted by the KS test against uniform
return(KSD)
}
if(missing(N)){
stop("If 'effect.estimation=TRUE', argument 'N' must be provided.")
}
if (length(N) != length(metaobject$TE)){
stop("N must be of same length as the number of studies contained in x.")
}
lower = metaobject$TE - (metaobject$seTE*1.96)
upper = metaobject$TE + (metaobject$seTE*1.96)
t_obs = ci.to.t(metaobject$TE, lower, upper, N)
df_obs = N-2
#Results will be stored in these vectors, create them first
loss.all=c()
di=c()
#Compute loss for effect sizes between d=c(dmin,dmax) in steps of .01
for (i in 0:((dmax-dmin)*100))
{
d=dmin+i/100 #effect size being considered
di=c(di,d) #add it to the vector (kind of silly, but kept for symmetry)
options(warn=-1) #turn off warning becuase R does not like its own pt() function!
loss.all=c(loss.all,loss(df_obs=df_obs,t_obs=t_obs,d_est=d))
#apply loss function so that effect size, store result
options(warn=0) #turn warnings back on
}
#find the effect leading to smallest loss in that set, that becomes the starting point in the optimize command
imin=match(min(loss.all),loss.all) #which i tested effect size lead to the overall minimum?
dstart=dmin+imin/100 #convert that i into a d.
#optimize around the global minimum
dhat=optimize(loss,c(dstart-.1,dstart+.1), df_obs=df_obs,t_obs=t_obs)
options(warn=-0)
#Plot results
plot(di,loss.all,xlab="Effect size\nCohen-d", ylab="Loss (D stat in KS test)",ylim=c(0,1), main="How well does each effect size fit? (lower is better)")
points(dhat$minimum,dhat$objective,pch=19,col="red",cex=2)
text(dhat$minimum,dhat$objective-.08,paste0("p-curve's estimate of effect size:\nd=",round(dhat$minimum,3)),col="red")
}
######################################################################################
# 10. Prepare Results for Return #####################################################
######################################################################################
# Get results
main.results = round(main.results, 3)
ktotal = round(main.results[1]) # Get the total number of inserted TEs
k.sign = round(main.results[2]) # Get the total number of significant TEs
k.025 = round(main.results[3]) # Get the number of p<0.25 TEs
skew.full.z = main.results[4] # Get the Z-score for the full curve skewness test
skew.full.p = main.results[5] # Get the p-value for the full curve skewness test
flat.full.z = main.results[6] # Get the Z-score for the full curve flatness test
flat.full.p = main.results[7] # Get the p-value for the full curve flatness test
skew.half.z = main.results[8] # Get the Z-score for the half curve skewness test
skew.half.p = main.results[9] # Get the p-value for the half curve skewness test
flat.half.z = main.results[10] # Get the Z-score for the half curve flatness test
flat.half.p = main.results[11] # Get the p-value for the half curve flatness test
skew.binomial.p = round(binomial[3], 3) # Get the skewness binomial p-value
flat.binomial.p = round(binomial[4], 3) # Get the flatness binomial p-value
# Make data.frame
skewness = c(skew.binomial.p, skew.full.z, skew.full.p, skew.half.z, skew.half.p)
flatness = c(flat.binomial.p, flat.full.z, flat.full.p, flat.half.z, flat.half.p)
colnames.df = c("pBinomial", "zFull", "pFull", "zHalf", "pHalf")
rownames.df = c("Right-skewness test", "Flatness test")
pcurveResults = rbind(skewness, flatness)
colnames(pcurveResults) = colnames.df
rownames(pcurveResults) = rownames.df
# Power results
power_results = round(power_results, 3)
powerEstimate = power_results[2]
powerLower = power_results[1]
powerUpper = power_results[3]
Power = as.data.frame(cbind(powerEstimate, powerLower, powerUpper))
rownames(Power) = ""
# Presence and absence of evidential value
# - If the half p-curve test is right-skewed with p<.05 or both the half and full test
# are right-skewed with p<.1, then p-curve analysis indicates the presence of evidential value
# - Evidential value is inadequate or absent if the 33% power test is p<.05 for the full p-curve
# or both the half p-curve and binomial 33% power test are p<.1
if (skew.half.p < 0.05 | (skew.half.p < 0.1 & skew.full.p < 0.1)){
presence.ev = "yes"
} else {
presence.ev = "no"
}
if (flat.full.p < 0.05 | (flat.half.p < 0.1 & flat.binomial.p < 0.1)){
absence.ev = "yes"
} else {
absence.ev = "no"
}
# Plot Data
PlotData = round(table_figure, 3)
# Input Data
table_calc[,1] = NULL
colnames(table_calc) = c("p", "ppSkewFull", "ppSkewHalf", "ppFlatFull", "ppFlatHalf", "zSkewFull", "zSkewHalf",
"zFlatFull", "zFlatHalf")
Input = cbind(metaobject$TE, round(table_calc,3))
rownames(Input) = paste(1:length(metaobject$TE), metaobject$studlab)
colnames(Input)[1] = "TE"
if (effect.estimation==TRUE){
dEstimate = round(dhat$minimum, 3)
return.list = list("pcurveResults" = pcurveResults,
"Power" = Power,
"PlotData" = PlotData,
"Input" = Input,
"EvidencePresent" = presence.ev,
"EvidenceAbsent" = absence.ev,
"kInput" = ktot,
"kAnalyzed" = k.sign,
"kp0.25" = k.025,
"dEstimate" = dEstimate,
"I2" = metaobject$I2,
"class.meta.object" = class(metaobject)[1])
class(return.list) = c("pcurve", "effect.estimation")
} else {
return.list = list("pcurveResults" = pcurveResults,
"Power" = Power,
"PlotData" = PlotData,
"Input" = Input,
"EvidencePresent" = presence.ev,
"EvidenceAbsent" = absence.ev,
"kInput" = ktot,
"kAnalyzed" = k.sign,
"kp0.25" = k.025,
"I2" = metaobject$I2,
"class.meta.object" = class(metaobject)[1])
class(return.list) = c("pcurve", "no.effect.estimation")
}
cat(" ", "\n")
invisible(return.list)
return.list
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.