Nothing
reitsma <- function(data, ...) UseMethod("reitsma")
reitsma.default <-
function(data = NULL, subset=NULL, formula = NULL,
TP="TP", FN="FN", FP="FP", TN="TN",
alphasens = 1, alphafpr = 1,
correction = 0.5, correction.control = "all",
method = "reml",
control = list(), ...)
{
call <- match.call()
mcall <- match.call(expand.dots = FALSE)
stopifnot(is.numeric(correction), 0 <= correction,
correction.control %in% c("all", "single", "none"),
0 <= alphasens, alphasens <= 2, 0 <= alphafpr, alphafpr <= 2,
method %in% c("fixed", "ml", "reml", "mm", "vc"),
is.numeric(TP) | (is.character(TP) & length(TP) == 1),
is.numeric(FP) | (is.character(FP) & length(FP) == 1),
is.numeric(TN) | (is.character(TN) & length(TN) == 1),
is.numeric(FN) | (is.character(FN) & length(FN) == 1))
if(is.null(formula)){formula <- cbind(tsens,tfpr)~1}
if(!is(formula, "formula")){stop("formula must be of class formula")}
if(!formula[2] == (cbind(tsens,tfpr)~1)[2]){stop("The left hand side of formula must be cbind(tsens, tfpr)")}
if(!is.null(data) & is.character(c(TP,FP,TN,FN))){
X <- as.data.frame(data)
origdata <- data
TP <- getElement(X,TP)
FN <- getElement(X,FN)
FP <- getElement(X,FP)
TN <- getElement(X,TN)
}
if(is.null(data) & !is.character(c(TP,FP,TN,FN))){
origdata <- data.frame(TP = TP, FN = FN, FP = FP, TN = TN)
}
freqdata <- cbind(TP,FN,FP,TN)
checkdata(freqdata)
N <- length(TP)
## apply continuity correction to _all_ studies if one contains zero
if(correction.control == "all"){if(any(c(TP,FN,FP,TN) == 0))
{TP <- TP + correction;
FN <- FN + correction;
FP <- FP + correction;
TN <- TN + correction}}
if(correction.control == "single"){
correction = ((((TP == 0)|(FN == 0))|(FP == 0))| (TN == 0)) *
correction
TP <- correction + TP
FN <- correction + FN
FP <- correction + FP
TN <- correction + TN
}
tpi <- TP
fni <- FN
fpi <- FP
tni <- TN
freqdata1 <- cbind(tpi,fni,fpi,tni)
number.of.pos <- TP + FN
number.of.neg <- FP + TN
sens<-TP/number.of.pos
fpr <- FP/number.of.neg
senstrafo <- function(x){return(talpha(alphasens)$linkfun(x))}
fprtrafo <- function(x){return(talpha(alphafpr)$linkfun(x))}
sensinvtrafo <- function(x){return(talpha(alphasens)$linkinv(x))}
fprinvtrafo <- function(x){return(talpha(alphafpr)$linkinv(x))}
trafo.sens <- senstrafo(sens)
trafo.fpr <- fprtrafo(fpr)
var.trafo.sens <- (sens*(1-sens)/number.of.pos) *
(jacobiantrafo(alphasens, sens)^2)
var.trafo.fpr <- (fpr*(1-fpr)/number.of.neg) *
(jacobiantrafo(alphafpr, fpr)^2)
# preparations for call to mvmeta.fit
mvmeta_S <- cbind(var.trafo.sens, rep(0,N), var.trafo.fpr)
mvmeta_data <- data.frame(X,tsens = trafo.sens, tfpr = trafo.fpr)
mvmeta_X <- model.matrix(formula,data=mvmeta_data)
mvmeta_y <- as.matrix(data.frame(tsens = trafo.sens, tfpr = trafo.fpr))
fit <- mvmeta.fit(mvmeta_X, mvmeta_y, mvmeta_S, method = method, control = control)
fit$call <- call
fit$formula <- formula
fit$control <- control
fit$method <- method
fit$alphasens <- alphasens
fit$alphafpr <- alphafpr
fit$data <- origdata
fit$freqdata <- as.data.frame(freqdata)
fit$freqdata1 <- as.data.frame(freqdata1)
class(fit) <- "reitsma"
## modify log likelihood, especially add jacobian
attr(fit$logLik, "df") <- 2*ncol(mvmeta_X) + 3
fit$logLik <- fit$logLik + sum(log(jacobiantrafo(alphasens,sens)) +
log(jacobiantrafo(alphafpr,fpr)))
fit$nobs <- 2*N
return(fit)
class(output) <- "reitsma"
return(output)
}
print.reitsma <- function (x, digits = 4, ...)
{
methodname <- c("reml", "ml", "fixed")
methodlabel <- c("REML", "ML", "Fixed")
cat("Call: ", paste(deparse(x$call), sep = "\n", collapse = "\n"),
"\n\n", sep = "")
cat("Fixed-effects coefficients:", "\n", sep = "")
table <- formatC(x$coefficients, digits = digits, format = "f")
table <- gsub("NA"," -",table)
print(table, quote = FALSE, right = TRUE, print.gap = 2)
cat("\n")
cat(x$dim$m, " studies, ",
x$df$fixed, " fixed and ", x$df$random, " random-effects parameters",
"\n", sep = "")
table <- c(x$logLik, AIC(x), BIC(x))
names(table) <- c("logLik", "AIC", "BIC")
table <- formatC(table, digits = digits, format = "f")
print(table, quote = FALSE, right = TRUE, print.gap = 2)
cat("\n")
}
summary.reitsma <-
function (object, level = 0.95, sroc.type = "ruttergatsonis", ...)
{
stopifnot(sroc.type %in% c("naive", "ruttergatsonis"))
ci.level <- level
if (ci.level <= 0 || ci.level >= 1)
stop("'ci.level' must be within 0 and 1")
coef <- object$coefficients
vcov <- object$vcov
dim <- object$dim
Psi <- object$Psi
lab <- object$lab
coef <- as.numeric(coef)
coef.se <- sqrt(diag(vcov))
zval <- coef/coef.se
zvalci <- qnorm((1 - ci.level)/2, lower.tail = FALSE)
pvalue <- 2 * (1 - pnorm(abs(zval)))
ci.lb <- coef - zvalci * coef.se
ci.ub <- coef + zvalci * coef.se
cilab <- paste(signif(ci.level, 2) * 100, "%ci.", c("lb",
"ub"), sep = "")
tabfixed <- cbind(coef, coef.se, zval, pvalue, ci.lb, ci.ub)
dimnames(tabfixed) <- list(if (dim$k > 1L) colnames(vcov) else lab$p,
c("Estimate", "Std. Error", "z", "Pr(>|z|)", cilab))
if(nrow(tabfixed) == 2){
tabfixed <- rbind(tabfixed,tabfixed)
tabfixed[3,] <- inv.trafo(object$alphasens, tabfixed[3,])
tabfixed[4,] <- inv.trafo(object$alphafpr, tabfixed[4,])
tabfixed[3:4,c(2,3,4)] <- NA
rownames(tabfixed)[3:4] <- c("sensitivity","false pos. rate")
}
corFixed <- vcov/outer(coef.se, coef.se)
ran.sd <- sqrt(diag(Psi))
corRandom <- Psi/outer(ran.sd, ran.sd)
# qstat <- unclass(qtest(object))
tau_a <- (tabfixed[1,2])^2
tau_b <- (tabfixed[2,2])^2
cor_sq <- (corRandom[1,2])^2
theta_a <- (tabfixed[1,1])
theta_b <- 0-(tabfixed[2,1])
K <- 1/as.numeric(dim$m)
k <- as.numeric(dim$m)
Nia <- 1/(object$freqdata$TP + object$freqdata$FP)
Nib <- 1/(object$freqdata$TN + object$freqdata$FN)
NiA <- sum(Nia)
NiB <- sum(Nib)
tau <- (1-cor_sq)*tau_a*tau_b
em <- ((exp((tau_a/2)+theta_a)+exp((tau_a/2)-theta_a)+2)*K*NiA)*((exp((tau_b/2)+theta_b)+exp((tau_b/2)-theta_b)+2)*K*NiB)
i2_z <- (tau^0.5)/((em^0.5)+(tau^0.5))
theta_i <- log((object$freqdata1$tpi/(object$freqdata1$tpi+object$freqdata1$fni)))/log((object$freqdata1$fpi/(object$freqdata1$fpi+object$freqdata1$tni)))
var_i <- ((1/(object$freqdata1$tpi+object$freqdata1$fni))*((1-(object$freqdata1$tpi/(object$freqdata1$tpi+object$freqdata1$fni)))/((object$freqdata1$tpi/(object$freqdata1$tpi+object$freqdata1$fni))*(log((object$freqdata1$tpi/(object$freqdata1$tpi+object$freqdata1$fni)))^2)))) + ((1/(object$freqdata1$fpi+object$freqdata1$tni))*((1-(object$freqdata1$fpi/(object$freqdata1$fpi+object$freqdata1$tni)))/((object$freqdata1$fpi/(object$freqdata1$fpi+object$freqdata1$tni))*(log((object$freqdata1$fpi/(object$freqdata1$fpi+object$freqdata1$tni)))^2))))
weight_i <- 1/var_i
var_i0 <- (((1-(object$freqdata1$tpi/(object$freqdata1$tpi+object$freqdata1$fni)))/((object$freqdata1$tpi/(object$freqdata1$tpi+object$freqdata1$fni))*(log((object$freqdata1$tpi/(object$freqdata1$tpi+object$freqdata1$fni)))^2))) + ((1-(object$freqdata1$fpi/(object$freqdata1$fpi+object$freqdata1$tni)))/((object$freqdata1$fpi/(object$freqdata1$fpi+object$freqdata1$tni))*(log((object$freqdata1$fpi/(object$freqdata1$fpi+object$freqdata1$tni)))^2))))
weight_i0 <- 1/var_i0
sq_weight_i <- weight_i*weight_i
sq_weight_i0 <- weight_i0*weight_i0
weighted_theta_i <- weight_i*log(theta_i)
sum_var_i <- sum(var_i)
sum_var_i0 <- sum(var_i0)
sum_weight_i <- sum(weight_i)
sum_weight_i0 <- sum(weight_i0)
sum_sq_weight_i <- sum(sq_weight_i)
sum_sq_weight_i0 <- sum(sq_weight_i0)
sum_fe_estimator <- sum(weighted_theta_i)/sum_weight_i
chisq_i <- weight_i*((log(theta_i)-sum_fe_estimator)^2)
chi_sq <- sum(chisq_i)
tau02 <- (chi_sq-(k-1)) / (sum_weight_i - (sum_sq_weight_i/sum_weight_i))
tau2 <- max(c(0,tau02))
var1 <- ((k-1)*sum_weight_i)/ ((sum_weight_i*sum_weight_i)-sum_sq_weight_i)
var2 <- (k/sum_weight_i)
var3 <- (sum_var_i/k)
var10 <- ((k-1)*sum_weight_i0)/((sum_weight_i0*sum_weight_i0)-sum_sq_weight_i0)
var20 <- (k/sum_weight_i0)
var30 <- (sum_var_i0/k)
i2_1 <- tau2/(tau2+var1)
i2_2 <- tau2/(tau2+var2)
i2_3 <- tau2/(tau2+var3)
i2_10 <- tau2/(tau2+var10)
i2_20 <- tau2/(tau2+var20)
i2_30 <- tau2/(tau2+var30)
i2 <- data.frame("Zhou"=i2_z,"HollingUnadjusted1"=i2_1,"HollingUnadjusted2"=i2_2,"HollingUnadjusted3"=i2_3,"HollingAdjusted1"=i2_10,"HollingAdjusted2"=i2_20,"HollingAdjusted3"=i2_30)
if(attr(object$logLik,"df")== 5){
AUC <- AUC.reitsma(object, sroc.type = sroc.type)
coef_hsroc <- calc_hsroc_coef(object)
}else{
AUC <- NULL
coef_hsroc <- NULL}
keep <- match(c("vcov", "Psi", "df.res", "rank", "logLik",
"converged", "dim", "df", "lab", "na.action", "call",
"terms", "method", "alphasens", "alphafpr"), names(object), 0L)
out <- c(list(coefficients = tabfixed), object[keep], list(AIC = AIC(object),
BIC = BIC(object),
corFixed = corFixed,
corRandom = corRandom,
# qstat = qstat,
ci.level = ci.level,
AUC = AUC,
i2=i2))
class(out) <- "summary.reitsma"
return(out)
}
print.summary.reitsma <- function (x, digits = 3, ...){
methodname <- c("reml", "ml", "fixed")
methodlabel <- c("REML", "ML", "Fixed")
cat("Call: ", paste(deparse(x$call), sep = "\n", collapse = "\n"),
"\n\n", sep = "")
cat("Bivariate diagnostic ",
ifelse(x$method == "fixed", "fixed", "random"),
"-effects meta-", ifelse(x$dim$p > 2,"regression", "analysis"), "\n", sep = "")
if (x$method != "fixed") {
cat("Estimation method: ", methodlabel[which(x$method ==
methodname)], "\n", sep = "")
# cat("Variance-covariance matrix Psi: ", "unstructured",
# "\n", sep = "")
}
cat("\n")
cat("Fixed-effects coefficients", "\n", sep = "")
signif <- symnum(x$coefficients[, "Pr(>|z|)"], corr = FALSE,
na = FALSE, cutpoints = c(0, 0.001, 0.01, 0.05, 0.1,
1),
symbols = c("***", "**", "*", ".", " "))
tabletot <- formatC(x$coefficients, digits = digits, format = "f")
tabletot <- gsub("NA"," -",tabletot)
tabletot <- cbind(tabletot, signif)
colnames(tabletot)[7] <- ""
print(tabletot, quote = FALSE, right = TRUE, print.gap = 1)
cat("---\nSignif. codes: ", attr(signif, "legend"), "\n\n")
if (!x$method == "fixed") {
cat("Variance components: between-studies Std. Dev and correlation matrix",
"\n", sep = "")
corRan <- x$corRandom
corRan[upper.tri(x$corRan)] <- NA
table <- cbind(`Std. Dev` = sqrt(diag(x$Psi)), corRan)
if (x$dim$k == 1L)
rownames(table) <- ""
table <- formatC(table, digits = digits, format = "f")
table[grep("NA", table)] <- "."
print(table, quote = FALSE, right = TRUE, na.print = "",
print.gap = 1)
cat("\n")
}
table <- c(x$logLik, x$AIC, x$BIC)
names(table) <- c("logLik", "AIC", "BIC")
table <- formatC(table, digits = digits, format = "f")
print(table, quote = FALSE, right = TRUE, print.gap = 1)
cat("\n")
if(!is.null(x$AUC)){
cat(c("AUC: ",as.character(round(x$AUC$AUC,3))))
cat("\n")
cat(c("Partial AUC (restricted to observed FPRs and normalized): ",as.character(round(x$AUC$pAUC,3))))
cat("\n")
}
if(!is.null(x$coef_hsroc)){
cat("\n")
cat("HSROC parameters \n")
table <- as.numeric(x$coef_hsroc)
names(table) <- names(x$coef_hsroc)
table <- formatC(table, digits = digits, format = "f")
print(table, quote = FALSE, right = TRUE, print.gap = 1)
cat("\n")
}
cat("\n")
cat("I2 estimates \n")
cat("Zhou and Dendukuri approach: ", round(((x$i2$Zhou)*100),1) ,"% \n")
cat("Holling sample size unadjusted approaches: ", round((min(c(x$i2$HollingUnadjusted1,x$i2$HollingUnadjusted2,x$i2$HollingUnadjusted3))*100),1),"-",round((max(c(x$i2$HollingUnadjusted1,x$i2$HollingUnadjusted2,x$i2$HollingUnadjusted3))*100),1) ,"% \n")
cat("Holling sample size adjusted approaches: ", round((min(c(x$i2$HollingAdjusted1,x$i2$HollingAdjusted2,x$i2$HollingAdjusted3))*100),1),"-",round((max(c(x$i2$HollingAdjusted1,x$i2$HollingAdjusted2,x$i2$HollingAdjusted3))*100),1) ,"%")
}
vcov.reitsma <- function(object, ...){object$vcov}
logLik.reitsma <- function(object, ...){object$logLik}
sroc.reitsma <- function(fit, fpr = 1:99/100, type = "ruttergatsonis",
return_function = FALSE, ...){
stopifnot(is.logical(return_function))
stopifnot(type %in% c("ruttergatsonis", "naive"))
if(type == "ruttergatsonis"){
coef_hsroc <- calc_hsroc_coef(fit)
}
if(type == "naive"){
f <- function(x){return(sroc2(fit, fpr=x, ...)[,2])}
if(!return_function){
return(cbind(fpr, f(fpr)))
}else{
return(f)
}
}
if(type == "ruttergatsonis"){
Lambda <- coef_hsroc$Lambda
Beta <- coef_hsroc$beta
f <- function(x){
return(inv.trafo(fit$alphasens, (Lambda*exp(-Beta/2) + exp(-Beta)*trafo(fit$alphafpr, x))))
}
if(!return_function){
return(cbind(fpr, f(fpr)))
}else{
return(f)
}
}
}
mcsroc.reitsma <- function(fit, fpr = 1:99/100, replications = 10000, lambda = 100, ...){
mcsroc(fit, replications = replications, lambda = lambda, ...)
}
ROCellipse.reitsma <- function(x, level = 0.95, add = FALSE, pch = 1, ...){
ROC.ellipse2(x, conf.level = level, add = add, pch = pch, ...)
}
crosshair.reitsma <- function(x, level = 0.95, length = 0.1, pch = 1, ...){
s <- summary(x, level = level)
ch <- rbind(inv.trafo(x$alphasens,s$coefficients["tsens.(Intercept)",5:6]),
inv.trafo(x$alphafpr,s$coefficients["tfpr.(Intercept)",5:6]))
pe <- c(inv.trafo(x$alphafpr,s$coefficients["tfpr.(Intercept)",1]),
inv.trafo(x$alphasens,s$coefficients["tsens.(Intercept)",1]))
# if(!add){
# plot(matrix(pe, ncol = 2), pch = pch, xlim = xlim, ylim = ylim,
# ylab = "Sensitivity", xlab = "False Positive Rate", ...)
# }else{
points(matrix(pe, ncol = 2), pch = pch, ...)
# }
arrows(ch[2,1], pe[2], x1 =ch[2,2], length = length, angle = 90, code = 3, ...)
arrows(pe[1], ch[1,1], y1 =ch[1,2], length = length, angle = 90, code = 3, ...)
return(invisible(NULL))
}
plot.reitsma <- function(x, extrapolate = FALSE, plotsumm = TRUE, level = 0.95,
ylim = c(0,1), xlim = c(0,1), pch = 1,
sroclty = 1, sroclwd = 1,
predict = FALSE, predlty = 3, predlwd = 1,
type = "ruttergatsonis",
...)
{
plot(c(2,2), ylim = ylim, xlim = xlim,
xlab = "False Positive Rate", ylab = "Sensitivity", ...)
if(length(coef(x)) == 2){
FP <- x$freqdata$FP
negatives <- FP + x$freqdata$TN
FPR <- FP/negatives
if(extrapolate){bound = c(0,1)}
if(!extrapolate){bound = c(min(FPR), max(FPR))}
srocmat <- sroc(x, type = type)
lines(srocmat[cut(srocmat[,1],bound, "withinbound") == "withinbound",],
lty = sroclty, lwd = sroclwd)
}else{
warning("Not plotting any SROC for meta-regression")
}
if(plotsumm){ROCellipse(x, level = level, add = TRUE, pch = pch, ...)}
if(predict){
alpha.sens <- x$alphasens
alpha.fpr <- x$alphafpr
mu <- x$coefficients["(Intercept)",]
Sigma <- x$Psi + vcov(x)
talphaellipse <- ellipse(Sigma, centre = mu, level = level)
predellipse <- matrix(0, ncol = 2, nrow = nrow(talphaellipse))
predellipse[,1] <- inv.trafo(alpha.fpr, talphaellipse[,2])
predellipse[,2] <- inv.trafo(alpha.sens, talphaellipse[,1])
lines(predellipse, lty = predlty, lwd = predlwd)
}
return(invisible(NULL))
}
confint.reitsma <- function (object, parm, level = 0.95, ...)
{
cf <- as.numeric(coef(object))
if(missing(parm)){parm <- colnames(vcov(object))}
a <- (1 - level)/2
a <- c(a, 1 - a)
pct <- paste(format(100 * a, trim = TRUE, scientific = FALSE, digits = 3),
"%")
fac <- qnorm(a)
ci <- array(NA, dim = c(length(parm), 2L), dimnames = list(parm,
pct))
ses <- sqrt(diag(vcov(object)))
ci[] <- cf + ses %o% fac
ci
}
AUC.reitsma <- function(x, fpr = 1:99/100,
sroc.type = "ruttergatsonis",
...){
stopifnot(sroc.type %in% c("naive", "ruttergatsonis"))
if(!attr(x$logLik,"df") == 5){stop("AUC can not be calculated for meta-regression")}
rsroc <- sroc.reitsma(x, type = sroc.type,
return_function = TRUE)
AUC <- AUC.default(rsroc, fpr = fpr, ...)
obsfprrange <- range(fpr(x$freqdata))
obsfprrange[1] <- max(0.01,obsfprrange[1])
obsfprrange[2] <- min(0.99,obsfprrange[2])
obsfpr <- seq(from = obsfprrange[1], to = obsfprrange[2], length.out = 99)
pAUC <- AUC.default(rsroc, fpr = obsfpr, ...)
names(pAUC) <- c("pAUC")
res <- c(AUC,pAUC)
attr(res, "sroc.type") <- sroc.type
return(res)
}
calc_hsroc_coef <- function(fit){
coef <- fit$coefficients
coef <- as.numeric(coef)
Psi <- fit$Psi
ran.sd <- sqrt(diag(Psi))
if(attr(fit$logLik,"df")== 5){
## HSROC parameters (formulae from Harbord et al)
Theta <- 0.5*(sqrt(ran.sd[2]/ran.sd[1])*coef[1] + sqrt(ran.sd[1]/ran.sd[2])*coef[2]) ##coef[2] is the fpr, so need to change sign!
Lambda <- sqrt(ran.sd[2]/ran.sd[1])*coef[1] - sqrt(ran.sd[1]/ran.sd[2])*coef[2] ##coef[2] is the fpr, so need to change sign!
sigma2theta <- 0.5*(ran.sd[1]*ran.sd[2] + Psi[1,2]) ##coef[2] is the fpr, so need to change sign of Psi[1,2] as well!
sigma2alpha <- 2*(ran.sd[1]*ran.sd[2] - Psi[1,2]) ##coef[2] is the fpr, so need to change sign of Psi[1,2] as well!
beta <- log(ran.sd[2]/ran.sd[1])
coef_hsroc <- list(Theta = Theta,
Lambda = Lambda,
beta = beta,
sigma2theta = sigma2theta,
sigma2alpha = sigma2alpha)
coef_hsroc <- lapply(coef_hsroc, function(x){attr(x, "names") <- NULL; x})
}else{
coef_hsroc = NULL
}
if(is.null(coef_hsroc)){
warning("Can only compute coefficients for SROC curves without covariates. Returning NULL.")
}
return(coef_hsroc)
}
anova.reitsma <- function(object, fit2, ...){
fit1 <- object
if(!(is(fit1, "reitsma") & is(fit2, "reitsma"))){
stop("This function is only for models of class 'reitsma'.")
}
if(!(fit1$method == "ml" & fit2$method == "ml")){
stop("Both models have to be fitted with maximum likelihood. Set method = 'ml' in reitsma.")
if(!(identical(fit1$data, fit2$data))){
stop("Both models have to be fitted to the same data set.")
}
if(!(fit1$nobs == fit2$nobs)){
stop("Both models have to be fitted to the same number of primary studies. Make sure missing data did not lead to the omission of studies.")
}
}
df1 <- fit1$df$fixed
df2 <- fit2$df$fixed
df_diff <- df1 - df2
if(df_diff == 0){stop("Models have the same degrees of freedom, so they cannot be nested!")}
logLik_diff <- fit1$logLik - fit2$logLik
if((df_diff > 0 & logLik_diff < 0) | (df_diff < 0 & logLik_diff > 0)){
stop("log-likelihoods not consistent, models probably not nested")}
out <- c(2*abs(logLik_diff), abs(df_diff),
pchisq(2*abs(logLik_diff), abs(df_diff), lower.tail = FALSE))
names(out) <- c("Chi-squared", "Df", "Pr(>Chi-squared)")
out <- list(statistic = out, fit1 = fit1, fit2 = fit2)
class(out) <- "anova.reitsma"
out
}
print.anova.reitsma <- function(x, digits = 4, ...){
cat("Likelihood-ratio test \n")
cat("Model 1: ")
print(x$fit1$formula)
cat("Model 2: ")
print(x$fit2$formula)
cat("\n")
ret <- data.frame(ChiSquared = x$statistic[1],
Df = x$statistic[2],
x$statistic[3])
names(ret)[3] <- "Pr(>ChiSquared)"
rownames(ret) <- ""
printCoefmat(ret,
digits = digits,
tst.ind = 1,
zap.ind = 2,
P.values = TRUE,
has.Pvalue = TRUE,
cs.ind = NULL)
return(invisible(x))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.