[Tests] Improve tests

src/EdererI.jl

@@ -15,15 +15,15 @@ To call this function:
 """
 const EdererI = NPNSEstimator{EdererIMethod}
 
-function Λ!(::Type{EdererIMethod}, num_excess, den_excess, num_pop, den_pop, num_variance, T, Δ, age, year, rate_preds, ratetable, grid)
+function Λ!(::Type{EdererIMethod}, num_excess, den_excess, num_pop, den_pop, num_variance, T, Δ, age, year, rate_preds, ratetable, grid, ∂t)
     Tmax= Int(maximum(T))
     for i in eachindex(age)
         Tᵢ = searchsortedlast(grid, T[i])
         Λₚ = 0.0
         rtᵢ = ratetable[rate_preds[i,:]...]
         for j in 1:Tmax
             λₚ          = daily_hazard(rtᵢ, age[i] + grid[j], year[i] + grid[j])
-            ∂Λₚ         = λₚ * (grid[j+1]-grid[j]) 
+            ∂Λₚ         = λₚ * ∂t[j]
             Λₚ         += ∂Λₚ
             Sₚ          = exp(-Λₚ)
             num_pop[j] += (Sₚ * ∂Λₚ)

src/EdererII.jl

@@ -15,13 +15,13 @@ To call this function:
 """
 const EdererII = NPNSEstimator{EdererIIMethod}
 
-function Λ!(::Type{EdererIIMethod}, num_excess, den_excess, num_pop, den_pop, num_variance, T, Δ, age, year, rate_preds, ratetable, grid)
+function Λ!(::Type{EdererIIMethod}, num_excess, den_excess, num_pop, den_pop, num_variance, T, Δ, age, year, rate_preds, ratetable, grid, ∂t)
     for i in eachindex(age)
         Tᵢ = searchsortedlast(grid, T[i])
         rtᵢ = ratetable[rate_preds[i,:]...]
         for j in 1:Tᵢ
             λₚ          = daily_hazard(rtᵢ, age[i] + grid[j], year[i] + grid[j])
-            ∂Λₚ         = λₚ * (grid[j+1]-grid[j])
+            ∂Λₚ         = λₚ * ∂t[j]
             den_excess[j] += 1
             den_pop[j] += 1
             num_pop[j] += ∂Λₚ

src/GraffeoTest.jl

@@ -86,6 +86,7 @@ struct GraffeoTest
         num_variance = zero(grid)
         den_pop      = zero(grid)
         den_excess   = zero(grid)
+        ∂t = [diff(grid)...,1.0]
 
         # Compute Pohar Perme numerator and denominators on each strata&group (s,g)
         for s in eachindex(stratas)
@@ -97,7 +98,7 @@ struct GraffeoTest
                 num_variance .= 0
                 den_pop      .= 0
                 den_excess   .= 0
-                Λ!(PoharPermeMethod, num_excess, den_excess, num_pop, den_pop, num_variance, T[idx], Δ[idx], age[idx], year[idx], rate_preds[idx,:], ratetable, grid)
+                Λ!(PoharPermeMethod, num_excess, den_excess, num_pop, den_pop, num_variance, T[idx], Δ[idx], age[idx], year[idx], rate_preds[idx,:], ratetable, grid, ∂t)
                 ∂N[s, g, :] .= num_excess.- num_pop
                 ∂V[s, g, :] .= num_variance
                 D[s, g, :]  .= den_excess
@@ -149,7 +150,7 @@ function StatsBase.fit(::Type{E}, formula::FormulaTerm, df::DataFrame, rt::RateT
     are_strata = [t <: FunctionTerm{typeof(Strata)} for t in types]
 
     strata = groupindices(groupby(df,terms[are_strata]))
-    group  = groupindices(groupby(df,terms[(!).(are_strata)]))
+    group  = groupindices(groupby(df,terms))
 
     resp = modelcols(apply_schema(formula,schema(df)).lhs,df)
     rate_predictors = _get_rate_predictors(rt,df)

src/Hakulinen.jl

@@ -9,7 +9,7 @@ To call this function:
 """
 const Hakulinen = NPNSEstimator{HakulinenMethod}
 
-function Λ!(::Type{HakulinenMethod}, num_excess, den_excess, num_pop, den_pop, num_variance, T, Δ, age, year, rate_preds, ratetable, grid)
+function Λ!(::Type{HakulinenMethod}, num_excess, den_excess, num_pop, den_pop, num_variance, T, Δ, age, year, rate_preds, ratetable, grid, ∂t)
     Date_max = maximum(T .+ year)
     for i in eachindex(age)
         Tᵢ = searchsortedlast(grid, T[i])
@@ -24,7 +24,7 @@ function Λ!(::Type{HakulinenMethod}, num_excess, den_excess, num_pop, den_pop,
         rtᵢ = ratetable[rate_preds[i,:]...]
         for j in 1:T2
             λₚ          = daily_hazard(rtᵢ, age[i] + grid[j], year[i] + grid[j])
-            ∂Λₚ         = λₚ * (grid[j+1]-grid[j]) 
+            ∂Λₚ         = λₚ * ∂t[j]
             Λₚ         += ∂Λₚ
             Sₚ          = exp(-Λₚ)
             num_pop[j] += (Sₚ * ∂Λₚ)

src/NPNSEstimator.jl

@@ -17,7 +17,7 @@ struct NPNSEstimator{Method} <: StatisticalModel
 end
 
 function mk_grid(times,prec)
-    M = maximum(times)+1
+    M = maximum(times)
     return unique(sort([(1:prec:M)..., times..., M]))
 end
 function Λ(::Type{M}, T, Δ, age, year, rate_preds, ratetable, grid) where M<:NPNSMethod
@@ -26,7 +26,8 @@ function Λ(::Type{M}, T, Δ, age, year, rate_preds, ratetable, grid) where M<:N
     num_variance = zero(grid)
     den_pop      = zero(grid)
     den_excess   = zero(grid)
-    Λ!(M, num_excess, den_excess, num_pop, den_pop, num_variance, T, Δ, age, year, rate_preds, ratetable, grid)
+    ∂t = [diff(grid)...,1.0]
+    Λ!(M, num_excess, den_excess, num_pop, den_pop, num_variance, T, Δ, age, year, rate_preds, ratetable, grid, ∂t)
     return num_excess ./ den_excess, num_pop ./ den_pop, num_variance ./ (den_excess.^2)
 end
 

src/PoharPerme.jl

@@ -18,15 +18,15 @@ References:
 """
 const PoharPerme = NPNSEstimator{PoharPermeMethod}
 
-function Λ!(::Type{PoharPermeMethod}, num_excess, den_excess, num_pop, den_pop, num_variance, T, Δ, age, year, rate_preds, ratetable, grid)
+function Λ!(::Type{PoharPermeMethod}, num_excess, den_excess, num_pop, den_pop, num_variance, T, Δ, age, year, rate_preds, ratetable, grid, ∂t)
     for i in eachindex(age)
         Tᵢ = searchsortedlast(grid, T[i])
         Λₚ = 0.0
         wₚ = 1.0
         rtᵢ = ratetable[rate_preds[i,:]...]
         for j in 1:Tᵢ
             λₚ          = daily_hazard(rtᵢ, age[i] + grid[j], year[i] + grid[j])
-            ∂Λₚ         = λₚ * (grid[j+1]-grid[j]) # λₚ * ∂t 
+            ∂Λₚ         = λₚ * ∂t[j]
             Λₚ         += ∂Λₚ
             wₚ          = exp(Λₚ)
             num_pop[j] += ∂Λₚ * wₚ