R/samediff.R

Defines functions Tau Delta I11 I11.tau I12 I22 llSameDiff llds0 llDeltaInf tauRootDeltaInf ll2SameDiff print.samediff vcov.samediff coef.samediff logLik.samediff profile.samediff plot.profile.samediff summary.samediff print.summary.samediff pval.prof pval.vcov SE.samediff CI.vcov deltaRoot tauRoot pll.delta pll.tau confint.samediff contour.samediff ell.delta ell.tau sddvn plot.samediff samediffSim samediffPwr

Documented in confint.samediff plot.profile.samediff plot.samediff profile.samediff samediffPwr samediffSim summary.samediff

samediff <-
  function (nsamesame, ndiffsame, nsamediff, ndiffdiff, VCOV = TRUE)
{
  call <- m <- match.call()
  m[[1]] <- as.name("list")
  m <- eval.parent(m) # evaluate the *list* of arguments
  ss <- m[[1]]; ds <- m[[2]]; sd <- m[[3]]; dd <- m[[4]]
  data <- unlist(m)
  for(i in data) {
    if(i != trunc(i) | i < 0)
      stop("data have to be non-negative integers")
  }
  if(length(data[data == 0]) > 2)
    stop("Not enough information in data,",
         " when more than two entries are zero")
  
  vcov <- array(NA, dim = c(2,2))
  name <- c("tau", "delta")
  dimnames(vcov) <- list(name, name)
  se <- rep(NA, 2)
  conv <- NULL
  
  ## If the fraction of same answers to same-samples is larger than the
  ## fraction of same-answers to diff-samples, simple estimation of
  ## delta may proceed: 
  isdelta <- ifelse(ss/(ds + ss) > sd/(sd + dd), TRUE, FALSE)

  ## Estimate tau, delta and their covariance matrix:
  if(isdelta && all(data != 0)) {
    tau <- Tau(ss, ds)
    psd <- sd/(sd + dd)
    delta <- uniroot(f = Delta, lower = 0, upper = 10, tau = tau,
                      psd = psd, tol = 1e-6)$root
    i11 <- I11(tau, delta, ss, ds, sd, dd)
    i12 <- I12(tau, delta, sd, dd)
    i22 <- I22(tau, delta, sd, dd)
    vcov <- solve(cbind(c(i11, i12), c(i12, i22)))
    dimnames(vcov) <- list(name, name)
    se <- sqrt(diag(vcov))
    case <- 0
    logLik <- llSameDiff(tau, delta, ss, ds, sd, dd)
  }
  else if(ss == 0 && sd == 0) {
###    warning("No information on delta") ## parse this warning?
    tau <- 0
    delta <- NA
    logLik <- 0
    case <- 0.1
  }
  else if(ds == 0 && dd == 0) {
###    warning("No information on delta") ## parse this warning?
    tau <- Inf
    delta <- NA
    case <- 1
    logLik <- 0
  }
  else if(ds == 0 && sd == 0) {
    tau <- Inf
    delta <- Inf
    case <- 1.2
    logLik <- 0
  }
  else if(ss == 0 && ds == 0){
    tau <- Inf
    delta <- NA
    case <- 1.12
    tmp <- optim(c(5, 5), fn = llds0, sd = sd, dd = dd,
                 method = "L-BFGS-B", hessian = FALSE,
                 lower = c(1e-4, 1e-4), upper = c(100, 100), 
                 control = list(fnscale = -1))
    logLik <- tmp$value
    conv <- tmp$convergence
  }
  else if(sd == 0 && dd == 0) {
    delta <- NA
    tmp <- optim(1, fn = llDeltaInf, ss = ss, ds = ds,
                 method = "BFGS", hessian = VCOV,
                 control = list(fnscale = -1))
    tau <- tmp$par
    if(VCOV) {
      vcov[1, 1] <- -(tmp$hessian)^(-1)
      se <- c(sqrt(vcov[1, 1]), NA)
    }
    conv <- tmp$convergence
    case <- 1.3
    logLik <- llDeltaInf(tau, ss, ds)
  }
  else if(ds == 0) {
    tau <- Inf
    delta <- Inf
    case <- 1.22
    tmp <- optim(c(5, 5), fn = llds0, sd = sd, dd = dd,
                 method = "L-BFGS-B", hessian = FALSE,
                 lower = c(1e-4, 1e-4), upper = c(100, 100), 
                 control = list(fnscale = -1))
    logLik <- tmp$value
    conv <- tmp$convergence
  }
  else if(sd == 0) {
    delta <- Inf
    tmp <- optim(1, fn = llDeltaInf, ss = ss, ds = ds,
                 method = "BFGS", hessian = VCOV,
                 control = list(fnscale = -1))
    tau <- tmp$par
    if(VCOV) {
      vcov[1, 1] <- -(tmp$hessian)^(-1)
      se <- c(sqrt(vcov[1, 1]), NA)
    }
    conv <- tmp$convergence
    case <- 3
    logLik <- llDeltaInf(tau, ss, ds)
  }
  else if(dd == 0 || ss == 0 || !isdelta) {
    delta <- 0
    tmp <- optim(1, fn = llSameDiff, delta = 1e-4, ss = ss, ds = ds,
                 sd = sd, dd = dd, method = "BFGS",
                 hessian = VCOV,
                 control = list(fnscale = -1))
    tau <- tmp$par
    if(VCOV) {
      vcov[1, 1] <- -(tmp$hessian)^(-1)
      se <- c(sqrt(vcov[1, 1]), NA)
    }
    conv <- tmp$convergence
    case <- 2
    logLik <- llSameDiff(tau, 1e-4, ss, ds, sd, dd)
  }
  else
    stop("This should never occur! Please contact package maintainer!") 
  
  coef <- c(tau, delta)
  names(se) <- names(coef) <- name  
  fit <- list(coef = coef, se = se, vcov = vcov, data = data,
              test = "same-different",
              call = call, convergence = conv, case = case,
              logLik = logLik)
  class(fit) <- "samediff"
  fit
}

Tau <- function(ss, ds) { # compute tau
  frac <- (2 * ss + ds)/(2 * (ds + ss))
  tau <- sqrt(2) * qnorm(frac)
  tau
}

Delta <- function(d, tau, psd)  # compute d-prime
  pnorm((tau - d)/sqrt(2)) - pnorm(-(tau + d)/sqrt(2)) - psd


I11 <- function(tau, d, ss, ds, sd, dd) { ## 11-element of Fisher I-matrix
  E1 <- (dd + sd)^3/(2 * dd * sd)
  E2 <- (dnorm((tau - d)/sqrt(2)) + dnorm((-tau - d)/sqrt(2)))^2
  E3 <- 2*(ds + ss)^3/(ds * ss)
  E4 <- (dnorm(tau/sqrt(2)))^2
  E1 * E2 + E3 * E4
}

I11.tau <- function(tau, ss, ds) { ## Fisher information for tau
  E3 <- 2 * (ds + ss)^3/(ds * ss)
  E4 <- (dnorm(tau/sqrt(2)))^2
  E3 * E4
}

I12 <- function(tau, d, sd, dd) { ## 12 = 21-element of Fisher I-matrix
  E1 <- -(dd + sd)^3/(2 * dd * sd)
  E2 <- (dnorm((tau - d)/sqrt(2)))^2 - (dnorm((-tau - d)/sqrt(2)))^2
  E1 * E2
}

I22 <- function(tau, d, sd, dd) { ## 22-element of Fisher I-matrix
  E1 <- (dd + sd)^3/(2 * dd * sd)
  E2 <- (dnorm((tau - d)/sqrt(2)) - dnorm((-tau - d)/sqrt(2)))^2
  E1 * E2
}

llSameDiff <- function(tau, delta, ss, ds, sd, dd, verbose = FALSE) {
  if(any(c(tau, delta) <= 0)) return(-Inf)
  d <- delta
  sqrt.2 <- sqrt(2)
  Pss <- 2 * pnorm(tau/sqrt.2) - 1
  Pds <- 1 - Pss
  Psd <- pnorm((tau - d)/sqrt.2) - pnorm((-tau - d)/sqrt.2)
  Pdd <- 1 - Psd
  da <- c(ss, ds, sd, dd)
  pDa <- c(Pss, Pds, Psd, Pdd)
  da[pDa == 0] <- 0
  pDa[pDa == 0] <- 1
  ll <- ss * log(Pss) + ds * log(Pds) + sd * log(Psd) +
    dd * log(Pdd)
  if(verbose)
    cat(ll, tau, delta, Pss, Pds, Psd, Pdd, "\n")
  ll
}

llds0 <- function(par, sd, dd, verbose = FALSE) {
  if(any(par <= 0)) return(-Inf)
  d <- par[1]
  tau <- par[2]
  sqrt.2 <- sqrt(2)
  Psd <- pnorm((tau - d)/sqrt.2) - pnorm((-tau - d)/sqrt.2)
  Pdd <- 1 - Psd
  ll <- sd * log(Psd) + dd * log(Pdd)
  if(verbose)
    cat(ll, tau, d, Psd, Pdd, "\n")
  ll
}

llDeltaInf <- function(tau, ss, ds, verbose = FALSE) {
  if(any(tau <= 0)) return(-Inf)
  sqrt.2 <- sqrt(2)
  Pss <- 2 * pnorm(tau/sqrt.2) - 1
  Pds <- 1 - Pss
  ll <- ss * log(Pss) + ds * log(Pds)
  if(verbose)
    cat(ll, tau, Pss, Pds, "\n")
  ll
}

tauRootDeltaInf <- function(tau, ss, ds, pll.max, lim) {
### Solve this function for tau to obtain a confidence limit, when
### hat.delta == Inf
  pll <- llDeltaInf(tau, ss, ds)
  exp(pll - pll.max) - lim
}

ll2SameDiff <- function(delta, tau, ss, ds, sd, dd)
  llSameDiff(tau, delta, ss, ds, sd, dd) ## reverse arguments

print.samediff <-
  function(x, digits = max(3, getOption("digits") - 3), ...)
{
  cat("\nCall:\n", deparse(x$call), "\n\n", sep = "")
  if (length(coef(x))) {
    cat("Coefficients:\n")
    print.default(format(coef(x), digits = digits), print.gap = 2, 
                  quote = FALSE)
  }
  else cat("No coefficients\n")
  cat("\n")
  invisible(x)
}

vcov.samediff <- function(object, ...) {
  object$vcov
}

coef.samediff <- function(object, ...) {
  object$coef
}

logLik.samediff <- function(object, ...) {
  co <- coef(object)
  nobs <- sum(object$data)
  val <- object$logLik
  names(val) <- NULL
  attr(val, "nobs") <- nobs
  attr(val, "df") <- ifelse(is.na(co[2]), 1, 2)
### df = 2 even if parameters are known to be either 0 or Inf?
  class(val) <- "logLik"
  val
}

profile.samediff <-
  function(fitted, which = 1:2, max = 2, numpts = 100,
           max.delta = 10, max.tau = 10, ...)
{
  co <- fitted$coef
  da <- fitted$data
  
  ## Profiling tau:
  if(fitted$case == 1.12) {
    M <- max.tau
    ttau <- seq(1e-4, M, length = numpts) ## no evaluation at zero
    ll <- double(numpts)
    for(i in 1:numpts)
      ## optimize likelihood for delta at each of the tau-values.
      ll[i] <- pll.tau(delta = c(1e-2, M+2), tau = ttau[i], ss = da[1],
                       ds = da[2], sd = da[3], dd = da[4])
  }
  else if(is.na(co[2]) || co[1] == Inf && co[2] == Inf) {
    M <- max.tau
    ttau <- seq(1e-4, M, length = numpts) ## no evaluation at zero
    if(da[2] == 0 && da[3] == 0 || da[2]==0 && da[4]==0) 
      ll <- llDeltaInf(ttau, ss = da[1], ds = da[2])
    else {
      M <- max.tau
      ttau <- seq(1e-4, M, length = numpts) ## no evaluation at zero
      ll <- double(numpts)
      for(i in 1:numpts)
        ## optimize likelihood for delta at each of the tau-values.
        ll[i] <- pll.tau(delta = c(1e-2, M+2), tau = ttau[i], ss = da[1],
                         ds = da[2], sd = da[3], dd = da[4])
    }
  }
  else {
    if(co[1] == Inf)
      M <- max.tau
    else
      M <- co[1] + max
    ttau <- seq(1e-4, M, length = numpts) ## no evaluation at zero
    ll <- double(numpts)
    for(i in 1:numpts)
      ## optimize likelihood for delta at each of the tau-values.
      ll[i] <- pll.tau(delta = c(1e-2, M), tau = ttau[i], ss = da[1],
                       ds = da[2], sd = da[3], dd = da[4])
  }
  Data <- data.frame("tau" = ttau, "plTau" = ll)
  
  ## Profiling delta:
  if(!is.na(co[2])) {
    if(co[2] == Inf)
      M <- max.delta
    else
      M <- co[2] + max
    
    delta <- seq(1e-4, M, length = numpts)
    for(i in 1:numpts) {
      ll[i] <- pll.delta(tau = c(1e-4, M), delta = delta[i],
                         ss = da[1], ds =  da[2], sd = da[3], dd = da[4])
    }
    Data$delta <- delta
    Data$plDelta <- ll
  }
  attr(Data, "logLik") <- fitted$logLik
  class(Data) <- c("profile.samediff", "data.frame")
  ## Return:
  Data
}
  
plot.profile.samediff <-
  function(x, which = 1:nc, level = c(0.99, 0.95),
           fig = TRUE, ...)
{
  max.ll <- attr(x, "logLik")
  nc <- ncol(x)/2
  if(!all(which %in% 1:2) && length(which) <= 2)
    stop("Only 1 and 2 allowed in the vector 'which'",
         " with maximum two elements") 
  oldpar <- par(mfrow = c(1, length(which)))
  on.exit(par(oldpar))
  lim <- sapply(level, function(x) exp(-qchisq(x, df=1)/2) )
  x$nplTau <- exp(x[,2] - max.ll)
  x$nplDelta <- exp(x[,4] - max.ll)
  
  if(fig == TRUE) {
    if(any(which == 1)) {
      plot(x[,1], x$nplTau, type = "l", las = 1, ylim = c(0, 1), 
           xlab = expression(tau), 
           ylab = "Normalized Profile Likelihood",
         main = "")
      abline(h = lim, col = "grey")
    }
    if(any(which == 2)) {
      plot(x[,3], x$nplDelta, type = "l", las = 1, ylim = c(0, 1), 
           xlab = expression(delta),
           ylab = "Normalized Profile Likelihood",
           main = "")
      abline(h = lim, col = "grey")
    }
  }
  class(x) <- c("nProfile.samediff", "data.frame")
  invisible(x)
}

summary.samediff <- function(object, profile = TRUE, ...) {
  co <- coef(object)
  if(profile){
    CI <- confint(object, ...)
    pval <- pval.prof(object, ...)
  }
  else {
    CI <- CI.vcov(object, ...)
    pval <- pval.vcov(object, ...)
  }
  object$table <-
    cbind("Estimate" = co, "Std. Error" = SE.samediff(object), CI,
          "P-value" = pval) ## Names?
  object$AIC <- AIC(object)
  class(object) <- "summary.samediff"
  object
}

print.summary.samediff <-
  function(x, ..., digits = max(3, getOption("digits") - 3),
           signif.stars = getOption("show.signif.stars"))
{
  cat("\nCall:\n")
  cat(paste(deparse(x$call), sep = "\n", collapse = "\n"), 
      "\n\n", sep = "")
  cat("Coefficients\n")
  printCoefmat(x$table, digits = digits,
               signif.stars=signif.stars, tst.ind=integer(0),
               cs.ind=1:4,
               P.values=TRUE, has.Pvalue=TRUE,)
  cat("\nLog Likelihood:", round(x$logLik, digits),
      "\tAIC:", round(x$AIC, digits), "\n") 
  invisible(x)
}

pval.prof <- function(object, ...) {
  co <- coef(object)
  da <- object$data
  mll <- c(logLik(object))
  
  if(co[1] != 0) {
    if(co[1] == Inf && object$case %in% c(1, 1.2)) 
      PllTau <- llDeltaInf(1e-4, ss=da[1], ds=da[2])
    else 
      PllTau <- pll.tau(, 1e-4, ss=da[1], ds=da[2], sd=da[3],
                        dd=da[4])
    pval.tau <- 1 - pnorm(sqrt(2 * (mll - PllTau)))
  }
  else
    pval.tau <- 1

  if(is.na(co[2]))
    pval.delta <- NA
  else if(co[2] == 0)
    pval.delta <- 1
  else  {
    ## If delta > 0 and potentially Inf:
    PllDelta <-
      pll.delta(, 1e-4, ss=da[1], ds=da[2], sd=da[3], dd=da[4])
    pval.delta <- 1 - pnorm(sqrt(2 * (mll - PllDelta)))
  }
    
  pval <- c(pval.tau, pval.delta)
  pval
}

pval.vcov <- function(object, ...) {
  co <- coef(object)
  SE <- SE.samediff(object)
  wald.statistic <- co/SE
  pval <- 1 - pnorm(wald.statistic)
  pval
}

SE.samediff <- function(object, ...) {
  SE <- rep(NA, 2)
  var <- vcov(object)
  if(!is.na(var[1, 1]))
    SE[1] <- sqrt(var[1, 1])
  if(!is.na(var[2, 2]))
    SE[2] <- sqrt(var[2, 2])
  SE
}

CI.vcov <- function(object, level = .95, ...) {
  co <- coef(object)
  SE <- SE.samediff(object)
  lim <- qnorm((1 - level)/2, lower.tail = FALSE)
  CI <- matrix(rep(NA, 4), 2)
  dimnames(CI) <- list(c("", ""), c("Lower", "Upper"))

  if(!is.na(SE[1])) {
    CI[1, 1] <- co[1] - lim * SE[1]
    CI[1, 2] <- co[1] + lim * SE[1]
  }
  if(!is.na(SE[2])) {
    CI[2, 1] <- co[2] - lim * SE[2]
    CI[2, 2] <- co[2] + lim * SE[2]
  }

  CI
}

deltaRoot <- function(delta, ss, ds, sd, dd, pll.max, lim) {
### Solve this function for delta to obtain a confidence limit
  pll <- pll.delta( , delta, ss, ds, sd, dd)
  exp(pll - pll.max) - lim
}

tauRoot <- function(tau, ss, ds, sd, dd, pll.max, lim) {
### Solve this function for tau to obtain a confidence limit
  pll <- pll.tau( , tau, ss, ds, sd, dd)
  exp(pll - pll.max) - lim
}

pll.delta <- function(tau = c(1e-4, 10), delta, ss, ds, sd, dd) {
### Profile log-likelihood for delta
  optimize(llSameDiff, tau, delta = delta, ss = ss,
           ds = ds, sd = sd, dd = dd,
           maximum = TRUE)$objective
}

pll.tau <- function(delta = c(1e-4, 10), tau, ss, ds, sd, dd) {
### Profile log-likelihood for tau
  optimize(ll2SameDiff, delta, tau = tau, ss = ss,
           ds = ds, sd = sd, dd = dd, maximum = TRUE)$objective
}

confint.samediff <-
  function(object, parm = c("tau", "delta"), level = 0.95,
           max = c(10, 10), ...)
{
  co <- coef(object)
  da <- object$data
  ss <- da[1]; ds <- da[2]; sd <- da[3]; dd <- da[4]
  lim <- exp(-qchisq(level, df = 1)/2)
  lim2 <- exp(-qchisq(level, df = .5)/2)
  CI <- matrix(rep(c(0, Inf), each = 2), 2)
  lev.names <- paste(level, "% ", c("tau", "delta"), sep = "")
  dimnames(CI) <- list(lev.names, c("Lower", "Upper"))
  co[co == 0] <- 1e-4
  if(is.na(co[2])){ 
    parm <- ifelse("tau" %in% parm, "tau", "")
    CI[2, ] <- rep(NA, 2)
  }
  pll.max <- c(object$logLik)
  
  if("tau" %in% parm) {
    ## TAU:  
    if(object$case == 1.2) 
      CI[1, 1] <- uniroot(tauRootDeltaInf, c(1e-4, max[1]), ss, ds,
                          pll.max, lim)$root
    else if(object$case == 1.12) {
      co[1] <- max[1]
      if(tauRoot(tau = 1e-4, ss, ds, sd, dd, pll.max, lim) < 0)
        ## Search for root below MLE
        CI[1, 1] <- uniroot(tauRoot, c(1e-4, co[1]), ss, ds, sd, dd,
                            pll.max, lim)$root
    }
    else if(is.na(co[2]) && object$case != 1.3) {
      if(co[1] == Inf)
        CI[1, 1] <- uniroot(tauRootDeltaInf, c(1e-4, max[1]), ss, ds,
                            pll.max, lim)$root
      else
        CI[1, 2] <- uniroot(tauRootDeltaInf, c(1e-4, max[1]), ss, ds,
                            pll.max, lim)$root
    }
    else {
      if(co[1] != Inf)
        ## Search for upper limit
        CI[1, 2] <- uniroot(tauRoot, c(co[1], max[1]), ss, ds, sd, dd,
                            pll.max, lim)$root
      else co[1] <- max[1]
      if(tauRoot(tau = 1e-4, ss, ds, sd, dd, pll.max, lim) < 0)
        ## Search for root below MLE
        CI[1, 1] <- uniroot(tauRoot, c(1e-4, co[1]), ss, ds, sd, dd,
                            pll.max, lim)$root
    }
  }
  if("delta" %in% parm){
    ## DELTA:
    if(co[2] != Inf)
      ## Search for upper limit
      CI[2, 2] <- uniroot(deltaRoot, c(co[2], max[2]), ss, ds, sd, dd,
                          pll.max, lim)$root
    else co[2] <- max[2]
    if(deltaRoot(delta = 1e-4, ss, ds, sd, dd, pll.max, lim) < 0)
      ## Search for root below MLE
      CI[2, 1] <- uniroot(deltaRoot, c(1e-4, co[2]), ss, ds, sd, dd,
                          pll.max, lim)$root
  }
  ## Return confidence interval:
  CI
}

contour.samediff <-
  function(x, norm = TRUE, level = c(.90, .95, .99), 
           lim = list(tau = c(1e-4, 10), delta = c(1e-4, 10),
             length = c(100, 100)), ...)
{
  da <- x$data
  dots <- list(...)
  Tau <- seq(lim$tau[1], lim$tau[2], length = lim$length[1])
  Delta <- seq(lim$delta[1], lim$delta[2], length = lim$length[2])
  ll <- outer(Tau, Delta, llSameDiff, ss=da[1], ds=da[2], sd=da[3],
              dd=da[4]) ## log-likelihood
  if(norm) { ## Normalized likelihood
    ll.max <- max(c(x$logLik), max(ll[!is.na(ll)])) 
    ll <- exp(ll - ll.max)
    contour(Tau, Delta, ll,
            xlab = expression(tau),
            ylab = expression(delta),
            levels = 1 - level,
            labels = level, ...)
  }
  else {
    contour(Tau, Delta, ll,
          xlab = expression(tau),
          ylab = expression(delta), ...)
  }
  invisible(list(tau = Tau, delta = Delta, ll = ll))
}
  
ell.delta <- function(tau.hat, delta, ss, ds, sd, dd) {
### Estimated log-likelihood for delta
  llSameDiff(tau.hat, delta, ss, ds, sd, dd)
}

ell.tau <- function(tau.hat, delta, ss, ds, sd, dd){
### Estimated log-likelihood for tau
  llSameDiff(tau.hat, delta, ss, ds, sd, dd)
}

sddvn <- 
  function(sn, n1, dn, n2)
{
#### d' and variance of d' for same-different method.
#### 
  s <- sn/n1
  d <- dn/n2
  p1 <- s + (1 - s)/2
  ta <- qnorm(p1) * sqrt(2)
  del <- function(delta, ta, d) 
    pnorm((ta - delta)/sqrt(2)) -
      pnorm(( - ta - delta)/sqrt(2)) - d
  delta <- uniroot(del, interval = c(0, 10), ta = ta, d = d) 
  delta <- delta[[1]]
  fi <- dnorm(qnorm(p1))
  fi <- fi^2
  btau <- (s * (1 - s))/(2 * n1 * fi)
  bvp2 <- (d * (1 - d))/n2
  a <- (ta - delta)/sqrt(2)
  b <- ( - ta - delta)/sqrt(2)
  dx <- - dnorm(a)/sqrt(2) + dnorm(b)/sqrt(2)
  dy <- dnorm(a)/sqrt(2) + dnorm(b)/sqrt(2)
  bdelta <- bvp2/dx^2 + ((dy^2) * btau)/(dx^2)
  dv <- c(delta, bdelta)
##  dv <- round(dv, 3)
  dv
}

plot.samediff <-
  function(x, main = TRUE, length = 1000,
           limits, fig = TRUE, ...)
{
  co <- coef(x)
  is.delta <- !is.na(co[2]) && co[2] != 0 && co[2] != Inf
  if(missing(limits)) {
    limits <- c(-4, 4)
    if(is.delta)
      limits[2] <- co[2] + 4
  }
  z <- seq(limits[1], limits[2], length.out = length)
  base.dist <- dnorm(z)
  object <- data.frame(z, base.dist)
  if(fig == TRUE) {
    main.txt <- ifelse(main,
                       paste("Distribution of sensory intensity for the",
                             "same-different test"), c("") )
    plot(z, base.dist, type="l", xlab = "Sensory Magnitude",
###       ylab = "Sensory Intensity",
         ylab = "", main = main.txt, las = 1, lty = 2, ...)
    if(is.delta) {
      object$delta.dist <- dnorm(z, mean = coef(x)[2])
      lines(z, object$delta.dist, col = "red", lty = 1, ...)
    }
  }
  invisible(object)
}

samediffSim <- function(n, tau, delta, Ns, Nd) {
  m <- match.call()
  m[[1]] <- as.name("list")
  eval.parent(m)
  pss <- 2 * pnorm(tau/sqrt(2)) - 1
  ss <- rbinom(n = n, size = Ns, prob = pss)
  ds <-  Ns - ss
  psd <- pnorm((tau - delta)/sqrt(2)) -
    pnorm((-tau - delta)/sqrt(2))
  sd <- rbinom(n = n, size = Nd, prob = psd)
  dd <- Nd - sd
  cbind(ss, ds, sd, dd)
}


samediffPwr <- function(n = 1000, tau, delta, Ns, Nd, alpha = 0.05) {
  m <- match.call()
  m[[1]] <- as.name("list")
  eval.parent(m)
  sds <- samediffSim(n, tau, delta, Ns, Nd)
  pval <- double(n)
  fun <- function(x) {
    sadi <- samediff(x[1], x[2], x[3], x[4])
    co <- coef(sadi)
    mll <- c(logLik(sadi))
    if(is.na(co[2]))
      pval.delta <- NA
    else if(co[2] == 0)
      pval.delta <- 1
    else  {
      ## If delta > 0 and potentially Inf:
      PllDelta <-
        pll.delta(, 1e-4, ss=x[1], ds=x[2], sd=x[3], dd=x[4])
      pval.delta <- 1 - pnorm(sqrt(2 * (mll - PllDelta)))
    }
    pval.delta
  }
  pval <- apply(sds, 1, fun)
  power <- mean(pval < alpha, na.rm = TRUE)
  power
}

Try the sensR package in your browser

Any scripts or data that you put into this service are public.

sensR documentation built on Feb. 11, 2020, 1:10 a.m.