deprecated: Wiener likelihood and criterion functions (deprecated)
In yeagle/RWiener: Wiener Process Distribution Functions

Description Usage Arguments Details References Examples

wiener_likelihood computes the log-likelihood for given parameter values and data. wiener_deviance computes the deviance. wiener_aic computes the AIC. wiener_bic computes the BIC. These functions can be very useful in combination with the optim funcion, to estimate parameters (see example below).

  wiener_likelihood(x, data)
  wiener_deviance(x, data)
  wiener_aic(x, data, loss=NULL)
  wiener_bic(x, data, loss=NULL)

`x`	vector with the four parameter values: alpha, tau, beta, delta.
`data`	dataframe with data. Needs a reaction time column and a accuracy/response column.
`loss`	Defaults to NULL, which means that the default computation is done by using -2wiener_likelihood(x,data) in the formula. If not NULL, this can be a function to replace the default -2wiener_likelihood(x,data) in the code and use a custom function instead.

The described functions are deprecated, but still fully supported. They are kept for backwards compatibility and to ensure one can reproduce the examples from Wabersich & Vandekerckhove (2014).

These functions are simple wrapper functions for the generic R functions reported below.

The User is encouraged to use the generic R functions instead: logLik.wdm, deviance.wdm, AIC.wdm, BIC.wdm. See the corresponding help pages for more information on these functions.

Wabersich, D., & Vandekerckhove, J. (2014). The RWiener package: An R package providing distribution functions for the Wiener diffusion model. The R Journal, 6(1), 49-56.

### Example 1: Parameter estimation
## generate random data
dat <- rwiener(100,2,.3,.5,0)

## compute likelihood
wiener_likelihood(c(2,.3,.5,0), dat)

## estimate parameters with optim
onm <- optim(c(1,.1,.1,1),wiener_deviance,data=dat, method="Nelder-Mead")
est <- optim(onm$par,wiener_deviance,data=dat, method="BFGS",hessian=TRUE)
est$par # parameter estimates

## the following code needs the MASS package 
## Not run: sqrt(diag(MASS::ginv(est$hessian))) # sd for parameters

### Example 2: Simple model comparison
## compare two models with deviance
wiener_deviance(c(3,.3,.5,0), dat)
wiener_deviance(c(3,.3,.5,0.5), dat)
## log-likelihood difference
wiener_likelihood(c(3,.3,.5,0), dat)-wiener_likelihood(c(3,.3,.5,0.5), dat)

## Not run: 
### Example 3: likelihood-ratio test and Wald test
## Suppose we have data from 2 conditions
dat1 <- rwiener(100,2,.3,.5,-.5)
dat2 <- rwiener(100,2,.3,.5,.5)
onm1 <- optim(c(1,.1,.1,1),wiener_deviance,data=dat1, method="Nelder-Mead")
est1 <- optim(onm1$par,wiener_deviance,data=dat1, method="BFGS",hessian=TRUE)
wiener_likelihood(est1$par,dat1)+wiener_likelihood(est1$par,dat2) # combined loglike
model_ll <- function(pars,delta,dat1,dat2) {
  wiener_likelihood(pars,dat1)+
    wiener_likelihood(c(pars[1:3],pars[4]+delta),dat2)
}
## likelihood-ratio test
## 0-model: delta=0; alt-model: delta=1
model_ll(est1$par,1,dat1,dat2)
## compute likelihood ratio
LR <- -2*model_ll(est1$par,0,dat1,dat2)+2*model_ll(est1$par,1,dat1,dat2)
## compare with critical X^2(1) quantile, alpha=0.05
LR > qchisq(0.95,1)
## get p-value from X^2(1)
pchisq(LR,1, lower.tail=FALSE)
## Wald-Test
## estimate parameter delta and test for significance
onm2 <- optim(c(1,.1,.1,1),wiener_deviance,data=dat2, method="Nelder-Mead")
est2 <- optim(onm2$par,wiener_deviance,data=dat2, method="BFGS",hessian=TRUE)
delta <- est2$par[4]-est1$par[4]
## the following code needs the MASS package 
est1.sd <- sqrt(diag(MASS::ginv(est1$hessian))) # sd for parameters
WT <- (est1$par[4]-(est1$par[4]+delta))/est1.sd[4]
## compare with critical quantile N(0,1), alpha=0.05
abs(WT) > qnorm(0.975)
## get p-value from N(0,1)
pnorm(WT)

## End(Not run)

### Example 4: Custom AIC loss function
many_drifts <- function(x,datlist) {
  l = 0
  for (c in 1:length(datlist)) {
    l = l + wiener_deviance(x[c(1,2,3,c+3)],datlist[[c]])
  }
  return(l)
}
dat1 <- rwiener(n=100, alpha=2, tau=.3, beta=.5, delta=0.5)
dat2 <- rwiener(n=100, alpha=2, tau=.3, beta=.5, delta=1)
datlist <- list(dat1,dat2)
wiener_aic(x=c(2,.3,.5,.5,1), data=datlist, loss=many_drifts)