#################################
### Main estimation functions:
#################################
dod <-
function(same, diff, d.prime0 = 0, conf.level = 0.95,
statistic = c("likelihood", "Pearson", "Wilcoxon", "Wald"),
alternative = c("difference", "similarity", "two.sided",
"less", "greater"), control=dodControl(), ...)
### Our main fitting function.
### ... parsed on to wilcox.test when statistic = "Wilcoxon".
{
## Match arg:
alt <- match.arg(alternative)
if(alt == "difference") alt <- "greater"
if(alt == "similarity") alt <- "less"
stat <- match.arg(statistic)
ctrl <- control
if(!inherits(control, "dodControl"))
stop("Specify 'control' with dodControl()")
## Test arg:
if(ctrl$test.args)
testdodargs(same=NULL, diff=NULL, conf.level=conf.level,
d.prime0=d.prime0,
return.data=FALSE, int.tol=ctrl$integer.tol)
if(isTRUE(all.equal(d.prime0, 0)) &&
!alt %in% c("difference", "greater"))
stop("'alternative' has to be 'difference' or 'greater' if 'd.prime0' is 0")
if(!isTRUE(all.equal(d.prime0, 0)) && stat == "Wilcoxon") {
stop("Wilcoxon statistic only available with d.prime0 = 0")
}
## Fit DOD model:
fit <- dod_fit(same, diff, control=ctrl, ...)
## Format data:
.data <- fit$data
nlev <- ncol(.data)
npar <- length(fit$coefficients)
same <- as.vector(.data[1, ])
diff <- as.vector(.data[2, ])
logLik <- fit$logLik
d.prime <- fit$d.prime
## Compute standard errors:
if(d.prime < 0.01 && ctrl$get.vcov && ctrl$do.warn) {
warning("d.prime < 0.01: standard errors are unavailable")
}
std.err <- rep(NA_real_, nlev)
if(!is.null(fit$vcov) && all(is.finite(fit$vcov)))
std.err <- sqrt(diag(fit$vcov))
names(std.err) <- names(fit$coefficients)
## Compute confidence intervals:
if(stat == "Wald") {
ci <- confint(fit, level=conf.level, method="Wald")
conf.method <- "Wald"
} else {
ci <- confint(fit, level=conf.level, method="likelihood")
conf.method <- "profile likelihood"
}
### FIXME: add convergence information to model, so it can be
### accessed and assessed.
## Collect estimates in tables:
coef.tab <- cbind(d.prime, std.err[npar], ci[1L], ci[2L])
rownames(coef.tab) <- "d.prime"
colnames(coef.tab) <- c("Estimates", "Std. Error", "Lower", "Upper")
boundary.coeftab <- rbind(fit$tau, std.err[-npar])
rownames(boundary.coeftab) <- c("Estimate", "Std. Error")
colnames(boundary.coeftab) <-
as.character(seq_len(ncol(boundary.coeftab)))
## Results list:
fit$coefficients <- coef.tab
res <- c(fit,
list(conf.level=conf.level, boundary.coeftab=boundary.coeftab,
statistic=stat, d.prime0=d.prime0, alternative=alt,
conf.method=conf.method))
res$call <- match.call(expand.dots=TRUE)
## Get p-value and value of test statistic:
res$stat.value <- res$p.value <- NA_real_
if(stat == "Wilcoxon") {
## Assumes that all response categories have at least one
## positive record:
wtest <- wilcox.test(rep.int(1:nlev, diff), rep.int(1:nlev, same),
alternative=alt, ...)
res$stat.value <- unname(wtest$statistic)
res$p.value <- unname(wtest$p.value)
} else if(stat == "Wald") {
if(is.finite(std.err[npar])) {
res$stat.value <- unname((d.prime - d.prime0) / std.err[npar])
res$p.value <- unname(normalPvalue(res$stat.value, alternative=alt))
}
} else if(stat %in% c("likelihood", "Pearson")) {
## Fit model under the null hypothesis where d.prime = d.prime0:
ctrl$get.vcov <- ctrl$get.grad <- ctrl$test.args <- FALSE
fit0 <- dod_fit(same, diff, d.prime=d.prime0, control=ctrl, ...)
logLik0 <- fit0$logLik
d.prime0 <- fit0$d.prime
## Assess convergence of model under the alternative:
if(logLik < logLik0 && abs(logLik - logLik0) > 1e-6 && ctrl$do.warn)
warning("Estimation of DOD model failed: likelihood did not increase")
if(stat == "likelihood") {
LR <- 2 * (logLik - logLik0) ## LR shold be positive
LR[LR < 0 && abs(LR) < 1e-4] <- 0
if(LR >= 0) {
res$stat.value <- unname(sign(d.prime - d.prime0) * sqrt(LR))
res$p.value <- unname(normalPvalue(res$stat.value, alternative=alt))
}
} else if(stat == "Pearson") {
## Pearson equivalent of LR test:
freq <- par2prob_dod(tau=fit$tau, d.prime=d.prime) * rowSums(.data)
## Probabilities and frequencies under the null model (d.prime0):
freq0 <- par2prob_dod(tau=fit0$tau, d.prime=fit0$d.prime) *
rowSums(.data)
### FIXME: Check that this is the right test statistic (or should 'freq' be in the
### denominator?)
X2 <- sum((freq - freq0)^2 / freq0)
res$stat.value <- unname(sign(d.prime - d.prime0) * sqrt(X2))
res$p.value <- unname(normalPvalue(res$stat.value, alternative=alt))
}
} else stop("'statistic' not recognized")
## Return:
class(res) <- "dod"
res
}
dod_fit <-
function(same, diff, tau=NULL, d.prime=NULL, control=dodControl(),
...)
{
## Return list:
res <- list(call = match.call(expand.dots=TRUE))
res$control <- control
if(!inherits(control, "dodControl"))
stop("Specify 'control' with dodControl()")
## Read data:
x <- readdoddata(same, diff)
## Test arg:
if(control$test.args)
x <- testdodargs(x$same, x$diff, tau=tau, d.prime0=d.prime,
int.tol=res$control$integer.tol,
return.data=TRUE)
same <- x$same
diff <- x$diff
nlev <- length(same)
.data <- rbind(same, diff)
colnames(.data) <- as.character(seq_len(ncol(.data)))
rownames(.data) <- c("same-pairs", "diff-pairs")
res$data <- .data
res$tau <- as.vector(tau)
res$d.prime <- as.vector(d.prime)
### 3 CASES to get parameters:
if(is.null(tau) && !is.null(d.prime) &&
isTRUE(all.equal(d.prime, 0))) {
### CASE 1:
## d-prime = 0 and tau=NULL:
## No optimization necessary.
res$tau <- dod_null_tau_internal(same=same, diff=diff)
res$d.prime <- d.prime
} else if(is.null(tau) || is.null(d.prime)) {
### CASE 2:
## We need to optimize either tau, d.prime or c(tau, d.prime):
if(is.null(d.prime) && is.null(tau)) { ## get c(tau, d.prime)
start <- c(cumsum(init_tau(nlev)), 1)
par.names <-
c(paste("tau", seq_len(nlev-1L), sep=""), "d.prime")
npar <- length(start)
objfun <- function(par) {
dod_nll_internal(tau=par[1:(npar - 1L)], d.prime=par[npar],
same=same, diff=diff)
}
case <- "21"
} else if(is.null(tau) && !is.null(d.prime)) { ## get tau
start <- cumsum(init_tau(nlev))
par.names <- paste("tau", seq_len(nlev-1L), sep="")
objfun <- function(par)
dod_nll_internal(tau=par, d.prime=d.prime,
same=same, diff=diff)
case <- "22"
} else if(!is.null(tau) && is.null(d.prime)) { ## get d.prime
start <- 1
par.names <- "d.prime"
objfun <- function(par)
dod_nll_internal(tau=tau, d.prime=par,
same=same, diff=diff)
case <- "23"
}
res$start <- setNames(start, par.names)
npar <- length(start) ## Perhaps not needed
## Optimization:
res$optRes <- nlminb(res$start, objfun, lower=0, control=control$optCtrl)
par <- as.vector(res$optRes$par)
## Convergence tests:
start.value <- objfun(res$start)
if(all(par > 1e-2) && control$get.grad){
g <- grad(objfun, x=par)
res$gradient <- setNames(g, par.names)
if(!all(is.finite(g)) && control$do.warn) {
warning("Cannot assess convergence: non-finite gradient")
} else {
if(max(abs(g)) > control$grad.tol && control$do.warn)
warning("Estimation failed with max(gradient) = ",
format(max(abs(g)), digits=2), " (grad.tol = ",
format(control$grad.tol, digits=2), ")")
if(control$get.vcov) {
## Hessian and vcov:
h <- res$vcov <- array(NA_real_, dim=c(npar, npar),
dimnames=list(par.names, par.names))
h <- try(hessian(objfun, x=par), silent=TRUE)
if((inherits(h, "try-error") || !all(is.finite(h)))
&& control$do.warn) {
warning("unable to compute Hessian")
} else {
ch <- try(chol(h), silent=TRUE)
if(inherits(ch, "try-error") && control$do.warn) {
warning("Model is ill-defined and may not have converged")
} else {
res$vcov[] <- chol2inv(ch)
}
dimnames(h) <- list(par.names, par.names)
res$Hessian <- h
}
}
}
}
## Match parameters:
if(case == "21") {
res$tau <- par[1:(npar-1L)]
res$d.prime <- par[npar]
} else if(case == "22") {
res$tau <- par
res$d.prime <- as.vector(d.prime)
} else { ## case == "23"
res$tau <- as.vector(tau)
res$d.prime <- par
}
}
### CASE 3:
## Both tau and d.prime set, so we just evaluate the likelihood
## No optimization necessary
res$logLik <- -dod_nll_internal(tau=res$tau, d.prime=res$d.prime,
same=same, diff=diff)
## Set parameter names:
res$coefficients <-
setNames(c(res$tau, res$d.prime),
c(paste("tau", seq_along(res$tau), sep=""), "d.prime"))
## Return fit:
class(res) <- "dod_fit"
res
}
#################################
### dodControl:
#################################
dodControl <- function(grad.tol = 1e-4,
## eigen.value.tol=1e-4,
integer.tol = 1e-8,
get.vcov = TRUE,
get.grad = TRUE,
test.args = TRUE,
do.warn=TRUE,
optCtrl=list())
{
stopifnot(length(grad.tol) == 1L,
grad.tol > 0,
is.finite(grad.tol))
stopifnot(length(integer.tol) == 1L,
integer.tol > 0,
is.finite(integer.tol))
stopifnot(length(get.vcov) == 1L,
is.logical(get.vcov))
stopifnot(length(get.grad) == 1L,
is.logical(get.grad))
stopifnot(length(test.args) == 1L,
is.logical(test.args))
stopifnot(length(do.warn) == 1L,
is.logical(do.warn))
## Return value:
res <- list(grad.tol=grad.tol,
integer.tol=integer.tol,
get.grad=get.grad,
get.vcov=get.vcov,
test.args=test.args,
do.warn=do.warn,
optCtrl=optCtrl)
class(res) <- "dodControl"
res
}
#################################
### Optimal tau estimation:
#################################
optimal_tau <-
function(d.prime, d.prime0 = 0, ncat=3,
method=c("equi.prob", "LR.max", "se.min"),
tau.start=NULL, equi.tol = 1e-4, grad.tol = 1e-2,
do.warn=TRUE)
### FIXME: Could add "Pearson.max" as an option here.
### Estimate optimal boundary parameters, tau in a DOD model with one
### of the following optimality criteria:
### 1) Equal category probabilities averaged over same-pairs and diff-pairs.
### 2) Minimum standard error of d.prime/Maximum Wald test.
### 3) Maximum LR statistics for the test: H0: d.prime = 0, HA: d.prime > 0.
{
## Define objective functions:
nlr_stat <- function(tpar) {
## Negative LR statistic.
tau <- cumsum(tpar)
## Limiting distribution of data given c(tau, d.prime):
data <- par2prob_dod(tau=tau, d.prime=d.prime) * 100
## logLiks at d.prime=d.prime and d.prime=0:
logLik <- -dod_nll_internal(tau=tau, d.prime=d.prime, same=data[1, ], diff=data[2, ])
if(isTRUE(all.equal(d.prime0, 0)))
logLik0 <- -dod_null_internal(same=data[1, ], diff=data[2, ])
else
logLik0 <- dod_fit(same=data[1, ], diff=data[2, ],
d.prime=d.prime0, control=ctrl)$logLik
-(logLik - logLik0) ## negative LR/2 statistic
}
min_se <- function(tpar) {
## function to minimize se(d.prime):
if(any(!is.finite(tpar)))
return(Inf)
tau <- cumsum(tpar)
data <- par2prob_dod(tau, d.prime) * 100
h <- try(hessian(dod_nll_all_internal, c(tau, d.prime), same=data[1, ],
diff=data[2, ]), silent=TRUE)
if(inherits(h, "try-error") || any(!is.finite(h)))
return(Inf)
vcov <- try(solve(h), silent=TRUE)
if(inherits(vcov, "try-error") || any(!is.finite(vcov)))
return(Inf)
sqrt(rev(diag(solve(h)))[1])
}
equi_prob <- function(tpar, const=1) {
tau <- cumsum(tpar)
prob <- colSums(par2prob_dod(tau, d.prime)) / 2
const + sum(((prob - 1/ncat) * 1e3)^2)
}
tau2tpar <- function(tau) c(tau[1], diff(tau))
## tpar2tau <- function(tpar) cumsum(tpar)
## Match and test arguments:
method <- match.arg(method)
stopifnot(is.numeric(d.prime),
length(d.prime) == 1L,
d.prime >= 0) ## OBS: strictly larger than 0!
stopifnot(is.numeric(ncat),
length(ncat) == 1L,
ncat >= 2)
stopifnot(length(do.warn) == 1L,
is.logical(do.warn))
if(abs(ncat - round(ncat)) > 1e-8 && do.warn) {
warning("Non-integer ncat rounded to ", round(ncat))
}
ncat <- as.integer(round(ncat))
## Match objective function:
objfun <- switch(method, "equi.prob" = equi_prob,
"se.min" = min_se,
"LR.max" = nlr_stat)
if(method == "LR.max") {
ctrl <- dodControl(test.args=FALSE, do.warn=FALSE,
get.grad=FALSE, get.vcov=FALSE,
integer.tol=10,
grad.tol=grad.tol)
}
## Get starting values:
if(!is.null(tau.start)) { ## 1, 2, 3 -> c(ts[1], diff(ts))
stopifnot(is.numeric(tau.start),
all(is.finite(tau.start)),
all(tau.start > 0))
if(any(diff(tau.start) <= 0))
stop("'tau.start' is not increasing")
if(length(tau.start) != ncat-1L)
stop("length(tau.start) is ", length(tau.start),
": was expecting length(tau.start) = ncat - 1 = ",
ncat - 1L)
start <- tau2tpar(tau.start)
} else
start <- init_tau(ncat)
## Estimate optimal tau:
fit <- nlminb(start=start, objective=objfun, lower=1e-4)
## if(fit$convergence != 0)
## warning("Estimation of tau failed to converge with nlminb code ",
## fit$convergence)
tau <- cumsum(fit$par)
prob <- par2prob_dod(tau, d.prime)
if(method == "equi.prob") {
diffs <- abs((colSums(prob) / 2) - 1/ncat)
if(max(diffs) > equi.tol && do.warn)
warning("Estimation of tau failed with max(diffs) = ",
format(max(diffs), digits=2), " (equi.tol = ",
format(equi.tol, digits=2), ")")
}
## Obtain and check gradient:
g <- grad(func=objfun, x=fit$par)
if(max(abs(g)) > grad.tol && do.warn)
warning("Estimation of tau failed with max(gradient) = ",
format(max(abs(g)), digits=2), " (grad.tol = ",
format(grad.tol, digits=2), ")")
## Results:
res <- list(tau=tau, prob=prob, tau.start=cumsum(start), gradient=g,
fit=fit, method=method, call=match.call(expand.dots=TRUE),
grad.tol=grad.tol, equi.tol=equi.tol)
res
}
#################################
### Misc estimation functions:
#################################
dod_null_tau_internal <- function(same, diff) {
### Compute tau-vector for the DOD null model (d.prime=0)
alldata <- rev(same + diff)
nlev <- length(alldata)
cumdata <- cumsum(alldata)
tau <- sapply(1:(nlev - 1), function(lev) {
-sqrt(2) * qnorm(cumdata[lev] / (2* sum(alldata)))
})
tau <- rev(tau)
tau
}
dod_null_tau <- function(same, diff) {
### Compute tau-vector for the DOD null model (d.prime=0)
x <- readdoddata(same, diff)
x <- testdodargs(x$same, x$diff, int.tol=1) ## don't warn with non-integers
dod_null_tau_internal(x$same, x$diff)
}
par2prob_dod_internal <- function(tau, d.prime) {
## Cumulative probilities:
gamma.same <- 2*pnorm(tau / sqrt(2)) - 1
gamma.diff <- pnorm((tau - d.prime)/sqrt(2)) -
pnorm((-tau - d.prime)/sqrt(2))
## Non-cumulative probabilities:
p.same <- c(gamma.same, 1) - c(0, gamma.same)
p.diff <- c(gamma.diff, 1) - c(0, gamma.diff)
## Make sure sum(p.same) == sum(p.diff) == 1:
p.same <- p.same / sum(p.same)
p.diff <- p.diff / sum(p.diff)
## Return prob matrix:
rbind(p.same, p.diff)
}
par2prob_dod <- function(tau, d.prime) {
### Computes the dod probability matrix from dod parameters tau and
### d.prime; rbind(same-pair probs, diff-pair probs);
### Result satisfies: 1) rowSums(<return_value>) == c(1, 1),
### 2) ncol(<return_value>) == length(tau) + 1L.
stopifnot(length(d.prime) == 1L,
is.numeric(d.prime),
d.prime >= 0)
stopifnot(is.numeric(tau),
length(tau) >= 1L,
all(tau > 0),
all(diff(tau) > 0))
par2prob_dod_internal(tau, d.prime)
}
init_tau <- function(ncat=4) c(1, rep(3/(ncat-1), ncat-2))
#################################
### Likelihood functions:
#################################
dod_null_internal <- function(same, diff) {
### Compute neg. logLik for DOD null model (with d.prime=0)
tau <- dod_null_tau(same, diff)
dod_nll_internal(tau=tau, d.prime=0, same=same, diff=diff)
}
dod_null <- function(same, diff, integer.tol=1e-8) {
### Like dod_null_internal, but also tests admissibility of
### arguments.
### Computes the neg. logLik for the DOD null model (with d.prime=0).
### same and diff can be factors or numeric vectors.
## Read same & diff:
x <- readdoddata(same, diff)
## Test same and diff args:
x <- testdodargs(x$same, x$diff, int.tol=integer.tol)
## Compute logLik:
dod_null_internal(x$same, x$diff)
}
dod_nll_internal <- function(tau, d.prime, same, diff) {
### Computes the nll of the dod Thurstonian model.
### Returns Inf if parameters are inadmissible.
prob <- par2prob_dod_internal(tau, d.prime)
## compute nll:
if(all(is.finite(prob)) && all(prob > 0))
## if(all(prob > 0))
-sum(same %*% log(prob[1, ])) -
sum(diff %*% log(prob[2, ]))
else
Inf
}
dod_nll <- function(tau, d.prime, same, diff, integer.tol=1e-8) {
### Computes the neg. logLik of the dod Thurstonian model.
### Returns an error if the parameters are inadmissible.
## Get case probabilities:
prob <- par2prob_dod(tau, d.prime)
### NOTE: tau and d.prime args are tested in par2prob_dod().
if(any(prob <= 0)) stop("Inadmissible parameters")
## Read same & diff:
x <- readdoddata(same, diff)
## Test same and diff args:
x <- testdodargs(x$same, x$diff, tau=tau, int.tol=integer.tol)
## Compute negative log-likelihood:
nll <- -sum(x$same %*% log(prob[1, ])) -
sum(x$diff %*% log(prob[2, ]))
nll
}
dod_nll_all_internal<- function(par, same, diff) {
### nll expressed as a function of all=c(tau, d.prime) parameters.
npar <- length(par)
dod_nll_internal(tau=par[1:(npar-1)], d.prime=par[npar], same, diff)
}
#################################
### Utils:
#################################
readdoddata <- function(same, diff) {
### Read same & diff arguments to dod() etc.
### Returns a list of same and diff numeric vectors of observations.
if(is.factor(same) && is.factor(diff)) {
stopifnot(nlevels(same) == nlevels(diff))
stopifnot(all(levels(same) == levels(diff)))
same <- as.vector(table(same))
diff <- as.vector(table(diff))
}
else {
same <- as.vector(same)
diff <- as.vector(diff)
}
list(same=same, diff=diff)
}
testdodargs <-
function(same=NULL, diff=NULL, tau=NULL, conf.level=NULL,
d.prime0=NULL, int.tol=1e-8, return.data=TRUE)
{
### Tests admissibility of same and diff arguments.
stopifnot(length(return.data) == 1L,
is.logical(return.data))
## Test integer.tol argument:
stopifnot(is.numeric(int.tol),
length(int.tol) == 1L,
int.tol >= 0)
## Test same & diff:
has.data <- !is.null(same) && !is.null(diff)
if(has.data) {
stopifnot(is.numeric(same), is.numeric(diff),
all(same >= 0),
all(diff >= 0),
length(same) == length(diff))
if(!is.null(tau))
stopifnot(length(same) == length(tau) + 1L)
.data <- rbind(same, diff)
.data <- .data[, colSums(.data) > 0, drop=FALSE]
### FIXME: warn if any(colSums(.data) <= 0)?
if(ncol(.data) <= 1L)
stop("Need counts in more than one response category", call.=FALSE)
if(any(abs(.data - round(.data)) > int.tol))
warning("non-integer counts in 'same' or 'diff'", call.=FALSE)
}
## conf.level:
if(!is.null(conf.level)) {
stopifnot(is.numeric(conf.level),
length(conf.level) == 1,
0 < conf.level,
conf.level < 1)
}
## d.prime0:
if(!is.null(d.prime0)) {
stopifnot(is.numeric(d.prime0),
length(d.prime0) == 1,
d.prime0 >= 0)
}
## Return modified data:
if(has.data && return.data)
list(same=as.vector(drop(.data[1, ])),
diff=as.vector(drop(.data[2, ])))
else
invisible()
}
#################################
### Print methods:
#################################
print.dod <-
function(x, digits = max(3, getOption("digits") - 3), ...)
{
cat(paste("\nResults for the Thurstonian model for the Degree-of-Difference method",
"\n\n"))
cat(paste("Confidence level for 2-sided ", x$conf.method,
" interval: ", round(100 * x$conf.level, 3), "%\n\n",
sep=""))
print(x$coefficients, digits=digits)
cat("\nBoundary coefficients:\n")
print(x$boundary.coeftab, digits=digits)
cat("\nData:\n")
print(as.data.frame(x$data))
cat("\nResults of discrimination test:\n")
if(x$statistic == "Wald")
cat(paste("Wald statistic = ", format(x$stat.value, digits),
", p-value = ", format.pval(x$p.value, digits=4), "\n",
sep=""))
if(x$statistic == "likelihood")
cat(paste("Likelihood Root statistic = ",
format(x$stat.value, digits),
", p-value = ", format.pval(x$p.value, digits=4), "\n",
sep=""))
if(x$statistic == "Pearson")
cat(paste("Signed Pearson root statistic = ",
format(x$stat.value, digits),
", p-value = ", format.pval(x$p.value, digits=4), "\n",
sep=""))
if(x$statistic == "Wilcoxon")
cat(paste("Wilcoxon statistic = ", x$stat.value,
"p-value =", format.pval(x$p.value, digits=4), "\n"))
cat("Alternative hypothesis: ")
alt <- x$alternative
cat(paste("d-prime is",
ifelse(alt == "two.sided", "different from",
paste(alt, "than")), format(x$d.prime0, digits),
"\n"))
invisible(x)
}
#################################
### Simulation functions:
#################################
get_tau <-
function(d.prime, tau=NULL, ncat=NULL,
method.tau = c("LR.max", "equi.prob", "se.min", "user.defined"),
d.prime0=0, ...)
### Internal function.
### Computes tau if not already present and tests that arguments
### conform.
{
method <- match.arg(method.tau)
## Test conformance of tau and method.tau:
if(method != "user.defined") {
if(!is.null(tau))
warning("'tau' is ignored when method.tau != 'user.defined' (method.tau was '",
method, "')")
stopifnot(length(ncat) == 1L,
is.numeric(ncat),
ncat >= 2L,
is.finite(ncat))
if(abs(ncat - round(ncat)) > 1e-6)
warning("non-integer 'ncat': ", format(ncat, digits=8),
" rounded to ", round(ncat))
ncat <- round(ncat)
tau <- optimal_tau(d.prime=d.prime, ncat=ncat, method=method,
do.warn=FALSE, d.prime0=d.prime0, ...)$tau ## grad.tol, equi.tol
} else { ## method == "user.defined"
if(is.null(tau))
stop("specify 'tau' when method.tau == 'user.defined'")
stopifnot(length(tau) >= 1L,
is.numeric(tau),
all(diff(tau) > 0),
all(is.finite(tau)))
### FIXME: Do we need more tests of tau here?
}
## Return:
as.vector(tau)
}
dodSim <-
function(d.prime, ncat=4, sample.size = c(100, 100),
method.tau = c("equi.prob", "LR.max", "se.min", "user.defined"),## max.X2?
tau = NULL, d.prime0 = 0, ...)
{
## Test d.prime:
stopifnot(length(d.prime) == 1L,
is.numeric(d.prime),
d.prime >= 0,
is.finite(d.prime))
stopifnot(length(d.prime0) == 1L,
is.numeric(d.prime0),
d.prime0 >= 0,
is.finite(d.prime0))
## Match and test tau-related args:
method <- match.arg(method.tau)
tau <- get_tau(d.prime=d.prime, tau=tau, ncat=ncat,
method.tau=method, d.prime0=d.prime0)
## Test sample.size arg:
stopifnot(is.numeric(sample.size),
all(sample.size >= 0),
length(sample.size) == 1L || length(sample.size) == 2L)
### NOTE: sample.size == 0 is allowed and runs with expected results.
size <- sample.size
if(max(abs(size - round(size))) > 1e-3)
warning("Non-integer 'sample.size' rounded to ", round(size))
size <- round(size)
if(length(size) == 1L) size <- c(size, size)
## Get prob and simulate data:
prob <- par2prob_dod(tau, d.prime)
### NOTE: tau and d.prime args are tested in par2prob_dod().
mat <- t(cbind(rmultinom(n=1, size=size[1], prob=prob[1, ]),
rmultinom(n=1, size=size[2], prob=prob[2, ])))
dimnames(mat) <- list(c("same-pairs", "diff-pairs"),
seq_len(ncol(mat)))
## Return:
mat
}
#################################
### profile likelihood CI
#################################
lroot <- function(d.prime, same, diff) {
## Profile likelihood root statistic for d.prime
## Currently not used.
ctrl <- dodControl(get.grad=FALSE, get.vcov=FALSE, do.warn=FALSE,
test.args=FALSE)
fit <- dod_fit(same, diff, control=ctrl)
logLik <- fit$logLik
mle <- fit$d.prime
logLik0 <- dod_fit(same, diff, d.prime=d.prime, control=ctrl)$logLik
lroot <- sign(mle - d.prime) * sqrt(2 * abs(logLik - logLik0))
lroot
}
get_ci <-
function(same, diff, alpha = 0.05, dod.fit=NULL,
tol.d.prime = 1e-3)
{
### Computes the profile likelihood CI for d-prime.
lroot_fun <- function(d.prime, lim) {
## Profile likelihood root statistic minus lim.
ll <- dod_fit(same=same, diff=diff, d.prime=d.prime,
control=ctrl)$logLik
lroot <- sign(mle - d.prime) * sqrt(2 * abs(dod.fit$logLik - ll))
lroot - lim
}
lim <- qnorm(p=1-alpha/2, lower.tail=TRUE)
ctrl <- dodControl(get.vcov=FALSE, get.grad=FALSE, do.warn=FALSE,
test.args=FALSE)
## Get fit if missing:
if(is.null(dod.fit)) {
dod.fit <- dod_fit(same, diff, control=ctrl)
} else {
stopifnot(class(dod.fit) %in% c("dod", "dod_fit"))
}
mle <- dod.fit$d.prime
## Get lower and upper limits for the lroot search:
lwr <- lroot_lwr(same, diff, mle=mle, lim=lim,
logLik=dod.fit$logLik, tol.d.prime=tol.d.prime)
if(length(lwr) == 2)
lwr <- uniroot(f=lroot_fun, interval=lwr, tol=1e-4, lim=lim)$root
upr <- lroot_upr(same, diff, mle=mle, lim=lim,
logLik=dod.fit$logLik, tol.d.prime=tol.d.prime)
if(length(upr) == 2)
upr <- uniroot(f=lroot_fun, interval=upr, tol=1e-4, lim=-lim)$root
## Return CI:
ci <- as.vector(c(lwr, upr))
ci
}
lroot_upr <- function(same, diff, mle, lim, logLik, tol.d.prime=1e-3, ...)
{
if(is.na(mle)) return(NA_real_)
if(mle == Inf) return(Inf)
ctrl <- dodControl(get.vcov=FALSE, get.grad=FALSE, do.warn=FALSE,
test.args=FALSE)
Upr <- 10
if(mle > Upr) {
warning("Cannot determine upper confidence bound")
return(NA_real_)
} else {
ll.Upr <- dod_fit(same=same, diff=diff, d.prime=Upr,
control=ctrl)$logLik
LR.Upr <- 2 * (logLik - ll.Upr)
if(LR.Upr < lim^2) {
warning("Cannot determine upper confidence bound")
return(NA_real_)
}
}
c(mle, Upr)
}
lroot_lwr <- function(same, diff, mle, lim, logLik, tol.d.prime=1e-3, ...)
{
ctrl <- dodControl(test.args=FALSE, get.vcov=FALSE,
get.grad=FALSE, do.warn=FALSE)
if(is.na(mle)) return(NA_real_)
if(mle < tol.d.prime) return(0)
LR0 <- 2 * (logLik + dod_null_internal(same, diff))
if(LR0 < lim^2) return(0)
## positive lwr exists:
Upr <- mle
if(mle == Inf) {
## determine upper bound for search
LR10 <- 2 * (logLik +
dod_fit(same=same, diff=diff, d.prime=10,
control=ctrl)$logLik)
if(LR10 > lim^2) {
warning("Cannot determine lower confindence bound")
return(NA_real_)
} else {
Upr <- 10
}
}
c(0, Upr)
}
confint.dod_fit <-
function(object, parm, level=0.95,
method=c("likelihood", "Wald"), ...)
{
x <- object
method <- match.arg(method)
stopifnot(length(level) == 1L,
is.numeric(level),
level > 0,
level < 1)
## test level.
Data <- x$data
if(method == "Wald") {
if(is.null(x$vcov)) {
## Refit model to get vcov/std.err.:
ctrl <- x$control
ctrl$get.grad <- ctrl$get.vcov <- TRUE
ctrl$test.args <- FALSE
x <- dod_fit(same=Data[1, ], diff=Data[2, ], control=ctrl)
}
std.err <- NA_real_
if(!is.null(x$vcov) && all(is.finite(x$vcov)))
std.err <- sqrt(rev(diag(x$vcov))[1L])
ci <- rep(NA_real_, 2L)
if(is.finite(std.err)) {
a <- 1 - level
ci <- x$d.prime + qnorm(1 - a/2) * c(-1, 1) * std.err
ci[1L] <- delimit(ci[1L], lower=0)
}
} else if(method == "likelihood") {
ci <- get_ci(same=Data[1, ], diff=Data[2, ], alpha=1 - level,
dod.fit=x)
} else
stop("'method' not recognized")
setNames(ci, nm=c("lower", "upper"))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.