# R/binom_gl.R In modelfree: Model-free estimation of a psychometric function

#### Documented in binom_gl

```binom_gl <-function( r, m, x, link, p, K,
initval) {
#
# THIS IS AN INTERNAL FUNCTION: USE BINOM_LIMS FOR BEST RESULTS
#
# Maximum likelihood estimates of the parameters of the psychometric function
# with guessing and lapsing rates. The estimated parameters for the
# linear part are in vector 'b' and the estimated rates are in 'guessing' and 'lapsing'.
#
# INPUT
#
# r    - number of successes at points x
# m    - number of trials at points x
# x    - stimulus levels
# p       - degree of the polynomial
# K       - power parameter for Weibull and reverse Weibull link
# initval - initial values for guessing and lapsing
#
# OUTPUT
#
# Object with 4 components:
# b - vector of estiamted coefficients for the linear part
# guessing - estimated guessing rate
# lapsing - estimated lapsing rate
# fit - glm object to be used in evaluation of fitted values

# LIKELIHOOD
likfun <- function( gl ) {

gl <- 1 / ( 1 + exp( -gl ) )
guessing <- min(gl[1], 1-gl[2] )
lapsing <- min( gl[2], 1-gl[1] )

# fit
fit <- glm( glmformula, data = glmdata, weights = m,
family = binomial( eval( call( linkfun, guessing, lapsing ) ) ) );

}
else {
fit <- glm( glmformula, data = glmdata, weights = m,
family = binomial( eval( call( linkfun, K, guessing, lapsing ) ) ) );
}

# FITTED PROBABILITIES
fitted <- as.numeric( predict( fit, type = "response" ) );
fitted[which( fitted <= guessing )] <- guessing + .Machine\$double.eps;
fitted[which( fitted >= 1-lapsing )] <- 1-lapsing - .Machine\$double.eps;

return( c( -( t( r ) %*% log( fitted ) + t( m - r ) %*%
log( 1 - fitted ) ) ) );
}

# MAIN PROGRAM
initval <- log( initval / ( 1 - initval ) );

# GLM settings
glmdata <- data.frame( cbind( r/m , m , x ) );
names( glmdata ) <- c( "resp", "m", "x" );

# formula
glmformula <- c( "resp ~ x" );
if( p > 1 ) {
for( pp in 2:p ) {
glmformula <- paste( glmformula, " + I(x^", pp,")", sep = "");
}
}
fit <- NULL;

# GLM fit
}else{
}

suppressWarnings( gl <- optim( initval, likfun )\$par );

guessing <- 1 / ( 1 + exp( -gl[1] ) );
lapsing <- 1 / ( 1 + exp( -gl[2] ) );

tmpglm <- glm( glmformula, data = glmdata, weights = m,
family = binomial( eval( call( linkfun, guessing, lapsing ) ) ) );
}
else {
tmpglm <- glm( glmformula, data = glmdata, weights = m,
family = binomial( eval( call( linkfun, K, guessing, lapsing ) ) ) );
}

b <- tmpglm\$coeff

tmpglm\$df.residual <- length(x) - (p + 1) - 2

fit <- tmpglm

value <- NULL
value\$guessing <- guessing
value\$lapsing <- lapsing
value\$b <- b
value\$fit <- fit

return( value );
}
```

## Try the modelfree package in your browser

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

modelfree documentation built on May 31, 2017, 4:28 a.m.