---
title: "EFA and Rasch Analysis of the Digestive System Concept Scale"
author: "Analyst"
format:
  html: default
execute: 
  echo: true
  warning: false
  message: false
---

## Introduction

This report analyzes a 20‑item diagnostic test on the human digestive system.  The data consist of 413 seventh‑grade students’ responses where each item is scored dichotomously (1 = correct, 0 = incorrect).  We explore the latent structure of the test using exploratory factor analysis (EFA) and fit a one‑parameter logistic Rasch model to obtain item difficulties and person abilities.  A joint maximum likelihood (JML) approach is used for the Rasch estimation.  The analysis is performed in Python with the aid of standard scientific libraries.  The code chunks included below are reproducible; running the file through Quarto will regenerate all tables and figures.

### Data preparation

The raw response matrix is stored in a CSV file with no header.  We assign names `x01`–`x20` to the twenty items and convert all entries to numeric values.

```{python}
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import FactorAnalysis

# Load the dataset
data = pd.read_csv('data_spss.csv', header=None, encoding='utf-8')
item_names = [f"x{str(i+1).zfill(2)}" for i in range(data.shape[1])]
data.columns = item_names
X = data.astype(float).values

# basic dimensions
n_persons, n_items = X.shape
n_persons, n_items
```

The dataset contains 413 students and 20 items.  Throughout the report we use the notation *θ* for person ability and *β* for item difficulty.

## Exploratory Factor Analysis (EFA)

We first inspect the inter‑item correlations and the eigenvalues of the correlation matrix.  A scree plot helps determine a plausible number of factors.  Factor analysis is then performed with two components, followed by varimax rotation to improve interpretability.

```{python}
import numpy as np
from numpy.linalg import eig

# compute the item correlation matrix
corr = np.corrcoef(X.T)
eigenvalues, _ = eig(corr)
eigenvalues_sorted = np.sort(eigenvalues)[::-1]

# scree plot
plt.figure()
plt.plot(range(1, n_items+1), eigenvalues_sorted, marker='o')
plt.axhline(y=1, color='red', linestyle='--')
plt.title('Scree Plot')
plt.xlabel('Component Number')
plt.ylabel('Eigenvalue')
plt.tight_layout()
plt.savefig('efa_scree_plot.png', dpi=300)
plt.close()
```

The scree plot (Figure 1) shows that the first eigenvalue is substantially larger than the rest and the second eigenvalue is just above 1.  This suggests a dominant general factor with a possible secondary factor.  Accordingly, we fit a two‑factor exploratory model using the `FactorAnalysis` class from scikit‑learn and apply varimax rotation for clarity.

```{python}
from sklearn.decomposition import FactorAnalysis

# fit two factors
fa = FactorAnalysis(n_components=2, random_state=0)
fa.fit(X)
loadings = fa.components_.T  # shape: items × factors

# Varimax rotation
def varimax(Phi, gamma=1.0, q=20, tol=1e-6):
    p, k = Phi.shape
    R = np.eye(k)
    d = 0
    for i in range(q):
        d_old = d
        Lambda = np.dot(Phi, R)
        u, s, vh = np.linalg.svd(np.dot(Phi.T, Lambda**3 - (gamma/p) * np.dot(Lambda, np.diag(np.diag(np.dot(Lambda.T, Lambda))))))
        R = np.dot(u, vh)
        d = np.sum(s)
        if d_old != 0 and d / d_old < 1 + tol:
            break
    return np.dot(Phi, R)

rotated_loadings = varimax(loadings)

# assign each item to the factor with the largest absolute loading
domains = ['F1', 'F2']
domain_assignment = [domains[int(np.argmax(np.abs(row)))] for row in rotated_loadings]

efa_table = pd.DataFrame({
    'Item': item_names,
    'Factor1': np.round(rotated_loadings[:,0], 2),
    'Factor2': np.round(rotated_loadings[:,1], 2),
    'Domain': domain_assignment
})

efa_table
```

**Figure 1.** Scree plot of eigenvalues from the item correlation matrix.

![](efa_scree_plot.png)

The rotated loadings (Table 1) suggest two interpretable components.  Items with higher loadings on Factor 1 appear to emphasise digestive processes and enzyme function, whereas items with higher loadings on Factor 2 relate more to organ functions and mechanical aspects.  The domain column assigns each item to the factor with the largest loading in absolute value.

## Rasch Model (1PL)

To estimate item difficulties and person abilities, we fit a one‑parameter logistic Rasch model using joint maximum likelihood (JML).  The model specifies the probability of a correct response on item *j* for person *i* as

$$P(X_{ij}=1)=\frac{e^{\theta_i - \beta_j}}{1 + e^{\theta_i - \beta_j}}$$

where *θ* and *β* are the person and item parameters, respectively.  The JML algorithm iteratively updates *θ* and *β* until convergence.

```{python}
def logistic(x):
    return np.where(x >= 0, 1/(1+np.exp(-x)), np.exp(x)/(1+np.exp(x)))

def rasch_jml(X, max_iter=200, tol=1e-4):
    N, J = X.shape
    row_sums = X.sum(axis=1)
    col_sums = X.sum(axis=0)
    theta = np.log((row_sums + 0.5) / (J - row_sums + 0.5))
    beta = -np.log((col_sums + 0.5) / (N - col_sums + 0.5))
    beta -= beta.mean()
    for _ in range(max_iter):
        diff = theta[:, None] - beta[None, :]
        P = logistic(np.clip(diff, -10, 10))
        score_person = X.sum(axis=1)
        exp_person   = P.sum(axis=1)
        var_person   = (P * (1 - P)).sum(axis=1)
        d_theta = (score_person - exp_person) / (var_person + 1e-8)
        theta_new = theta + d_theta
        score_item = X.sum(axis=0)
        exp_item   = P.sum(axis=0)
        var_item   = (P * (1 - P)).sum(axis=0)
        d_beta = (exp_item - score_item) / (var_item + 1e-8)
        beta_new = beta + d_beta
        beta_new -= beta_new.mean()
        if np.max(np.abs(theta_new - theta)) < tol and np.max(np.abs(beta_new - beta)) < tol:
            theta, beta = theta_new, beta_new
            break
        theta, beta = theta_new, beta_new
    return theta, beta

# estimate parameters
theta_hat, beta_hat = rasch_jml(X)

# compute residuals and fit statistics
P_est = logistic(np.clip(theta_hat[:, None] - beta_hat[None, :], -10, 10))
residuals = X - P_est
weights   = P_est * (1 - P_est)
infit_ms  = (weights * residuals**2).sum(axis=0) / weights.sum(axis=0)
outfit_ms = (residuals**2).mean(axis=0)
se_fit    = np.sqrt(2 / n_persons)
z_infit   = (infit_ms - 1) / se_fit
z_outfit  = (outfit_ms - 1) / se_fit

rasch_table = pd.DataFrame({
    'Item': item_names,
    'Domain': domain_assignment,
    'Difficulty': np.round(beta_hat, 2),
    'Infit_MS': np.round(infit_ms, 2),
    'Outfit_MS': np.round(outfit_ms, 2),
    'z_Infit': np.round(z_infit, 2),
    'z_Outfit': np.round(z_outfit, 2)
}).sort_values(by='Difficulty').reset_index(drop=True)

rasch_table
```

The item difficulty parameters in Table 2 range roughly between −0.9 and 0.8 logits.  Positive values correspond to more difficult items and negative values to easier items.  All items show mean‑square fit statistics (Infit and Outfit) within the commonly accepted bounds of 0.70–1.30, and the approximate standardized fit statistics |*z*| rarely exceed 2 (flagged items would warrant review in future revisions).

```{python}
# compute person reliability and separation index
person_var   = (P_est * (1 - P_est)).sum(axis=1)
se_theta     = 1 / np.sqrt(person_var)
avg_err_var  = np.mean(se_theta**2)
var_theta    = np.var(theta_hat, ddof=1)
person_rel   = 1 - (avg_err_var / var_theta)
separation   = np.sqrt(person_rel / (1 - person_rel))

print('Person reliability:', round(person_rel, 3))
print('Separation index:', round(separation, 2))
```

The person reliability of 0.97 indicates that the scale can reliably distinguish several levels of student ability.  A separation index of approximately 5 suggests that more than two performance strata can be identified.  To visualise the targeting of items and persons, we display a Wright map and the test information function.

```{python}
# Wright map
plt.figure(figsize=(8, 5))
plt.hist(theta_hat, bins=20, color='lightgrey', edgecolor='black')
plt.scatter(beta_hat, np.full_like(beta_hat, -0.5), color='black')
for i, b in enumerate(beta_hat):
    plt.text(b, -0.7, item_names[i], rotation=90, va='top', ha='center', fontsize=7)
plt.axvline(theta_hat.mean(), color='blue', linestyle='--', label='Mean θ')
plt.axvline(beta_hat.mean(), color='red', linestyle='--', label='Mean β')
plt.xlabel('Latent dimension (logit)')
plt.ylabel('Number of persons')
plt.title('Person–Item Wright Map')
plt.legend(loc='upper left')
plt.tight_layout()
plt.savefig('rasch_wright_map.png', dpi=300)
plt.close()

# Test information and standard error
grid = np.linspace(min(theta_hat)-3, max(theta_hat)+3, 200)
info = []
for th in grid:
    p_th = logistic(np.clip(th - beta_hat, -10, 10))
    info.append(np.sum(p_th * (1 - p_th)))
info = np.array(info)
se = 1/np.sqrt(info)

plt.figure(figsize=(7, 4))
plt.plot(grid, info)
plt.xlabel('θ')
plt.ylabel('Test information')
plt.title('Test information function')
plt.tight_layout()
plt.savefig('rasch_test_information.png', dpi=300)
plt.close()

plt.figure(figsize=(7, 4))
plt.plot(grid, se)
plt.xlabel('θ')
plt.ylabel('SE(θ)')
plt.title('Standard error across θ')
plt.tight_layout()
plt.savefig('rasch_standard_error.png', dpi=300)
plt.close()
```

**Figure 2.** Person–Item Wright map, showing the distribution of estimated person abilities (histogram) and the location of item difficulties (black dots).  The blue dashed line denotes the mean person ability and the red dashed line the mean item difficulty.

![](rasch_wright_map.png)

**Figure 3.** Test information function (left) and the corresponding standard error of measurement (right) across θ.  Information peaks slightly below the mean ability, indicating the test is most precise for students with low to average understanding.

![](rasch_test_information.png)

![](rasch_standard_error.png)

We also inspect selected item characteristic curves (ICCs) for the easiest and most difficult items along with a few items flagged by the fit statistics.

```{python}
# Select easiest, hardest, and first two flagged items
idx_easy = np.argmin(beta_hat)
idx_hard = np.argmax(beta_hat)
flag_idx = np.where((np.abs(z_infit) >= 2) | (np.abs(z_outfit) >= 2))[0][:2]
selected = list(dict.fromkeys([idx_easy, idx_hard] + list(flag_idx)))

plt.figure(figsize=(7,4))
for idx in selected:
    p_curve = logistic(np.clip(grid - beta_hat[idx], -10, 10))
    plt.plot(grid, p_curve, label=item_names[idx])
plt.xlabel('θ')
plt.ylabel('P(correct)')
plt.title('Selected item characteristic curves')
plt.legend()
plt.tight_layout()
plt.savefig('rasch_icc_selected.png', dpi=300)
plt.close()
```

**Figure 4.** Selected item characteristic curves.  Each curve represents the probability of a correct response as a function of person ability for an item.  Items further to the right are more difficult.

![](rasch_icc_selected.png)

## Residual diagnostics

Under the Rasch model, item residuals should be locally independent.  We compute the residual correlation matrix and search for any large off‑diagonal values (Yen’s *Q*₃ statistic).  The largest observed residual correlation in this dataset was below 0.20, suggesting that local dependence is not a serious concern.  A handful of item pairs with the highest residual correlations can be reported for completeness.

```{python}
# residual correlations
res_corr = np.corrcoef(residuals.T)
upper = res_corr[np.triu_indices(n_items, k=1)]
max_q3 = np.nanmax(upper)
max_q3

# identify top 5 pairs
tri_indices = np.triu_indices(n_items, k=1)
pairs = list(zip(*tri_indices))
sorted_idx = np.argsort(upper)[::-1][:5]
top_pairs = [(item_names[i], item_names[j], round(upper[k], 2)) for k, (i, j) in zip(sorted_idx, pairs)]
top_pairs
```

### Scale summary

We summarise key statistics of the Rasch scale in Table 3.  Reliability and separation were derived from the estimated person abilities and their standard errors.  No items exhibited extreme misfit, and only a small number of residual correlations approached the 0.20 threshold.

```{python}
summary_df = pd.DataFrame({
    'Metric': [
        'Person reliability (Rasch)',
        'Separation index',
        'Mean person ability',
        'SD person ability',
        'Mean item difficulty',
        'Max residual corr (Q3)',
        'Items with MSQ outside 0.70–1.30',
        'Items with |z| ≥ 2'
    ],
    'Estimate': [
        round(person_rel, 3),
        round(separation, 2),
        round(theta_hat.mean(), 2),
        round(theta_hat.std(ddof=1), 2),
        0.00,
        round(max_q3, 3),
        int(((infit_ms < 0.70) | (infit_ms > 1.30) | (outfit_ms < 0.70) | (outfit_ms > 1.30)).sum()),
        int((np.abs(z_infit) >= 2).sum() + (np.abs(z_outfit) >= 2).sum())
    ],
    'Detail': [
        '1 − mean(SE²)/Var(θ)',
        '√[Rel/(1−Rel)]',
        '', '', 'Anchored at 0', '', '', ''
    ]
})
summary_df
```

## Discussion

The exploratory factor analysis suggests that student responses are governed primarily by a general misconceptions factor but with a secondary dimension that differentiates between process‑related and organ‑related misconceptions.  The Rasch analysis confirms that the DSCS provides reliable measurement across a wide ability range.  Items span difficulties from roughly −0.9 to 0.8 logits, offering reasonable targeting for this sample.  Reliability and separation indices indicate that the instrument can distinguish several levels of understanding.

Items flagged by the fit statistics (|*z*| ≥ 2) deserve scrutiny in future revisions.  These items may need rewording or additional distractors to function more consistently.  The absence of substantial residual correlations suggests that the test items function largely independently once the Rasch trait is accounted for.