Add 5PL model (#14)

* add 5pl irf  function

* add 5PL

* add 5PL to docs

* add precompilation for 5PL

* rename has_asymmetry -> has_stiffness

* fix doctests

* add signatures to precompile

* check response types

.gitignore

@@ -24,3 +24,5 @@ docs/site/
 Manifest.toml
 
 .vscode
+
+dev/

Project.toml

@@ -5,7 +5,9 @@ version = "0.1.0"
 
 [deps]
 AbstractItemResponseModels = "0ab3451c-659c-47cd-a7a9-a2d579e209dd"
+DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
@@ -15,9 +17,9 @@ SimpleUnPack = "ce78b400-467f-4804-87d8-8f486da07d0a"
 AbstractItemResponseModels = "0.2"
 DocStringExtensions = "0.9"
 LogExpFunctions = "0.3"
+PrecompileTools = "1"
 Reexport = "1"
 SimpleUnPack = "1"
-PrecompileTools = "1"
 julia = "1.8"
 
 [extras]

docs/src/api.md

@@ -11,6 +11,7 @@ OneParameterLogisticPlusGuessingModel
 TwoParameterLogisticModel
 ThreeParameterLogisticModel
 FourParameterLogisticModel
+FiveParameterLogisticModel
 PartialCreditModel
 GeneralizedPartialCreditModel
 RatingScaleModel

src/ItemResponseFunctions.jl

@@ -1,17 +1,23 @@
 module ItemResponseFunctions
 
+using AbstractItemResponseModels: Dichotomous, Nominal, checkresponsetype
 using DocStringExtensions: SIGNATURES, TYPEDEF, METHODLIST
 using LogExpFunctions: logistic, cumsum!, softmax!
 using Reexport: @reexport
 using SimpleUnPack: @unpack
 
+using DifferentiationInterface
+import ForwardDiff
+
 # AbstractItemResponseModels interface extensions
 @reexport import AbstractItemResponseModels:
     ItemResponseModel, irf, iif, expected_score, information
 
-import AbstractItemResponseModels: response_type, Dichotomous
+import AbstractItemResponseModels: response_type
 
 export DichotomousItemResponseModel,
+    FivePL,
+    FiveParameterLogisticModel,
     FourPL,
     FourParameterLogisticModel,
     GPCM,
@@ -42,6 +48,6 @@ include("expected_score.jl")
 include("information.jl")
 include("scoring_functions.jl")
 
-include("precompile.jl")
+# include("precompile.jl")
 
 end

src/iif.jl

@@ -23,39 +23,39 @@ julia> iif(TwoPL, 0.0, (a = 1.3, b = 0.2))
 ### 3 Parameter Logistic Model
 ```jldoctest
 julia> iif(ThreePL, 0.0, (a = 1.5, b = 0.5, c = 0.15))
-0.3162871805861735
+0.3162871805861734
 
 ```
 ### 4 Parameter Logistic Model
 ```jldoctest
 julia> iif(FourPL, 0.0, (a = 2.1, b = -1.5, c = 0.15, d = 0.9))
-0.033871578233130605
+0.03387157823313065
 ```
 
 """
-function iif(M::Type{FourPL}, theta, beta::NamedTuple, y = 1)
-    @unpack a, b, c, d = beta
-    prob = irf(M, theta, beta, y)
-
-    num = a^2 * (prob - c)^2 * (d - prob)^2
-
-    # return early to avoid NaN cases when 0/0
-    num == 0 && return num
-
-    denum = (d - c)^2 * prob * (1 - prob)
-    info = num / denum
-
+function iif(M::Type{<:DichotomousItemResponseModel}, theta, beta::NamedTuple)
+    info = zero(theta)
+    for y in 0:1
+        info += iif(M, theta, beta, y)
+    end
     return info
 end
 
-function iif(M::Type{<:DichotomousItemResponseModel}, theta, beta, y = 1)
-    pars = merge_pars(M, beta)
-    return iif(FourPL, theta, pars, y)
+function iif(M::Type{<:DichotomousItemResponseModel}, theta, beta::NamedTuple, y)
+    checkresponsetype(response_type(M), y)
+    return _iif(M, theta, beta, y)
 end
 
-function iif(M::Type{OnePL}, theta::Real, beta::Real, y = 1)
-    prob = irf(M, theta, beta, y)
-    return prob * (1 - prob)
+iif(M::Type{OnePL}, theta, beta::Real, y) = iif(M, theta, (; b = beta), y)
+iif(M::Type{OnePL}, theta, beta::Real) = iif(M, theta, (; b = beta))
+
+function _iif(M::Type{<:DichotomousItemResponseModel}, theta, beta, y)
+    adtype = AutoForwardDiff()
+    f = x -> irf(M, x, beta, y)
+    prob, deriv = value_and_derivative(f, adtype, theta)
+    iszero(prob) && return 0.0  # early return to avoid NaNs
+    deriv2 = second_derivative(f, adtype, theta)
+    return deriv^2 / prob - deriv2
 end
 
 function iif(M::Type{GPCM}, theta, beta; scoring_function::F = identity) where {F}
@@ -92,6 +92,7 @@ function iif(
     y;
     scoring_function::F = identity,
 ) where {F}
+    checkresponsetype(response_type(M), y)
     prob = irf(M, theta, beta, y)
     return prob * iif(M, theta, beta; scoring_function)
 end

src/information.jl

@@ -33,7 +33,7 @@ julia> information(TwoPL, 0.0, betas)
 julia> betas = fill((; a = 1.5, b = 0.5, c = 0.2), 4);
 
 julia> information(ThreePL, 0.0, betas)
-1.102180659985265
+1.1021806599852657
 ```
 
 ### 4 Parameter Logistic Model
@@ -52,7 +52,7 @@ function information(
     info = zero(T)
 
     for beta in betas
-        info += iif(M, theta, beta, 1)
+        info += iif(M, theta, beta)
     end
 
     return info

src/irf.jl

@@ -89,35 +89,41 @@ julia> irf(RSM, 0.0, beta, 3)
 ```
 
 """
-function irf(M::Type{FourPL}, theta::Real, beta::NamedTuple, y = 1)
-    @unpack a, b, c, d = beta
-    prob = c + (d - c) * logistic(a * (theta - b))
-    return ifelse(y == 1, prob, 1 - prob)
-end
-
 function irf(M::Type{<:DichotomousItemResponseModel}, theta, beta, y = 1)
+    checkresponsetype(response_type(M), y)
     pars = merge_pars(M, beta)
-    return irf(FourPL, theta, pars, y)
+    return _irf(M, theta, pars, y)
 end
 
 function irf(M::Type{OnePL}, theta::Real, beta::Real, y = 1)
+    checkresponsetype(response_type(M), y)
     prob = logistic(theta - beta)
     return ifelse(y == 1, prob, 1 - prob)
 end
 
-function irf(M::Type{GPCM}, theta, beta)
+function _irf(M::Type{<:DichotomousItemResponseModel}, theta, beta::NamedTuple, y)
+    @unpack a, b, c, d, e = beta
+    prob = c + (d - c) * logistic(a * (theta - b))^e
+    return ifelse(y == 1, prob, 1 - prob)
+end
+
+# polytomous models
+function irf(M::Type{<:PolytomousItemResponseModel}, theta, beta)
+    pars = merge_pars(M, beta)
+    return _irf(GPCM, theta, pars)
+end
+
+function _irf(M::Type{GPCM}, theta, beta)
     @unpack t = beta
     probs = similar(t, length(t) + 1)
     return irf!(M, probs, theta, beta)
 end
 
-function irf(M::Type{<:PolytomousItemResponseModel}, theta, beta)
-    pars = has_discrimination(M) ? beta : merge(beta, (; a = 1.0))
-    return irf(GPCM, theta, pars)
+function irf(M::Type{<:PolytomousItemResponseModel}, theta, beta, y)
+    checkresponsetype(response_type(M), y)
+    return irf(M, theta, beta)[y]
 end
 
-irf(M::Type{<:PolytomousItemResponseModel}, theta, beta, y) = irf(M, theta, beta)[y]
-
 """
     $(SIGNATURES)
 
@@ -146,15 +152,15 @@ julia> probs
  0.13454527815807202
 ```
 """
-function irf!(M::Type{GPCM}, probs, theta, beta)
+function irf!(M::Type{<:PolytomousItemResponseModel}, probs, theta, beta)
+    return _irf!(GPCM, probs, theta, beta)
+end
+
+function _irf!(M::Type{GPCM}, probs, theta, beta)
     @unpack a, b, t = beta
     probs[1] = 0.0
     @. probs[2:end] = a * (theta - b + t)
     cumsum!(probs, probs)
     softmax!(probs, probs)
     return probs
 end
-
-function irf!(M::Type{<:PolytomousItemResponseModel}, probs, theta, beta)
-    return irf!(GPCM, probs, theta, beta)
-end

src/model_types.jl

@@ -14,6 +14,12 @@ An abstract type representing an item response model with polytomous responses.
 """
 abstract type PolytomousItemResponseModel <: ItemResponseModel end
 
+response_type(::Type{<:PolytomousItemResponseModel}) = Nominal
+
+has_stiffness(::Type{<:PolytomousItemResponseModel}) = false
+has_lower_asymptote(::Type{<:PolytomousItemResponseModel}) = false
+has_upper_asymptote(::Type{<:PolytomousItemResponseModel}) = false
+
 """
     $(TYPEDEF)
 
@@ -35,6 +41,7 @@ const OnePL = OneParameterLogisticModel
 has_discrimination(::Type{OnePL}) = false
 has_lower_asymptote(::Type{OnePL}) = false
 has_upper_asymptote(::Type{OnePL}) = false
+has_stiffness(::Type{OnePL}) = false
 
 """
     $(TYPEDEF)
@@ -57,6 +64,7 @@ const OnePLG = OneParameterLogisticPlusGuessingModel
 has_discrimination(::Type{OnePLG}) = false
 has_lower_asymptote(::Type{OnePLG}) = true
 has_upper_asymptote(::Type{OnePLG}) = false
+has_stiffness(::Type{OnePLG}) = false
 
 """
     $(TYPEDEF)
@@ -79,6 +87,7 @@ const TwoPL = TwoParameterLogisticModel
 has_discrimination(::Type{TwoPL}) = true
 has_lower_asymptote(::Type{TwoPL}) = false
 has_upper_asymptote(::Type{TwoPL}) = false
+has_stiffness(::Type{TwoPL}) = false
 
 """
     $(TYPEDEF)
@@ -102,6 +111,7 @@ const ThreePL = ThreeParameterLogisticModel
 has_discrimination(::Type{ThreePL}) = true
 has_lower_asymptote(::Type{ThreePL}) = true
 has_upper_asymptote(::Type{ThreePL}) = false
+has_stiffness(::Type{ThreePL}) = false
 
 """
     $(TYPEDEF)
@@ -126,6 +136,31 @@ const FourPL = FourParameterLogisticModel
 has_discrimination(::Type{FourPL}) = true
 has_lower_asymptote(::Type{FourPL}) = true
 has_upper_asymptote(::Type{FourPL}) = true
+has_stiffness(::Type{FourPL}) = false
+
+"""
+    $(TYPEDEF)
+
+An abstract representation of a 5 Parameter Logistic Model with an item response function
+given by
+
+``P(Y_{ij}=1|\\theta_i,\\boldsymbol{\\beta}_j) = c_j + (d_j - c_j)\\cdot\\mathrm{logistic}(a_j(\\theta_i - b_j)^e_j)``
+
+The item parameters `beta` must be a destructurable object with the following fields:
+
+- `a`: the item discrimination
+- `b`: the item difficulty (location)
+- `c`: the lower asymptote
+- `d`: the upper asymptote
+- `e`: the item asymmetry
+"""
+abstract type FiveParameterLogisticModel <: DichotomousItemResponseModel end
+const FivePL = FiveParameterLogisticModel
+
+has_discrimination(::Type{FivePL}) = true
+has_lower_asymptote(::Type{FivePL}) = true
+has_upper_asymptote(::Type{FivePL}) = true
+has_stiffness(::Type{FivePL}) = true
 
 """
     $(TYPEDEF)
@@ -149,8 +184,6 @@ abstract type GeneralizedPartialCreditModel <: PolytomousItemResponseModel end
 const GPCM = GeneralizedPartialCreditModel
 
 has_discrimination(::Type{GPCM}) = true
-has_lower_asymptote(::Type{GPCM}) = false
-has_upper_asymptote(::Type{GPCM}) = false
 
 """
     $(TYPEDEF)
@@ -173,8 +206,6 @@ abstract type PartialCreditModel <: PolytomousItemResponseModel end
 const PCM = PartialCreditModel
 
 has_discrimination(::Type{PCM}) = false
-has_lower_asymptote(::Type{PCM}) = false
-has_upper_asymptote(::Type{PCM}) = false
 
 """
     $(TYPEDEF)
@@ -197,8 +228,6 @@ abstract type RatingScaleModel <: PolytomousItemResponseModel end
 const RSM = RatingScaleModel
 
 has_discrimination(::Type{RSM}) = false
-has_lower_asymptote(::Type{RSM}) = false
-has_upper_asymptote(::Type{RSM}) = false
 
 """
     $(TYPEDEF)
@@ -222,5 +251,3 @@ abstract type GeneralizedRatingScaleModel <: PolytomousItemResponseModel end
 const GRSM = GeneralizedRatingScaleModel
 
 has_discrimination(::Type{GRSM}) = true
-has_lower_asymptote(::Type{GRSM}) = false
-has_upper_asymptote(::Type{GRSM}) = false

src/precompile.jl

@@ -1,14 +1,16 @@
 using PrecompileTools: @setup_workload, @compile_workload
 
 @setup_workload begin
-    models = [OnePL, OnePLG, TwoPL, ThreePL, FourPL, PCM, GPCM, RSM, GRSM]
-    beta = (a = 1.0, b = 0.0, c = 0.0, d = 1.0, t = zeros(3))
+    models = [OnePL, OnePLG, TwoPL, ThreePL, FourPL, FivePL, PCM, GPCM, RSM, GRSM]
+    beta = (a = 1.0, b = 0.0, c = 0.0, d = 1.0, e = 1.0, t = zeros(3))
     betas = fill(beta, 3)
 
     @compile_workload begin
         for model in models
             irf(model, 0.0, beta, 1)
+            irf(model, 0.0, beta)
             iif(model, 0.0, beta, 1)
+            iif(model, 0.0, beta)
             expected_score(model, 0.0, betas)
             information(model, 0.0, betas)
         end

src/utils.jl

@@ -18,5 +18,9 @@ function merge_pars(M::Type{<:ItemResponseModel}, beta)
         pars = merge(pars, (; d = 1.0))
     end
 
+    if !has_stiffness(M)
+        pars = merge(pars, (; e = 1.0))
+    end
+
     return pars
 end

test/models/five_parameter_logistic.jl

@@ -0,0 +1,8 @@
+@testset "FiveParameterLogisticModel" begin
+    T = FivePL
+
+    @test has_discrimination(T) == true
+    @test has_lower_asymptote(T) == true
+    @test has_upper_asymptote(T) == true
+    @test has_stiffness(T) == true
+end