Nothing
#' Compute a test of trend in prevalences based on a likelihood-ratio statistic
#'
#' This function takes a series of point estimates and their associated standard errors and
#' computes the p-value for the test of a monotone decrease in the
#' population prevalences (in sequence order).
#' The p-value for a monotone increase is
#' also reported. An optional plot of the estimates and the null distribution of the test statistics is provided.
#' More formally, let the \eqn{K} population prevalences in sequence order be \eqn{p_1, \ldots, p_K}.
#' We test the null hypothesis:\cr
#' \deqn{H_0 : p_1 = \ldots = p_K}
#' vs
#' \deqn{H_1 : p_1 \ge p_2 \ldots \ge p_K}
#' with at least one equality strict. The alternatie hypothesis is for a monotone decreasing trend.
#' A likelihood ratio statistic for this test has
#' been derived (Bartholomew 1959).
#' The null distribution of the likelihood ratio statistic is very complex
#' but can be determined by a simple Monte Carlo process.\cr
#' Alternatively, we can test the null hypothesis:\cr
#' \deqn{H_0 : p_1 \ge p_2 \ldots \ge p_K}
#' vs
#' \deqn{H_1 : \overline{H_0}}
#' The null distribution of the likelihood ratio statistic is very complex
#' but can be determined by a simple Monte Carlo process.
#' In both cases we also test for:\cr
#' \deqn{H : p_1 \le p_2 \ldots \le p_K}
#' that is, a monotonically increasing trend.
#' The function requires the isotone library.
#'
#' @aliases LRT.trend
#' @aliases LRT.trend.null
#' @aliases LRT.trend.test
#' @param data A two row matrix or data.frame of prevalence estimates and
#' their standard errors. The first row is the prevalence estimates and the
#' second are the standard errors. The column are the comparison groups in the
#' order (e.g., time) there are to be assessed. The row names of \code{data}
#' should be "estimate" and "sigma". This is
#' @param variables A character vector of column names it select from \code{data}.
#' @param null A character string indicating the null hypothesis to use. The value \code{"monotone"} uses the various
#' monotone hypotheses as the nulls. If not \code{"monotone"}, the null is chosen to be that of equality of the means
#' over all periods.
#' @param confidence.level The confidence level for the confidence intervals. The default is 0.95 for 95\%.
#' @param number.of.bootstrap.samples The number of Monte Carlo draws to
#' determine the null distribution of the likelihood ratio statistic.
#' @param plot A character vector of choices, a subset of \code{estimates}, \code{distributions}.
#' If \code{estimates} is given then a plot of the estimates and nominal 95\% confidence bands (as error bars) is produced.
#' If \code{distributions} is given then a plot is produced of the null distributions of
#' the likelihood
#' ratio statistic with the observed likelihood ratio statistics plotted as a vertical dashed line.
#' @param seed The value of the random number seed. Preset by default to allow reproducibility.
#' @return
#' A list with components
#' \itemize{ \item\code{pvalue.increasing}: The p-value for the test of a monotone increase in population prevalence.
#' \item\code{pvalue.decreasing}: The p-value for the test of a monotone decrease in population prevalence.
#' \item\code{L}: The value of the likelihood-ratio statistic.
#' \item\code{x}: The passed vector of prevalence estimates in the order (e.g., time).
#' \item\code{sigma} The passed vector of standard error estimates corresponding to \code{x}.
#' }
#'
#' @author Mark S. Handcock
#' @references Bartholomew, D. J. (1959). A test of homogeneity for ordered alternatives. Biometrika 46 36-48.
#' @keywords survey manip
#' @examples
#'
#' d <- t(data.frame(estimate=c(0.16,0.15,0.3), sigma=c(0.04,0.04,0.1)))
#' colnames(d) <- c("time_1","time_2","time_3")
#' LRT.trend.test(d,number.of.bootstrap.samples=1000)
#' @export
LRT.trend.test <- function(data,variables=colnames(data),null="monotone",confidence.level=0.95,
number.of.bootstrap.samples=5000,plot=NULL,seed=1){
LRT.fn <- ifelse(null=="monotone",LRT.trend.null,LRT.trend)
out <- LRT.fn(x=as.numeric(data["estimate",variables]),
sigma=as.numeric(data["sigma",variables]),
number.of.bootstrap.samples=number.of.bootstrap.samples,
confidence.level=confidence.level,
plot=plot,seed=seed)
invisible(out)
}
# Distribution
#' @export
LRT.trend <- function(x,sigma,number.of.bootstrap.samples=5000,confidence.level=0.95,plot=NULL,seed=1,smooth=0.35){
# require(isotone, quietly=TRUE, warn.conflicts=FALSE)
# if(plot){
# suppressMessages(require(locfit, quietly = TRUE))
# }
a <- 1/sigma^2
xbar <- sum(a*x)/sum(a)
k <- length(x)
set.seed(seed)
r <- matrix(stats::rnorm(k*number.of.bootstrap.samples),ncol=k)
r <- sweep(r,2,sigma,"*");r <- sweep(r,2,xbar,"+")
L <- t(apply(r,1,LRT.value.trend,sigma=sigma))
obsL <- LRT.value.trend(x,sigma)
options("scipen"=100)
cat(sprintf("In a test of the null hypothesis of equal proportions against the alternative that at least two are unequal, the p-value is %s.\n",format.pval(mean(obsL[3] <= L[,3]),eps=0.0001)))
if(mean(obsL[3] <= L[,3]) < (1-confidence.level)){
cat(sprintf("We reject the null hypothesis of equal proportions in favor of the alternative that at least two are unequal (at the %s%% level).\n\n",format(100*(1-confidence.level))))
}else{
cat(sprintf("We do not reject the null hypothesis of equal proportions (at the %s%% level).\n\n",format(100*(1-confidence.level))))
}
if(all(diff(x) >= 0) | all(diff(x) <= 0)){
cat(sprintf("In a test of the null hypothesis of equal proportions against the alternative of a monotone trend (either up or down), the p-value is %s.\n",format.pval(mean(obsL[4] <= L[,4]),eps=0.0001)))
if(mean(obsL[4] <= L[,4]) < (1-confidence.level)){
cat(sprintf("We reject the null hypothesis of equal proportions in favor of the alternative of a monotone trend (at the %s%% level).\n\n",format(100*(1-confidence.level))))
}else{
cat(sprintf("We do not reject the null hypothesis of equal proportions (at the %s%% level).\n\n",format(100*(1-confidence.level))))
}
}else{
cat(sprintf("We do not accept the hypothesis of a monotone trend (at the %s%% level).\n\n",format(100*(1-confidence.level))))
}
if(all(diff(x) >= 0)){
cat(sprintf("In a test of the null hypothesis of equal proportions against the alternative of an increasing trend, the p-value is %s.\n",format.pval(mean(obsL[1] <= L[,1]),eps=0.0001)))
if(mean(obsL[1] <= L[,1]) < (1-confidence.level)){
cat(sprintf("We reject the null hypothesis of equal proportions in favor of the alternative of an increasing trend (at the %s%% level).\n\n",format(100*(1-confidence.level))))
}else{
cat(sprintf("We do not reject the null hypothesis of equal proportions (at the %s%% level).\n\n",format(100*(1-confidence.level))))
}
}else{
cat(sprintf("We do not accept the hypothesis of an increasing trend (at the %s%% level).\n\n",format(100*(1-confidence.level))))
}
if(all(diff(x) <= 0)){
cat(sprintf("In a test of the null hypothesis of equal proportions against the alternative of an decreasing trend, the p-value is %s.\n",format.pval(mean(obsL[2] <= L[,2]),eps=0.0001)))
if(mean(obsL[2] <= L[,2]) < (1-confidence.level)){
cat(sprintf("We reject the null hypothesis of equal proportions in favor of the alternative of a decreasing trend (at the %s%% level).\n\n",format(100*(1-confidence.level))))
}else{
cat(sprintf("We do not reject the null hypothesis of equal proportions (at the %s%% level).\n\n",format(100*(1-confidence.level))))
}
}else{
cat(sprintf("We do not accept the hypothesis of a decreasing trend (at the %s%% level).\n\n",format(100*(1-confidence.level))))
}
#
if("estimates" %in% plot){
yminus <- x - 1.96*sigma
yplus <- x + 1.96*sigma
errorbar(x=seq_along(x),y=x,yminus=yminus,
yplus=yplus,xlab="sequence",ylab="estimate",
main="Trend of estimates")
}
if("distributions" %in% plot){
binn <- 500
maxl <- stats::quantile(L[,3],0.95)
gpdf=bgk_kde(L[,3],n=2^(ceiling(maxl/log(2))),MIN=0,MAX=1.1*max(L[,3]))
maxl <- max(1.1*obsL[3],maxl)
r <- seq(0, maxl, length = binn + 1)[-1] - 0.5*maxl/binn
gpdf <- stats::spline(x=gpdf[1,],y=gpdf[2,],xout=r,method="natural")$y
gpdf[r>max(L[,3])] <- 0
gpdf[gpdf<0] <- 0
scalef <- binn/sum(gpdf)
gpdf <- gpdf * scalef
if(mean(obsL[3] <= L[,3]) < (1-confidence.level)){
mainm <- paste("Null Distribution of the Overall Test Statistic",
sprintf("\nWe reject the null hypothesis of equal proportions (at the %s%% level).",format(100*(1-confidence.level))))
}else{
mainm <- paste("Null Distribution of the Overall Test Statistic",
sprintf("\nWe do not reject the null hypothesis of equal proportions (at the %s%% level).",format(100*(1-confidence.level))))
}
plot(x=r,y=gpdf,xlim=c(0,maxl), type="l",
ylab="Density",
xlab="Likelihood Ratio Statistic",main=mainm,cex.main=1,sub="The vertical line is the observed likelihood ratio statistic")
graphics::abline(v=obsL[3],lty=2)
maxl <- stats::quantile(L[,1],0.95)
gpdf=bgk_kde(L[,1],n=2^(ceiling(maxl/log(2))),MIN=0,MAX=1.1*max(L[,1]))
maxl <- max(1.1*obsL[1],maxl)
r <- seq(0, maxl, length = binn + 1)[-1] - 0.5*maxl/binn
gpdf <- stats::spline(x=gpdf[1,],y=gpdf[2,],xout=r,method="natural")$y
gpdf[r>max(L[,1])] <- 0
gpdf[gpdf<0] <- 0
scalef <- binn/sum(gpdf)
gpdf <- gpdf * scalef
if(all(diff(x) >= 0)){
if(mean(obsL[1] <= L[,1]) < (1-confidence.level)){
mainm <- paste("Null Distribution of the Increasing Test Statistic",
sprintf("\nWe reject the null hypothesis of equal proportions (at the %s%% level).",format(100*(1-confidence.level))))
}else{
mainm <- paste("Null Distribution of the Increasing Test Statistic",
sprintf("\nWe do not reject the null hypothesis of equal proportions (at the %s%% level).",format(100*(1-confidence.level))))
}
plot(x=r,y=gpdf,xlim=c(0,maxl), type="l",
ylab="Density",
xlab="Likelihood Ratio Statistic",main=mainm,cex.main=1,sub="The vertical line is the observed likelihood ratio statistic")
graphics::abline(v=obsL[1],lty=2)
}
maxl <- stats::quantile(L[,2],0.95)
gpdf=bgk_kde(L[,2],n=2^(ceiling(maxl/log(2))),MIN=0,MAX=1.1*max(L[,2]))
maxl <- max(1.1*obsL[2],maxl)
r <- seq(0, maxl, length = binn + 1)[-1] - 0.5*maxl/binn
gpdf <- stats::spline(x=gpdf[1,],y=gpdf[2,],xout=r,method="natural")$y
gpdf[r>max(L[,2])] <- 0
gpdf[gpdf<0] <- 0
scalef <- binn/sum(gpdf)
gpdf <- gpdf * scalef
if(all(diff(x) <= 0)){
if(mean(obsL[2] <= L[,2]) < (1-confidence.level)){
mainm <- paste("Null Distribution of the Decreasing Test Statistic",
sprintf("\nWe reject the null hypothesis of equal proportions (at the %s%% level).",format(100*(1-confidence.level))))
}else{
mainm <- paste("Null Distribution of the Decreasing Test Statistic",
sprintf("\nWe do not reject the null hypothesis of equal proportions (at the %s%% level).",format(100*(1-confidence.level))))
}
plot(x=r,y=gpdf,xlim=c(0,maxl), type="l",
ylab="Density",
xlab="Likelihood Ratio Statistic",main=mainm,cex.main=1,sub="The vertical line is the observed likelihood ratio statistic")
graphics::abline(v=obsL[2],lty=2)
}
maxl <- stats::quantile(L[,4],0.95)
gpdf=bgk_kde(L[,4],n=2^(ceiling(maxl/log(2))),MIN=0,MAX=1.1*max(L[,4]))
maxl <- max(1.1*obsL[4],maxl)
r <- seq(0, maxl, length = binn + 1)[-1] - 0.5*maxl/binn
gpdf <- stats::spline(x=gpdf[1,],y=gpdf[2,],xout=r,method="natural")$y
gpdf[r>max(L[,4])] <- 0
gpdf[gpdf<0] <- 0
scalef <- binn/sum(gpdf)
gpdf <- gpdf * scalef
if(all(diff(x) >= 0) | all(diff(x) <= 0)){
mainm <- paste("Null Distribution of the Monotone Test Statistic",
sprintf("\nWe do not reject the null hypothesis of equal proportions (at the %s%% level).",format(100*(1-confidence.level))))
plot(x=r,y=gpdf,xlim=c(0,maxl), type="l",
ylab="Density",
xlab="Likelihood Ratio Statistic",main=mainm,cex.main=1,sub="The vertical line is the observed likelihood ratio statistic")
graphics::abline(v=obsL[4],lty=2)
}
}
invisible(list(pvalue.increasing=mean(obsL[1] <= L[,1]),
pvalue.decreasing=mean(obsL[2] <= L[,2]),
pvalue.overall=mean(obsL[3] <= L[,3]),
L=L,x=x,sigma=sigma))
}
#' Compute a test of trend in prevalences based on a likelihood-ratio statistic
#'
#' This function takes a series of point estimates and their associated standard errors and
#' computes the p-value for the test of a monotone decrease in the
#' population prevalences (in sequence order).
#' The p-value for a monotone increase is
#' also reported.
#' More formally, let the \eqn{K} population prevalences in sequence order be \eqn{p_1, \ldots, p_K}.
#' We test the null hypothesis:\cr
#' \deqn{H_0 : p_1 = \ldots = p_K}
#' vs
#' \deqn{H_1 : p_1 \ge p_2 \ldots \ge p_K}
#' with at least one equality strict. A likelihood ratio statistic for this test has
#' been derived (Bartholomew 1959).
#' The null distribution of the likelihood ratio statistic is very complex
#' but can be determined by a simple Monte Carlo process.\cr
#' We also test the null hypothesis:\cr
#' \deqn{H_0 : p_1 \ge p_2 \ldots \ge p_K}
#' vs
#' \deqn{H_1 : \overline{H_0}}
#' The null distribution of the likelihood ratio statistic is very complex
#' but can be determined by a simple Monte Carlo process.
#' The function requires the isotone library.
#'
#' @aliases LRT.value.trend
#' @param x A vector of prevalence estimates in the order (e.g., time).
#' @param sigma A vector of standard error estimates corresponding to \code{x}.
#' @return
#' A list with components
#' \itemize{ \item\code{pvalue.increasing}: The p-value for the test of a monotone increase in population prevalence.
#' \item\code{pvalue.decreasing}: The p-value for the test of a monotone decrease in population prevalence.
#' \item\code{L}: The value of the likelihood-ratio statistic.
#' \item\code{x}: The passed vector of prevalence estimates in the order (e.g., time).
#' \item\code{sigma} The passed vector of standard error estimates corresponding to \code{x}.
#' }
#'
#' @author Mark S. Handcock
#' @references Bartholomew, D. J. (1959). A test of homogeneity for ordered alternatives. Biometrika 46 36-48.
#' @keywords survey manip
#' @examples
#'
#' \dontrun{
#' x <- c(0.16,0.15,0.3)
#' sigma <- c(0.04,0.04,0.1)
#' LRT.value.trend(x,sigma)
#' }
#' @export
LRT.value.trend <- function(x,sigma){
a <- 1/sigma^2
# xbar <- mean(x)
xbar <- sum(a*x)/sum(a)
s <- seq_along(x)
Atot <- cbind(s[-length(s)], s[-1])
fit.ls1 <- isotone::activeSet(Atot, "LS", y = x, weights = a)
# m=fit.ls1$x
# LRT.increasing <- sum(a*(x-xbar)^2)-sum(a*(x-m)^2)
LRT.increasing <- sum(a*(x-xbar)^2)-fit.ls1$fval
Atot <- cbind(s[-1],s[-length(s)])
fit.ls2 <- isotone::activeSet(Atot, "LS", y = x, weights = a)
#m=fit.ls2$x
#LRT.decreasing <- sum(a*(x-xbar)^2)-sum(a*(x-m)^2)
LRT.decreasing <- sum(a*(x-xbar)^2)-fit.ls2$fval
LRT.overall <- sum(a*(x-xbar)^2)
LRT.monotone <- sum(a*(x-xbar)^2)-min(fit.ls1$fval,fit.ls2$fval,na.rm=TRUE)
H0.monotone <- min(fit.ls1$fval,fit.ls2$fval,na.rm=TRUE)
H0.increasing <- fit.ls1$fval
H0.decreasing <- fit.ls2$fval
return(c(LRT.increasing,LRT.decreasing,LRT.overall,LRT.monotone,H0.increasing,H0.decreasing,H0.monotone))
}
#
LRT.trend.null <- function(x,sigma,number.of.bootstrap.samples=5000,confidence.level=0.95,plot=NULL,seed=1,smooth=0.35){
# require(isotone, quietly=TRUE, warn.conflicts=FALSE)
# if(plot){
# suppressMessages(require(locfit, quietly = TRUE))
# }
a <- 1/sigma^2
xbar <- sum(a*x)/sum(a)
k <- length(x)
set.seed(seed)
s <- seq_along(x)
# Find best increasing means
Atot <- cbind(s[-length(s)], s[-1])
fit.ls1 <- isotone::activeSet(Atot, "LS", y = x, weights = a)
m.inc=fit.ls1$x
r <- matrix(stats::rnorm(k*number.of.bootstrap.samples),ncol=k)
r <- sweep(r,2,sigma,"*");r <- sweep(r,2,m.inc,"+")
L.inc <- t(apply(r,1,LRT.value.trend,sigma=sigma))
# Find best increasing means
Atot <- cbind(s[-1],s[-length(s)])
fit.ls2 <- isotone::activeSet(Atot, "LS", y = x, weights = a)
m.dec=fit.ls2$x
r <- matrix(stats::rnorm(k*number.of.bootstrap.samples),ncol=k)
r <- sweep(r,2,sigma,"*");r <- sweep(r,2,m.dec,"+")
L.dec <- t(apply(r,1,LRT.value.trend,sigma=sigma))
L <- L.inc
L[,6] <- L.dec[,6]
L[,7] <- pmin(L.inc[,7],L.dec[,7])
obsL <- LRT.value.trend(x,sigma)
#
options("scipen"=100)
cat(sprintf("In a test of the null hypothesis of an increasing trend in proportions against the complementary hypothesis, the p-value is %s.\n",format.pval(mean(obsL[5] <= L[,5]),eps=0.0001)))
if(mean(obsL[5] <= L[,5]) < (1-confidence.level)){
cat(sprintf("We reject the null hypothesis of an increasing trend (at the %s%% level).\n\n",format(100*(1-confidence.level))))
}else{
cat(sprintf("We do not reject the null hypothesis of an increasing trend (at the %s%% level).\n\n",format(100*(1-confidence.level))))
}
cat(sprintf("In a test of the null hypothesis of an decreasing trend in proportions against the complementary hypothesis, the p-value is %s.\n",format.pval(mean(obsL[6] <= L[,6]),eps=0.0001)))
if(mean(obsL[6] <= L[,6]) < (1-confidence.level)){
cat(sprintf("We reject the null hypothesis of an decreasing trend (at the %s%% level).\n\n",format(100*(1-confidence.level))))
}else{
cat(sprintf("We do not reject the null hypothesis of an decreasing trend (at the %s%% level).\n\n",format(100*(1-confidence.level))))
}
cat(sprintf("In a test of the null hypothesis of a monotone trend in proportions against the complementary hypothesis, the p-value is %s.\n",format.pval(mean(obsL[7] <= L[,7]),eps=0.0001)))
if(mean(obsL[7] <= L[,7]) < (1-confidence.level)){
cat(sprintf("We reject the null hypothesis of a monotone trend (at the %s%% level).\n\n",format(100*(1-confidence.level))))
}else{
cat(sprintf("We do not reject the null hypothesis of a monotone trend (at the %s%% level).\n\n",format(100*(1-confidence.level))))
}
#
#
#
if("estimates" %in% plot){
yminus <- x - 1.96*sigma
yplus <- x + 1.96*sigma
errorbar(x=seq_along(x),y=x,yminus=yminus,
yplus=yplus,xlab="sequence",ylab="estimate",
main="Trend of estimates")
}
if("distributions" %in% plot){
binn <- 500
if(var(L[,7]) > 0 ){
maxl <- stats::quantile(L[,7],0.95)
gpdf=bgk_kde(L[,7],n=2^(ceiling(maxl/log(2))),MIN=0,MAX=1.1*max(L[,7]))
maxl <- max(1.1*obsL[7],maxl)
r <- seq(0, maxl, length = binn + 1)[-1] - 0.5*maxl/binn
gpdf <- stats::spline(x=gpdf[1,],y=gpdf[2,],xout=r,method="natural")$y
gpdf[r>max(L[,7])] <- 0
gpdf[gpdf<0] <- 0
scalef <- binn/sum(gpdf)
gpdf <- gpdf * scalef
if(mean(obsL[7] <= L[,7]) < (1-confidence.level)){
mainm <- paste("Monotone Null Distribution of the Overall Test Statistic",
sprintf("\nWe reject the null hypothesis of a monotone trend (at the %s%% level).",format(100*(1-confidence.level))))
}else{
mainm <- paste("Monotone Null Distribution of the Overall Test Statistic",
sprintf("\nWe do not reject the null hypothesis of a monotone trend (at the %s%% level).",format(100*(1-confidence.level))))
}
plot(x=r,y=gpdf,xlim=c(0,maxl), type="l",
ylab="Density",
xlab="Likelihood Ratio Statistic",main=mainm,cex.main=1,sub="The vertical line is the observed likelihood ratio statistic")
graphics::abline(v=obsL[7],lty=2)
}
if(var(L[,5]) > 0){
maxl <- stats::quantile(L[,5],0.95)
gpdf=bgk_kde(L[,5],n=2^(ceiling(maxl/log(2))),MIN=0,MAX=1.1*max(L[,5]))
maxl <- max(1.1*obsL[5],maxl)
r <- seq(0, maxl, length = binn + 1)[-1] - 0.5*maxl/binn
gpdf <- stats::spline(x=gpdf[1,],y=gpdf[2,],xout=r,method="natural")$y
gpdf[r>max(L[,5])] <- 0
gpdf[gpdf<0] <- 0
scalef <- binn/sum(gpdf)
gpdf <- gpdf * scalef
if(mean(obsL[5] <= L[,5]) < (1-confidence.level)){
mainm <- paste("Monotone Null Distribution of the Increasing Test Statistic",
sprintf("\nWe reject the null hypothesis of an increasing trend (at the %s%% level).",format(100*(1-confidence.level))))
}else{
mainm <- paste("Monotone Null Distribution of the Increasing Test Statistic",
sprintf("\nWe do not reject the null hypothesis of an increasing trend (at the %s%% level).",format(100*(1-confidence.level))))
}
plot(x=r,y=gpdf,xlim=c(0,maxl), type="l",
ylab="Density",
xlab="Likelihood Ratio Statistic",main=mainm,cex.main=1,sub="The vertical line is the observed likelihood ratio statistic")
graphics::abline(v=obsL[5],lty=2)
}
if(var(L[,6]) > 0){
maxl <- stats::quantile(L[,6],0.95)
gpdf=bgk_kde(L[,6],n=2^(ceiling(maxl/log(2))),MIN=0,MAX=1.1*max(L[,6]))
maxl <- max(1.1*obsL[6],maxl)
r <- seq(0, maxl, length = binn + 1)[-1] - 0.5*maxl/binn
gpdf <- stats::spline(x=gpdf[1,],y=gpdf[2,],xout=r,method="natural")$y
gpdf[r>max(L[,6])] <- 0
gpdf[gpdf<0] <- 0
scalef <- binn/sum(gpdf)
gpdf <- gpdf * scalef
if(mean(obsL[6] <= L[,6]) < (1-confidence.level)){
mainm <- paste("Monotone Null Distribution of the Decreasing Test Statistic",
sprintf("\nWe reject the null hypothesis of an decreasing trend (at the %s%% level).",format(100*(1-confidence.level))))
}else{
mainm <- paste("Monotone Null Distribution of the Decreasing Test Statistic",
sprintf("\nWe do not reject the null hypothesis of an decreasing trend (at the %s%% level).",format(100*(1-confidence.level))))
}
plot(x=r,y=gpdf,xlim=c(0,maxl), type="l",
ylab="Density",
xlab="Likelihood Ratio Statistic",main=mainm,cex.main=1,sub="The vertical line is the observed likelihood ratio statistic")
graphics::abline(v=obsL[6],lty=2)
}
}
invisible(list(pvalue.increasing=mean(obsL[5] <= L[,5]),
pvalue.decreasing=mean(obsL[6] <= L[,6]),
pvalue.overall=mean(obsL[7] <= L[,7]),
L=L,x=x,sigma=sigma))
}
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.