# R/qdm.R In qdm: Fitting a Quadrilateral Dissimilarity Model to Same-Different Judgments

#### Documented in persp.qdmpredict.qdmqdm

```# qdm.R
#
# last mod: Aug/21/2014, NU

## QDM function
qdmfun <- function(x, y, p, response = c("logistic", "guessing", "gumbel",
"gompertz", "weibull", "cauchy", "shepardA", "shepardAneg",
"shepardB", "shepardBneg", "shepardD", "shepardDneg",
"shepardE", "shepardEneg", "shepardF", "shepardFneg"),
bias = 0) {

s <- p[1] + p[2]*x + p[3]*x^2 + 2*p[4]*abs((x-bias) - y) + p[5]*y + p[6]*y^2

response <- match.arg(response)
switch(EXPR = response,
logistic = 1/(1 + exp(-p[7]*s - p[8])),
guessing = p[9] + (1 - p[9])*1/(1 + exp(-p[7]*s - p[8])),
gumbel = exp(-exp((-p[7]*s - p[8])/p[9])),
gompertz = p[7]*exp(-p[8]*exp(-p[9]*s)),
weibull = pweibull(s, shape=p[7], scale=p[8]),
cauchy = pcauchy(s, location=p[7], scale=p[8]),
shepardA = 1 - s*p[7] + (s*p[7])*log(abs(s*p[7])),
shepardAneg = 1 - (1 - s*p[7] + (s*p[7])*log(abs(s*p[7]))),
shepardB = 1 - (s*p[7])^2 + p[8]*(s*p[7])*log(abs(s*p[7])),
shepardBneg = 1 - (1 - (s*p[7])^2 + p[8]*(s*p[7])*log(abs(s*p[7]))),
shepardD = 1 - p[8]*(s*p[7]) + (s*p[7])^2,
shepardDneg = 1 - (1 - p[8]*(s*p[7]) + (s*p[7])^2),
shepardE = 1 - p[8]*s*p[7] + p[8]*(s*p[7])^2 - (s*p[7])^3,
shepardEneg = 1 - (1 - p[8]*s*p[7] + p[8]*(s*p[7])^2 - (s*p[7])^3),
shepardF = exp(-p[7]*s - p[8]),
shepardFneg = 1 - exp(-p[7]*s - p[8])
)
}

## QDM objective function
objfun <- function(p, psi, estimfun = c("minchi2", "ols", "wls"), ...){
yhat <- outer(psi\$x, psi\$y, function(x, y) qdmfun(x, y, p, ...))

estimfun <- match.arg(estimfun)
out <- switch(EXPR = estimfun,
minchi2 =
sum(na.omit(as.vector((psi\$freq - psi\$ntrials*yhat)^2 /
(psi\$ntrials*yhat * (1 - yhat))))),
ols =
sum((psi\$prob - yhat)^2),
wls = {
weights <- 1 / (psi\$prob*(1 - (psi\$prob - 0.01)))
sum(na.omit(as.vector(weights*((psi\$freq - psi\$ntrials*yhat)^2))))
}
)
#yhat[is.na(yhat)] <- 0.5
if (any(yhat < 0) | any(yhat > 1)) Inf else out
}

## QDM user interface
qdm <- function(psi, start, respfun = c("logistic", "guessing", "gumbel",
"gompertz", "weibull", "cauchy", "shepardA", "shepardAneg",
"shepardB", "shepardBneg", "shepardD", "shepardDneg",
"shepardE", "shepardEneg", "shepardF", "shepardFneg"),
bias = 0, estimfun = c("minchi2", "ols", "wls"),
optimizer = c("optim", "nlm"), optimargs = list()){

stopifnot(class(psi) == "psi")
respfun <- match.arg(respfun)
estimfun <- match.arg(estimfun)

## Optimize
optimizer <- match.arg(optimizer)
if (optimizer == "nlm") {
optArgs <- list(f=objfun, p=start, psi=psi, estimfun=estimfun,
response=respfun, bias=bias)
optArgs <- c(optArgs, as.list(optimargs))
optimout <- do.call(nlm, optArgs)
coefficients <- optimout\$estimate
} else {  # optim()
optArgs <- list(par=start, fn=objfun, psi=psi, estimfun=estimfun,
response=respfun, bias=bias)
optArgs <- c(optArgs, as.list(optimargs))
optimout <- do.call(optim, optArgs)
coefficients <- optimout\$par
}

retval = list(optimout=optimout, coefficients=coefficients, psi=psi,
respfun=respfun, bias=bias)
class(retval) <- "qdm"
retval
}
# TODO Give warning if optimizer does not give new parameters?

print.qdm <- function(x, digits = max(3L, getOption("digits") - 3L),
...){
cat("Coefficients:\n")
print.default(format(coef(x), digits = digits), print.gap = 2L,
quote = FALSE)
cat("\n")
invisible(x)
}

## FIX ME: args to qdmfun
predict.qdm <- function(object, x = object\$psi\$x, y = object\$psi\$y,
respfun = object\$respfun, bias = object\$bias, ...){
yhat <- outer(x, y, function(x, y) qdmfun(x, y, coef(object), respfun, bias))
dimnames(yhat) <- list(x, y)
yhat
}

persp.qdm <- function(x, col="gray", zlim=0:1, phi=10, theta=-25,
xlab="OA1", ylab="OA2", zlab="Predictions", ...){
persp(x\$psi\$x, x\$psi\$y, predict(x), col=col, zlim=zlim, phi=phi,
theta=theta, xlab=xlab, ylab=ylab, zlab=zlab, ...)
}
```

## Try the qdm package in your browser

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

qdm documentation built on May 30, 2017, 4:05 a.m.