2
$\begingroup$

I am trying to implement a 2-PL dichotomous IRT Model for my dataset from scratch in Python. Here is my code so far:


def get_data_matrix(num_students, num_items):
    data = np.random.randint(0, 2, size=(num_students, num_items))
    return data



# Define the logistic function
def logistic_func(theta, b, a):
    return 1 / (1 + np.exp(-a * (theta - b)))


# Define the log-likelihood function for the 2PL model
def log_likelihood(params, response_matrix):
    num_students, num_questions = response_matrix.shape
    theta = params[:num_students]
    b = params[num_students:num_students + num_questions]
    a = params[num_students + num_questions:]

    ll = 0
    for i in range(num_students):
        for j in range(num_questions):
            ll += response_matrix[i, j] * np.log(logistic_func(theta[i], b[j], a[j]) + 1e-9)
            ll += (1 - response_matrix[i, j]) * np.log(1 - logistic_func(theta[i], b[j], a[j]) + 1e-9)

    return -ll

# Define the CMLE function for estimating theta
def cmle(response_matrix, b_estimated, a_estimated):
    num_students, num_items = response_matrix.shape

    def objective_function(theta):
        ll = 0
        for i in range(num_students):
           for j in range(num_items):
                ll += response_matrix[i, j] * np.log(logistic_func(theta[i], b_estimated[j], a_estimated[j]) + 1e-9)
                ll += (1 - response_matrix[i, j]) * np.log(1 - logistic_func(theta[i], b_estimated[j], a_estimated[j]) + 1e-9)

        return -ll

    initial_theta = np.zeros(num_students)  # Initial theta values
    bounds = [(-3, 3)] * num_students  # Bounds for theta estimation

    result = minimize(objective_function, initial_theta, bounds=bounds, method='L-BFGS-B')

    theta_estimated = result.x

    return theta_estimated




if __name__ == "__main__":

    # Generate example data
    num_students = 400
    num_questions = 20
    
    # Assuming you have responses and design_matrix as your data inputs
    # This is a random dataset -- will be replaced by the original data
    data = get_matrix(num_students, num_questions)

    initial_theta = np.zeros(num_students)
    initial_b = np.zeros(num_questions)
    initial_a = np.ones(num_questions)
    initial_params = np.concatenate((initial_theta, initial_b, initial_a))

    # Define the bounds for the parameters
    bounds = [(-3, 3)] * num_students + [(-3, 3)] * num_questions + [(-3, 3)] * num_questions

    # Perform maximum likelihood estimation
    result = minimize(log_likelihood, initial_params, args=(data,), bounds=bounds, method='L-BFGS-B')


    # Extract the estimated parameters
    estimated_b = result.x[num_students:num_students + num_questions]
    estimated_a = result.x[num_students + num_questions:]

    # Fit a CMLE again to estimate theta values again
    estimated_theta = cmle(data, estimated_b, estimated_a)

    # Compute the infit and outfit values
    data_succ_proportions = np.sum(data, axis=0)/num_students

    predicted_data = logistic_func(estimated_theta.reshape(-1, 1), estimated_b, estimated_a)
    # Compute the information matrix
    information_matrix = np.zeros_like(data, dtype=float)
    for i in range(num_questions):
        for j in range(num_students):
            p = predicted_data[j][i]
            information_matrix[j][i] = p * (1 - p)

    outfit = np.sum((data - predicted_data) ** 2/ information_matrix, axis=0) / (num_students)
    infit = np.sum((data - predicted_data) ** 2, axis=0) / np.sum(information_matrix, axis=0)



    # Compute the range of theta values
    theta_range = np.linspace(-4, 4, num=100)

    # Plot the ICC for each item
    for j in range(num_questions):
        p = logistic_func(theta_range, estimated_b[j], estimated_a[j])
        print("Item " + str(j) + ":" + str(estimated_b[j]) + " " + str(estimated_a[j]))
        plt.plot(theta_range, p, label=f"Item {j + 1}")


    plt.xlabel("Theta")
    plt.ylabel("Probability")
    plt.title("Item Characteristic Curves (ICC)")
    # Legend customization
    legend = plt.legend(fontsize='x-small')
    legend.get_frame().set_linewidth(0.5)
    legend.get_frame().set_edgecolor('black')
    legend.get_frame().set_alpha(0.7)
    plt.show()

    # Print the infit and outfit values
    print("Data success:")
    for item_idx, value in enumerate(data_succ_proportions):
        print(f"Item {item_idx + 1}: {value}")

    print("Infit:")
    for item_idx, value in enumerate(infit):
        print(f"Item {item_idx + 1}: {value}")

    print("\nOutfit:")
    for item_idx, value in enumerate(outfit):
        print(f"Item {item_idx + 1}: {value}")

However, when I use a in-built library in a software like Jamovi to compute the same 2-PL IRT model on my dataset, I get a very different set of b_estimates, a_estimates and even InFit and OutFit values. For each item, the b_estimates at least exhibit similar trends but the infit and outfit values do not even show the same trends.

My questions are the following:

  • Am I implementing 2-PL dichotomous IRT Model correctly and computing the infit and outfit stats correctly?
  • On JAMOVI (snowirt module) they say the Marginal Maximum Likelihood estimate was used to obtain the result tables. Is my implementation also for Marginal Maximum Likelihood estimate?
$\endgroup$

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.