GitHub

LSAT: Item Response

This is the classic LSAT example from Volume 1 of the BUGS examples (see also the OpenBUGS write-up). Section 6 of the Law School Aptitude Test (LSAT) is a 5-item multiple choice test; students score 1 on each item for a correct answer and 0 otherwise, giving 32 possible response patterns. Bock and Lieberman (1970) present data on this test for 1000 students, recorded as the frequency of each of the R = 32 response patterns.

The question is how to separate the difficulty of each test item from the ability of each student. The data are analysed with the one-parameter Rasch model, a foundational item response model: the probability that student j answers item k correctly follows a logistic function of an item difficulty parameter alpha[k] and a latent ability theta[j], with abilities assumed normally distributed in the student population. The scale parameter beta (constrained to be positive) governs the spread of the ability distribution. Because the location of the difficulties is only identified relative to the mean ability (fixed at zero), the model also computes centred difficulties a[k] = alpha[k] - mean(alpha), which can be compared with the marginal maximum likelihood estimates of Bock and Aitkin (1981).

Model

Writing $p_{jk}$ for the probability that student $j$ answers item $k$ correctly, the model is

\[\begin{aligned} r_{jk} &\sim \text{Bernoulli}(p_{jk}) \\ \operatorname{logit}(p_{jk}) &= \beta\,\theta_j - \alpha_k, \qquad j = 1,\dots,1000;\ k = 1,\dots,5 \\ \theta_j &\sim \text{Normal}(0, 1) \end{aligned}\]

with vague normal priors on the item difficulties $\alpha_k$ and a flat prior on $(0, 1000)$ for $\beta$.

The data arrive as 32 aggregated response patterns rather than individual answers, so the first part of the program expands them: using the cumulative pattern counts culm, it assigns each of the 1000 students the binary response vector of their pattern. This deterministic data transformation is carried over directly from the original BUGS program. As in BUGS generally, the link-function form logit(p[j, k]) <- ... is written in Julia by applying the inverse link (logistic) on the right-hand side.

using JuliaBUGS

lsat = @bugs begin
    # Calculate individual (binary) responses to each test from Multinomial data
    for j in 1:culm[1]
        for k in 1:T
            r[j, k] = response[1, k]
        end
    end

    for i in 2:R
        for j in (culm[i - 1] + 1):culm[i]
            for k in 1:T
                r[j, k] = response[i, k]
            end
        end
    end

    # Rasch model
    for j in 1:N
        for k in 1:T
            p[j, k] = logistic(beta * theta[j] - alpha[k])
            r[j, k] ~ dbern(p[j, k])
        end
        theta[j] ~ dnorm(0, 1)
    end

    # Priors
    for k in 1:T
        alpha[k] ~ dnorm(0, 0.0001)
        a[k] = alpha[k] - mean(alpha[:])
    end
    beta ~ dunif(0, 1000)
end
BUGSModelDef:
begin
    for j = 1:culm[1]
        for k = 1:T
            r[j, k] = response[1, k]
        end
    end
    for i = 2:R
        for j = culm[i - 1] + 1:culm[i]
            for k = 1:T
                r[j, k] = response[i, k]
            end
        end
    end
    for j = 1:N
        for k = 1:T
            p[j, k] = logistic(beta * theta[j] - alpha[k])
            r[j, k] ~ dbern(p[j, k])
        end
        theta[j] ~ dnorm(0, 1)
    end
    for k = 1:T
        alpha[k] ~ dnorm(0, 0.0001)
        a[k] = alpha[k] - mean(alpha[:])
    end
    beta ~ dunif(0, 1000)
end

Data

The data are supplied as a NamedTuple. N = 1000 is the number of students, T = 5 the number of test items, and R = 32 the number of distinct response patterns. Each row of response is one pattern of five 0/1 answers, and culm gives the cumulative number of students whose answers match patterns up to and including that row.

data = (
    N = 1000,
    R = 32,
    T = 5,
    culm = [3, 9, 11, 22, 23, 24, 27, 31, 32, 40, 40, 56, 56, 59, 61, 76, 86, 115,
        129, 210, 213, 241, 256, 336, 352, 408, 429, 602, 613, 674, 702, 1000],
    response = [0 0 0 0 0
                0 0 0 0 1
                0 0 0 1 0
                0 0 0 1 1
                0 0 1 0 0
                0 0 1 0 1
                0 0 1 1 0
                0 0 1 1 1
                0 1 0 0 0
                0 1 0 0 1
                0 1 0 1 0
                0 1 0 1 1
                0 1 1 0 0
                0 1 1 0 1
                0 1 1 1 0
                0 1 1 1 1
                1 0 0 0 0
                1 0 0 0 1
                1 0 0 1 0
                1 0 0 1 1
                1 0 1 0 0
                1 0 1 0 1
                1 0 1 1 0
                1 0 1 1 1
                1 1 0 0 0
                1 1 0 0 1
                1 1 0 1 0
                1 1 0 1 1
                1 1 1 0 0
                1 1 1 0 1
                1 1 1 1 0
                1 1 1 1 1]
)

model = lsat(data)
BUGSModel (parameters are in transformed (unconstrained) space, with dimension 1006):

  Model parameters:
    alpha[5], alpha[4], alpha[3], alpha[2], alpha[1]
    beta
    theta[1000], theta[999], theta[998], theta[997], theta[996], theta[995], theta[994], theta[993], theta[992], theta[991], theta[990], theta[989], theta[988], theta[987], theta[986], theta[985], theta[984], theta[983], theta[982], theta[981], theta[980], theta[979], theta[978], theta[977], theta[976], theta[975], theta[974], theta[973], theta[972], theta[971], theta[970], theta[969], theta[968], theta[967], theta[966], theta[965], theta[964], theta[963], theta[962], theta[961], theta[960], theta[959], theta[958], theta[957], theta[956], theta[955], theta[954], theta[953], theta[952], theta[951], theta[950], theta[949], theta[948], theta[947], theta[946], theta[945], theta[944], theta[943], theta[942], theta[941], theta[940], theta[939], theta[938], theta[937], theta[936], theta[935], theta[934], theta[933], theta[932], theta[931], theta[930], theta[929], theta[928], theta[927], theta[926], theta[925], theta[924], theta[923], theta[922], theta[921], theta[920], theta[919], theta[918], theta[917], theta[916], theta[915], theta[914], theta[913], theta[912], theta[911], theta[910], theta[909], theta[908], theta[907], theta[906], theta[905], theta[904], theta[903], theta[902], theta[901], theta[900], theta[899], theta[898], theta[897], theta[896], theta[895], theta[894], theta[893], theta[892], theta[891], theta[890], theta[889], theta[888], theta[887], theta[886], theta[885], theta[884], theta[883], theta[882], theta[881], theta[880], theta[879], theta[878], theta[877], theta[876], theta[875], theta[874], theta[873], theta[872], theta[871], theta[870], theta[869], theta[868], theta[867], theta[866], theta[865], theta[864], theta[863], theta[862], theta[861], theta[860], theta[859], theta[858], theta[857], theta[856], theta[855], theta[854], theta[853], theta[852], theta[851], theta[850], theta[849], theta[848], theta[847], theta[846], theta[845], theta[844], theta[843], theta[842], theta[841], theta[840], theta[839], theta[838], theta[837], theta[836], theta[835], theta[834], theta[833], theta[832], theta[831], theta[830], theta[829], theta[828], theta[827], theta[826], theta[825], theta[824], theta[823], theta[822], theta[821], theta[820], theta[819], theta[818], theta[817], theta[816], theta[815], theta[814], theta[813], theta[812], theta[811], theta[810], theta[809], theta[808], theta[807], theta[806], theta[805], theta[804], theta[803], theta[802], theta[801], theta[800], theta[799], theta[798], theta[797], theta[796], theta[795], theta[794], theta[793], theta[792], theta[791], theta[790], theta[789], theta[788], theta[787], theta[786], theta[785], theta[784], theta[783], theta[782], theta[781], theta[780], theta[779], theta[778], theta[777], theta[776], theta[775], theta[774], theta[773], theta[772], theta[771], theta[770], theta[769], theta[768], theta[767], theta[766], theta[765], theta[764], theta[763], theta[762], theta[761], theta[760], theta[759], theta[758], theta[757], theta[756], theta[755], theta[754], theta[753], theta[752], theta[751], theta[750], theta[749], theta[748], theta[747], theta[746], theta[745], theta[744], theta[743], theta[742], theta[741], theta[740], theta[739], theta[738], theta[737], theta[736], theta[735], theta[734], theta[733], theta[732], theta[731], theta[730], theta[729], theta[728], theta[727], theta[726], theta[725], theta[724], theta[723], theta[722], theta[721], theta[720], theta[719], theta[718], theta[717], theta[716], theta[715], theta[714], theta[713], theta[712], theta[711], theta[710], theta[709], theta[708], theta[707], theta[706], theta[705], theta[704], theta[703], theta[702], theta[701], theta[700], theta[699], theta[698], theta[697], theta[696], theta[695], theta[694], theta[693], theta[692], theta[691], theta[690], theta[689], theta[688], theta[687], theta[686], theta[685], theta[684], theta[683], theta[682], theta[681], theta[680], theta[679], theta[678], theta[677], theta[676], theta[675], theta[674], theta[673], theta[672], theta[671], theta[670], theta[669], theta[668], theta[667], theta[666], theta[665], theta[664], theta[663], theta[662], theta[661], theta[660], theta[659], theta[658], theta[657], theta[656], theta[655], theta[654], theta[653], theta[652], theta[651], theta[650], theta[649], theta[648], theta[647], theta[646], theta[645], theta[644], theta[643], theta[642], theta[641], theta[640], theta[639], theta[638], theta[637], theta[636], theta[635], theta[634], theta[633], theta[632], theta[631], theta[630], theta[629], theta[628], theta[627], theta[626], theta[625], theta[624], theta[623], theta[622], theta[621], theta[620], theta[619], theta[618], theta[617], theta[616], theta[615], theta[614], theta[613], theta[612], theta[611], theta[610], theta[609], theta[608], theta[607], theta[606], theta[605], theta[604], theta[603], theta[602], theta[601], theta[600], theta[599], theta[598], theta[597], theta[596], theta[595], theta[594], theta[593], theta[592], theta[591], theta[590], theta[589], theta[588], theta[587], theta[586], theta[585], theta[584], theta[583], theta[582], theta[581], theta[580], theta[579], theta[578], theta[577], theta[576], theta[575], theta[574], theta[573], theta[572], theta[571], theta[570], theta[569], theta[568], theta[567], theta[566], theta[565], theta[564], theta[563], theta[562], theta[561], theta[560], theta[559], theta[558], theta[557], theta[556], theta[555], theta[554], theta[553], theta[552], theta[551], theta[550], theta[549], theta[548], theta[547], theta[546], theta[545], theta[544], theta[543], theta[542], theta[541], theta[540], theta[539], theta[538], theta[537], theta[536], theta[535], theta[534], theta[533], theta[532], theta[531], theta[530], theta[529], theta[528], theta[527], theta[526], theta[525], theta[524], theta[523], theta[522], theta[521], theta[520], theta[519], theta[518], theta[517], theta[516], theta[515], theta[514], theta[513], theta[512], theta[511], theta[510], theta[509], theta[508], theta[507], theta[506], theta[505], theta[504], theta[503], theta[502], theta[501], theta[500], theta[499], theta[498], theta[497], theta[496], theta[495], theta[494], theta[493], theta[492], theta[491], theta[490], theta[489], theta[488], theta[487], theta[486], theta[485], theta[484], theta[483], theta[482], theta[481], theta[480], theta[479], theta[478], theta[477], theta[476], theta[475], theta[474], theta[473], theta[472], theta[471], theta[470], theta[469], theta[468], theta[467], theta[466], theta[465], theta[464], theta[463], theta[462], theta[461], theta[460], theta[459], theta[458], theta[457], theta[456], theta[455], theta[454], theta[453], theta[452], theta[451], theta[450], theta[449], theta[448], theta[447], theta[446], theta[445], theta[444], theta[443], theta[442], theta[441], theta[440], theta[439], theta[438], theta[437], theta[436], theta[435], theta[434], theta[433], theta[432], theta[431], theta[430], theta[429], theta[428], theta[427], theta[426], theta[425], theta[424], theta[423], theta[422], theta[421], theta[420], theta[419], theta[418], theta[417], theta[416], theta[415], theta[414], theta[413], theta[412], theta[411], theta[410], theta[409], theta[408], theta[407], theta[406], theta[405], theta[404], theta[403], theta[402], theta[401], theta[400], theta[399], theta[398], theta[397], theta[396], theta[395], theta[394], theta[393], theta[392], theta[391], theta[390], theta[389], theta[388], theta[387], theta[386], theta[385], theta[384], theta[383], theta[382], theta[381], theta[380], theta[379], theta[378], theta[377], theta[376], theta[375], theta[374], theta[373], theta[372], theta[371], theta[370], theta[369], theta[368], theta[367], theta[366], theta[365], theta[364], theta[363], theta[362], theta[361], theta[360], theta[359], theta[358], theta[357], theta[356], theta[355], theta[354], theta[353], theta[352], theta[351], theta[350], theta[349], theta[348], theta[347], theta[346], theta[345], theta[344], theta[343], theta[342], theta[341], theta[340], theta[339], theta[338], theta[337], theta[336], theta[335], theta[334], theta[333], theta[332], theta[331], theta[330], theta[329], theta[328], theta[327], theta[326], theta[325], theta[324], theta[323], theta[322], theta[321], theta[320], theta[319], theta[318], theta[317], theta[316], theta[315], theta[314], theta[313], theta[312], theta[311], theta[310], theta[309], theta[308], theta[307], theta[306], theta[305], theta[304], theta[303], theta[302], theta[301], theta[300], theta[299], theta[298], theta[297], theta[296], theta[295], theta[294], theta[293], theta[292], theta[291], theta[290], theta[289], theta[288], theta[287], theta[286], theta[285], theta[284], theta[283], theta[282], theta[281], theta[280], theta[279], theta[278], theta[277], theta[276], theta[275], theta[274], theta[273], theta[272], theta[271], theta[270], theta[269], theta[268], theta[267], theta[266], theta[265], theta[264], theta[263], theta[262], theta[261], theta[260], theta[259], theta[258], theta[257], theta[256], theta[255], theta[254], theta[253], theta[252], theta[251], theta[250], theta[249], theta[248], theta[247], theta[246], theta[245], theta[244], theta[243], theta[242], theta[241], theta[240], theta[239], theta[238], theta[237], theta[236], theta[235], theta[234], theta[233], theta[232], theta[231], theta[230], theta[229], theta[228], theta[227], theta[226], theta[225], theta[224], theta[223], theta[222], theta[221], theta[220], theta[219], theta[218], theta[217], theta[216], theta[215], theta[214], theta[213], theta[212], theta[211], theta[210], theta[209], theta[208], theta[207], theta[206], theta[205], theta[204], theta[203], theta[202], theta[201], theta[200], theta[199], theta[198], theta[197], theta[196], theta[195], theta[194], theta[193], theta[192], theta[191], theta[190], theta[189], theta[188], theta[187], theta[186], theta[185], theta[184], theta[183], theta[182], theta[181], theta[180], theta[179], theta[178], theta[177], theta[176], theta[175], theta[174], theta[173], theta[172], theta[171], theta[170], theta[169], theta[168], theta[167], theta[166], theta[165], theta[164], theta[163], theta[162], theta[161], theta[160], theta[159], theta[158], theta[157], theta[156], theta[155], theta[154], theta[153], theta[152], theta[151], theta[150], theta[149], theta[148], theta[147], theta[146], theta[145], theta[144], theta[143], theta[142], theta[141], theta[140], theta[139], theta[138], theta[137], theta[136], theta[135], theta[134], theta[133], theta[132], theta[131], theta[130], theta[129], theta[128], theta[127], theta[126], theta[125], theta[124], theta[123], theta[122], theta[121], theta[120], theta[119], theta[118], theta[117], theta[116], theta[115], theta[114], theta[113], theta[112], theta[111], theta[110], theta[109], theta[108], theta[107], theta[106], theta[105], theta[104], theta[103], theta[102], theta[101], theta[100], theta[99], theta[98], theta[97], theta[96], theta[95], theta[94], theta[93], theta[92], theta[91], theta[90], theta[89], theta[88], theta[87], theta[86], theta[85], theta[84], theta[83], theta[82], theta[81], theta[80], theta[79], theta[78], theta[77], theta[76], theta[75], theta[74], theta[73], theta[72], theta[71], theta[70], theta[69], theta[68], theta[67], theta[66], theta[65], theta[64], theta[63], theta[62], theta[61], theta[60], theta[59], theta[58], theta[57], theta[56], theta[55], theta[54], theta[53], theta[52], theta[51], theta[50], theta[49], theta[48], theta[47], theta[46], theta[45], theta[44], theta[43], theta[42], theta[41], theta[40], theta[39], theta[38], theta[37], theta[36], theta[35], theta[34], theta[33], theta[32], theta[31], theta[30], theta[29], theta[28], theta[27], theta[26], theta[25], theta[24], theta[23], theta[22], theta[21], theta[20], theta[19], theta[18], theta[17], theta[16], theta[15], theta[14], theta[13], theta[12], theta[11], theta[10], theta[9], theta[8], theta[7], theta[6], theta[5], theta[4], theta[3], theta[2], theta[1]

  Variable sizes and types:
    alpha: size = (5,), type = Vector{Float64}
    p: size = (1000, 5), type = Matrix{Float64}
    culm: size = (32,), type = Vector{Int64}
    N: type = Int64
    R: type = Int64
    r: size = (1000, 5), type = Matrix{Int64}
    T: type = Int64
    a: size = (5,), type = Vector{Float64}
    beta: type = Float64
    response: size = (32, 5), type = Matrix{Int64}
    theta: size = (1000,), type = Vector{Float64}

All of the classic examples ship with the package under JuliaBUGS.BUGSExamples, bundling the model definition, data, initial values, and reference results (this one is JuliaBUGS.BUGSExamples.VOLUME_1.lsat).

Sampling

To draw posterior samples, construct the model with gradient support and run the No-U-Turn sampler:

using AbstractMCMC, AdvancedHMC, ADTypes, Mooncake, FlexiChains
using LogDensityProblems

model = lsat(data; adtype=AutoMooncake(; config=nothing))

n_samples, n_adapts = 2000, 1000
D = LogDensityProblems.dimension(model)
chain = AbstractMCMC.sample(
    model, NUTS(0.8), n_samples;
    chain_type=VNChain, n_adapts=n_adapts,
    init_params=rand(D), discard_initial=n_adapts,
)
summarystats(chain)

BUGS-style initial values for this example are available as JuliaBUGS.BUGSExamples.VOLUME_1.lsat.inits and can be applied with initialize!(model, inits).

Results

Reference posterior summaries are not bundled with the package for this example (the reference_results field of JuliaBUGS.BUGSExamples.VOLUME_1.lsat is empty). The published results — posterior summaries for the centred item difficulties a[1]–a[5] and the scale parameter beta — are shown in the OpenBUGS write-up of this example. A correctly converged chain's summarystats output should match those published values up to Monte Carlo error.

See also: gallery overview and the getting-started tutorial.