em_bme_3plm.c File Reference

Functions to estimate the item parameters by MMLE/EM-BME. More...

#include "libirt.h"
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_multiroots.h>
#include <gsl/gsl_linalg.h>

Data Structures
struct	bme_3plm_struct
	Used to passed extra parameter to bme_3plm_fdfdf2. More...

Functions
double	prob_3plm (double ability, double slope, double threshold, double asymptote)
	Compute the 3 parameter logistic function.
void	probs_3plm (gsl_vector *slopes, gsl_vector *thresholds, gsl_vector *asymptotes, gsl_vector *quad_points, gsl_matrix *probs)
	Compute the response functions for a 3PLM.
int	bme_3plm_fdfdf2 (const gsl_vector *par_3plm, void *params, double *f, gsl_vector *df, gsl_matrix *df2)
	Compute the log likelihood, gradient and Hessian of the item parameter.
int	bme_3plm_fdf (const gsl_vector *par_3plm, void *params, double *f, gsl_vector *df)
	Compute the log likelihood and gradient of the item parameter.
int	bme_3plm_dfdf2 (const gsl_vector *par_3plm, void *params, gsl_vector *df, gsl_matrix *df2)
	Compute the log likelihood gradient and Hessian of the item parameter.
double	bme_3plm_f (const gsl_vector *par_3plm, void *params)
	Compute the log likelihood of the item parameter.
int	bme_3plm_df (const gsl_vector *par_3plm, void *params, gsl_vector *df)
	Compute the log likelihood gradient of the item parameter.
int	bme_3plm_df2 (const gsl_vector *par_3plm, void *params, gsl_matrix *df2)
	Compute the log likelihood Hessian of the item parameter.
int	bme_3plm (int nbr_par, int max_iter, double prec, gsl_matrix *quad_freqs, gsl_vector *quad_points, gsl_vector *quad_sizes, int slope_prior, double slope_mean, double slope_dev, int thresh_prior, double thresh_mean, double thresh_dev, int asymp_prior, double asymp_mean, double asymp_weight, gsl_vector *slopes, gsl_vector *thresholds, gsl_vector *asymptotes, gsl_vector *slopes_stddev, gsl_vector *thresh_stddev, gsl_vector *asymp_stddev, gsl_vector_int *ignore, gsl_vector_int *notconverge, double *mllk)
	Does the maximization step of the EM algorithm to estimate the parameters of a 2PLM or 3PLM by MMLE or more generally (if a prior is used) by BME.
int	em_bme_3plm (int nbr_par, int max_em_iter, int max_nr_iter, double prec, gsl_matrix_int *patterns, gsl_vector *counts, gsl_vector *quad_points, gsl_vector *quad_weights, int slope_prior, double slope_mean, double slope_dev, int thresh_prior, double thresh_mean, double thresh_dev, int asymp_prior, double asymp_mean, double asymp_weight, gsl_vector *slopes, gsl_vector *thresholds, gsl_vector *asymptotes, gsl_vector *slopes_stddev, gsl_vector *thresh_stddev, gsl_vector *asymp_stddev, gsl_vector_int *ignore, int *nbr_notconverge, gsl_vector_int *notconverge, int adjust_weights)
	Estimate the parameters of a 2PLM or 3PLM by MMLE or more generally (if a prior is used) by BME.

Detailed Description

Functions to estimate the item parameters by MMLE/EM-BME.

Author

Stephane Germain germs.nosp@m.te@g.nosp@m.mail..nosp@m.com

The overall objectif is to find the ICC (item characteristic curves) maximizing the ML (marginal likelihood). An EM (expectation-maximization) iterative algorithm is used.

1. A grid of ability levels has to be fixed. Something like 32 values from -4 to 4 will do. Those are called the quadrature classes in the code and can be generated with the function "quadrature".

1. A first approximation of the ICC has to be available. Almost anything will do. For example the slopes can be initialized to 1, the thresholds to 0, and the asymptotes (if used) to 0.1.

1. For each pattern (a response vector from one subject) the a posteriori probabilities of being in each quadrature classes is computed by the function "posteriors".

1. The expected number of subject in each quadrature classes (quad_sizes), and for each item the expected number of subject in each quadrature classes having a success at this item (quad_freqs) are computed by the function "frequencies".

1. Once these quantities are assumed to be known, the log-likelihood can be maximized independantly for each item. The maximization is done by a root finding algorithm in the function "bme_3plm". A variant of Newton-Raphson from the gsl library is used. For that we must have a function giving the firsts (gradient) and seconds (hessian) derivatives of the log-likelihood by the item parameters, those are computed in "bme_3plm_fdfdf2".

1. Steps 3-5 are repeated until convergence is achieved.

Function Documentation

double prob_3plm	(	double	ability,
		double	slope,
		double	threshold,
		double	asymptote
	)

Compute the 3 parameter logistic function.

Parameters

[in]	ability
[in]	slope
[in]	threshold
[in]	asymptote

Returns

The 3PLM evaluated at ability.

void probs_3plm	(	gsl_vector *	slopes,
		gsl_vector *	thresholds,
		gsl_vector *	asymptotes,
		gsl_vector *	quad_points,
		gsl_matrix *	probs
	)

Compute the response functions for a 3PLM.

Parameters

[in]	slopes	A vector(items) with the slope parameters of each item.
[in]	thresholds	A vector(items) with the threshold parameters of each item.
[in]	asymptotes	A vector(items) with the asymptote parameters of each item. Can be NULL for a 2PLM.
[in]	quad_points	A vector(classes) with the middle points of each quadrature class.
[out]	probs	A matrix(items x classes) with the response functions.

Warning

The memory for probs should be allocated before.

Todo:

Stddev of the probs

int bme_3plm_fdfdf2	(	const gsl_vector *	par_3plm,
		void *	params,
		double *	f,
		gsl_vector *	df,
		gsl_matrix *	df2
	)

Compute the log likelihood, gradient and Hessian of the item parameter.

Parameters

[in]	par_3plm	The 3, 2 or 1 parameters to the 3PLM.
[in]	params	The extra parameter to passes to the function.
[out]	f	The log likelihood.
[out]	df	The gradient of the log likelihood.
[out]	df2	The Hessian of the log likelihood.

This function is not used directly by the root finding functions, but by others functions that comply with the gsl.

Returns

GSL_SUCCESS for success.

int bme_3plm_fdf	(	const gsl_vector *	par_3plm,
		void *	params,
		double *	f,
		gsl_vector *	df
	)

Compute the log likelihood and gradient of the item parameter.

Parameters

[in]	par_3plm	The 3, 2 or 1 parameters to the 3PLM.
[in]	params	The extra parameter to passes to the function.
[out]	f	The log likelihood.
[out]	df	The gradient of the log likelihood.

This function is just a wrapper around bme_3plmfdfdf2 to be used by the root finding functions in the gsl.

int bme_3plm_dfdf2	(	const gsl_vector *	par_3plm,
		void *	params,
		gsl_vector *	df,
		gsl_matrix *	df2
	)

Compute the log likelihood gradient and Hessian of the item parameter.

Parameters

[in]	par_3plm	The 3, 2 or 1 parameters to the 3PLM.
[in]	params	The extra parameter to passes to the function.
[out]	df	The gradient of the log likelihood.
[out]	df2	The Hessian of the log likelihood.

This function is just a wrapper around bme_3plmfdfdf2 to be used by the root finding functions in the gsl.

Returns

GSL_SUCCESS for success.

double bme_3plm_f	(	const gsl_vector *	par_3plm,
		void *	params
	)

Compute the log likelihood of the item parameter.

Parameters

[in]	par_3plm	The 3, 2 or 1 parameters to the 3PLM.
[in]	params	The extra parameter to passes to the function.

This function is just a wrapper around bme_3plmfdfdf2 to be used by the root finding functions in the gsl.

Returns

The log likelihood.

int bme_3plm_df	(	const gsl_vector *	par_3plm,
		void *	params,
		gsl_vector *	df
	)

Compute the log likelihood gradient of the item parameter.

Parameters

[in]	par_3plm	The 3, 2 or 1 parameters to the 3PLM.
[in]	params	The extra parameter to passes to the function.
[out]	df	The gradient of the log likelihood.

This function is just a wrapper around bme_3plmfdfdf2 to be used by the root finding functions in the gsl.

Returns

GSL_SUCCESS for success.

int bme_3plm_df2	(	const gsl_vector *	par_3plm,
		void *	params,
		gsl_matrix *	df2
	)

Compute the log likelihood Hessian of the item parameter.

Parameters

[in]	par_3plm	The 3, 2 or 1 parameters to the 3PLM.
[in]	params	The extra parameter to passes to the function.
[out]	df2	The Hessian of the log likelihood.

This function is just a wrapper around bme_3plmfdfdf2 to be used by the root finding functions in the gsl.

Returns

GSL_SUCCESS for success.

int bme_3plm	(	int	nbr_par,
		int	max_iter,
		double	prec,
		gsl_matrix *	quad_freqs,
		gsl_vector *	quad_points,
		gsl_vector *	quad_sizes,
		int	slope_prior,
		double	slope_mean,
		double	slope_dev,
		int	thresh_prior,
		double	thresh_mean,
		double	thresh_dev,
		int	asymp_prior,
		double	asymp_mean,
		double	asymp_weight,
		gsl_vector *	slopes,
		gsl_vector *	thresholds,
		gsl_vector *	asymptotes,
		gsl_vector *	slopes_stddev,
		gsl_vector *	thresh_stddev,
		gsl_vector *	asymp_stddev,
		gsl_vector_int *	ignore,
		gsl_vector_int *	notconverge,
		double *	mllk
	)

Does the maximization step of the EM algorithm to estimate the parameters of a 2PLM or 3PLM by MMLE or more generally (if a prior is used) by BME.

Parameters

[in]	nbr_par	The number of parameter (1, 2 or 3).
[in]	max_iter	The maximum number of Newton iterations performed for each item.
[in]	prec	The desired precision of each parameter estimate.
[in]	quad_freqs	A matrix(items x classes) with the expected number of subjects in the class having a success at the item.
[in]	quad_points	A vector(classes) with the middle points of each quadrature class.
[in]	quad_sizes	A vector (classes) with the expected number of subjects in the class.
[in]	slope_prior	Controls whether to use a lognormal(slope_mean, slope_dev) prior on the slope.
[in]	slope_mean	The mean of the lognormal prior on the slope.
[in]	slope_dev	The standard deviation of the lognormal prior on the slope.
[in]	thresh_prior	Controls whether to use a normal(thresh_mean, thresh_dev) prior on the threshold
[in]	thresh_mean	The mean of the normal prior on the threshold.
[in]	thresh_dev	The standard deviation of the normal prior on the threshold.
[in]	asymp_prior	Controls whether to use a beta(asymp_mean, asymp_weight) prior on the asymptote
[in]	asymp_mean	The mean of the beta prior on the asymptote. The alpha and beta parameter are alpha = mean * weight + 1 and beta = (1 - mean) * weight + 1).
[in]	asymp_weight	The weight of the beta prior on the asymptote. The alpha and beta parameter are alpha = mean * weight + 1 and beta = (1 - mean) * weight + 1).
[in,out]	slopes	A vector(items) with the estimated slope parameters. They should be initialize first.
[out]	slopes_stddev	A vector(items) with the standard errors of the estimated slopes.
[in,out]	thresholds	A vector(items) with the estimated threshold parameters. They should be initialize first.
[out]	thresh_stddev	A vector(items) with the standard errors of the estimated thresholds.
[in,out]	asymptotes	A vector(items) with the estimated asymptote parameters. They should be initialize first.
[out]	asymp_stddev	A vector(items) with the standard errors of the estimated asymptotes.
[in]	ignore	A vector(items) of ignore flag.
[out]	notconverge	A vector(items) of flag set for the items that didn't converged.
[out]	mllk	The maximum log likelihood.

Returns

The number of item that did not converged.

Warning

The memory for the output parameters should be allocated before.

int em_bme_3plm	(	int	nbr_par,
		int	max_em_iter,
		int	max_nr_iter,
		double	prec,
		gsl_matrix_int *	patterns,
		gsl_vector *	counts,
		gsl_vector *	quad_points,
		gsl_vector *	quad_weights,
		int	slope_prior,
		double	slope_mean,
		double	slope_dev,
		int	thresh_prior,
		double	thresh_mean,
		double	thresh_dev,
		int	asymp_prior,
		double	asymp_mean,
		double	asymp_weight,
		gsl_vector *	slopes,
		gsl_vector *	thresholds,
		gsl_vector *	asymptotes,
		gsl_vector *	slopes_stddev,
		gsl_vector *	thresh_stddev,
		gsl_vector *	asymp_stddev,
		gsl_vector_int *	ignore,
		int *	nbr_notconverge,
		gsl_vector_int *	notconverge,
		int	adjust_weights
	)

Estimate the parameters of a 2PLM or 3PLM by MMLE or more generally (if a prior is used) by BME.

Parameters

[in]	nbr_par	The number of parameter (1, 2 or 3).
[in]	max_em_iter	The maximum number of EM iterations. At least 20 iteration are made.
[in]	max_nr_iter	The maximum number of Newton iterations performed for each item at each EM iteration.
[in]	prec	The relative change in the likelihood to stop the EM algorithm. This value divided by 10 is also the desired precision of each parameter estimate.
[in]	patterns	A matrix(patterns x items) of binary responses.
[in]	counts	A vector(patterns) with the count of each pattern. If NULL the counts are assumed to be all 1.
[in]	quad_points	A vector(classes) with the middle points of each quadrature class.
[in]	quad_weights	A vector(classes) with the prior weights of each quadrature class.
[in]	slope_prior	Controls whether to use a lognormal(slope_mean, slope_dev) prior on the slope.
[in]	slope_mean	The mean of the lognormal prior on the slope.
[in]	slope_dev	The standard deviation of the lognormal prior on the slope.
[in]	thresh_prior	Controls whether to use a normal(thresh_mean, thresh_dev) prior on the threshold
[in]	thresh_mean	The mean of the normal prior on the threshold.
[in]	thresh_dev	The standard deviation of the normal prior on the threshold.
[in]	asymp_prior	Controls whether to use a beta(asymp_mean, asymp_weight) prior on the asymptote
[in]	asymp_mean	The mean of the beta prior on the asymptote. The alpha and beta parameter are alpha = mean * weight + 1 and beta = (1 - mean) * weight + 1).
[in]	asymp_weight	The weight of the beta prior on the asymptote. The alpha and beta parameter are alpha = mean * weight + 1 and beta = (1 - mean) * weight + 1).
[in,out]	slopes	A vector(items) with the estimated slope parameters. They should be initialize first.
[out]	slopes_stddev	A vector(items) with the standard errors of the estimated slopes.
[in,out]	thresholds	A vector(items) with the estimated threshold parameters. They should be initialize first.
[out]	thresh_stddev	A vector(items) with the standard errors of the estimated thresholds.
[in,out]	asymptotes	A vector(items) with the estimated asymptote parameters. They should be initialize first.
[out]	asymp_stddev	A vector(items) with the standard errors of the estimated asymptotes.
[in]	ignore	A vector(items) of ignore flag.
[out]	nbr_notconverge	The number of items that didn't converged.
[out]	notconverge	A vector(items) of flag set for the items that didn't converged.
[in]	adjust_weights	Controls whether adjust the quadrature weights after each iteration.

Returns

1 if the relative change in the maximum log likelihood was less than prec else 0.

Warning

The memory for the output parameters should be allocated before.

Todo:

Check for the biggest change in the EM iterations and add a stop criterion.