From 97c4de6af9e4d69b5b8e2abb351ed2f0acf22d55 Mon Sep 17 00:00:00 2001
From: PharmCat <13901158+PharmCat@users.noreply.github.com>
Date: Sat, 5 Nov 2022 01:49:51 +0300
Subject: [PATCH 01/14] Type III table

---
 Project.toml |   2 +
 src/ftest.jl | 115 +++++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 117 insertions(+)

diff --git a/Project.toml b/Project.toml
index 83b8d844..0bc32a67 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,6 +6,7 @@ version = "1.8.1"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
@@ -20,6 +21,7 @@ CSV = "0.7, 0.8"
 CategoricalArrays = "0.8, 0.9, 0.10"
 DataFrames = "0.22, 1"
 Distributions = "0.16, 0.17, 0.18, 0.19, 0.20, 0.21, 0.22, 0.23, 0.24, 0.25"
+PrettyTables = "1, 2"
 RDatasets = "0.5, 0.6, 0.7"
 Reexport = "0.1, 0.2, 1.0"
 SpecialFunctions = "0.6, 0.7, 0.8, 0.9, 0.10, 1, 2.0"
diff --git a/src/ftest.jl b/src/ftest.jl
index 00a19437..8ca63975 100644
--- a/src/ftest.jl
+++ b/src/ftest.jl
@@ -229,3 +229,118 @@ function show(io::IO, ftr::FTestResult{N}) where N
     end
     print(io, '─'^totwidth)
 end
+
+
+##############################################
+# Group effects table
+# Baset on F-statistics
+# The s×p full row rank matrix. The rows are estimable functions. s≥1
+
+using PrettyTables
+struct GroupEffectsTable
+    name
+    df
+    f
+    p
+end
+
+"""
+θ + A * B * A'
+
+Change θ (only upper triangle). B is symmetric.
+"""
+function mulαβαtinc!(θ::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix)
+    axb  = axes(B, 1)
+    sa   = size(A, 1)
+    for m ∈ 1:sa
+        for n ∈ m:sa
+            for j ∈ axb
+                @inbounds for i ∈ axb
+                    θ[m, n] +=  A[m, i] * B[i, j] * A[n, j]
+                end
+            end
+        end
+    end
+    θ
+end
+"""
+    lcontrast(obj, i::Int)
+
+L-contrast matrix for `i` fixed effect.
+"""
+function lcontrast(obj, i::Int)
+    n = obj.mf.schema.schema.count
+    if i > n || n < 1 error("Factor number out of range 1-$(n)") end
+    inds = findall(x -> x==i, obj.mm.assign)
+    if typeof(obj.mf.f.rhs.terms[i]) <: CategoricalTerm
+        mxc   = zeros(size(obj.mf.f.rhs.terms[i].contrasts.matrix, 1), size(obj.mm.m, 2))
+        mxcv  = view(mxc, :, inds)
+        mxcv .= obj.mf.f.rhs.terms[i].contrasts.matrix
+        mx    = zeros(size(obj.mf.f.rhs.terms[i].contrasts.matrix, 1) - 1, size(obj.mm.m, 2))
+        for i = 2:size(obj.mf.f.rhs.terms[i].contrasts.matrix, 1)
+            mx[i-1, :] .= mxc[i, :] - mxc[1, :]
+        end
+    else
+        mx = zeros(length(inds), size(obj.mm.m, 2))
+        for i = 1:length(inds)
+            mx[i, inds[i]] = 1
+        end
+    end
+    mx
+end
+
+"""
+"""
+function typeiii(obj)
+    V           = vcov(obj) 
+    replace!(V, NaN => 0)
+    B           = coef(obj)   
+    c           = obj.mf.schema.schema.count
+    d           = Vector{Int}(undef, 0)
+    fac         = Vector{String}(undef, c)
+    F           = Vector{Float64}(undef,c)
+    df          = Vector{Float64}(undef, c)
+    ndf         = Vector{Float64}(undef, c)
+    pval        = Vector{Float64}(undef, c)
+    for i = 1:c
+        if typeof(obj.mf.f.rhs.terms[i]) <: InterceptTerm{true}
+            fac[i] = "(Intercept)"
+        elseif typeof(obj.mf.f.rhs.terms[i]) <: InterceptTerm{false}
+            push!(d, i)
+            fac[i] = ""
+            continue
+        else
+            fac[i] = string(obj.mf.f.rhs.terms[i].sym)
+        end
+        L       = lcontrast(obj, i)
+        for c = 1:length(B)
+            if isnan(B[i]) || iszero(B[i])
+                L[i, :] .= 0
+            end
+        end
+
+        θ = zeros(size(L, 1), size(L, 1))
+        mulαβαtinc!(θ, L, V)
+
+
+        F[i]    = (B' * L' * pinv(Symmetric(θ)) * L * B)/rank(L)
+
+        df[i]  = rank(L)
+
+        pval[i] = ccdf(FDist(df[i], dof_residual(obj)), F[i])
+  
+    end
+    if length(d) > 0
+        deleteat!(fac, d)
+        deleteat!(F, d)
+        deleteat!(df, d)
+        deleteat!(pval, d)
+    end
+    GroupEffectsTable(fac, df, F, pval)
+
+end
+
+function Base.show(io::IO, obj::GroupEffectsTable)
+    mx = hcat(obj.name, obj.df, obj.f, obj.p)
+    PrettyTables.pretty_table(io, mx; tf = PrettyTables.tf_compact, header = ["Name", "DF" ,"F" ,"Pval"])
+end
\ No newline at end of file

From 4693bb18dfefe93c4ff6ce41a6fac813bce23110 Mon Sep 17 00:00:00 2001
From: PharmCat <13901158+PharmCat@users.noreply.github.com>
Date: Sat, 5 Nov 2022 04:52:19 +0300
Subject: [PATCH 02/14] update

---
 src/ftest.jl | 100 ++++++++++++++++++++++++++++++++++++++++++---------
 1 file changed, 84 insertions(+), 16 deletions(-)

diff --git a/src/ftest.jl b/src/ftest.jl
index 8ca63975..53e7ac5b 100644
--- a/src/ftest.jl
+++ b/src/ftest.jl
@@ -232,16 +232,15 @@ end
 
 
 ##############################################
-# Group effects table
+# Tests of Between-Subjects Effects
 # Baset on F-statistics
-# The s×p full row rank matrix. The rows are estimable functions. s≥1
+# L: The s×p full row rank matrix. The rows are estimable functions. s≥1 where p number of coefs
 
-using PrettyTables
 struct GroupEffectsTable
-    name
-    df
-    f
-    p
+    name::Vector{String}
+    df::Vector{Float64}
+    f::Vector{Float64}
+    p::Vector{Float64}
 end
 
 """
@@ -263,6 +262,65 @@ function mulαβαtinc!(θ::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix
     end
     θ
 end
+"""
+a' * B * a
+
+"""
+function mulαtβα(a::AbstractVector, B::AbstractMatrix{T}) where T
+    if length(a) != size(B, 2)::Int  || size(B, 2)::Int  != size(B, 1)::Int  error("Dimention error") end
+    axbm  = axes(B, 1)
+    axbn  = axes(B, 2)
+    c = zero(T)
+    for i ∈ axbm
+        ct = zero(T)
+        for j ∈ axbn
+            @inbounds  ct += B[i, j] * a[j]
+        end
+        @inbounds c += ct * a[i]
+    end
+    c
+end
+
+# See SPSS (GLM/UNIANOVA) and SAS (PROC GLM) documentation 
+# https://www.ibm.com/docs/en/spss-statistics/29.0.0?topic=effects-tests-between-subjects
+# L is a s×p matrix corresponding to plan-matrix of Factor
+# p - number of columns - coefs number
+# s - number of levels for this factor in the model
+# For Example
+# If you have model matrix with Intercept and two factors A and B with 3 and 4 levels
+# with Dummy coding you will have:
+# 
+# I A2 A3 B2 B3 B4
+# 1 1  0  1  0  0
+# 1 1  0  1  0  0
+# 1 1  0  1  0  0
+# 1 1  0  1  0  0
+# 1 0  1  1  0  0
+# 1 0  1  0  1  0
+# 1 0  1  0  1  0
+# 1 0  1  0  1  0
+# 1 0  1  0  1  0
+# 1 0  0  0  0  1
+# 1 0  0  0  0  1
+# 1 0  0  0  0  1
+# 1 0  0  0  0  0
+# 1 0  0  0  0  0
+#
+# Then yoy wil have L matrix fo intercept:
+#
+# 1 0  0  0  0  0 
+#
+# For A:
+#
+# 0 1  0  0  0  0
+# 0 0  1  0  0  0
+#
+# For B:
+#
+# 0 0  0  1  0  0
+# 0 0  0  0  1  0
+# 0 0  0  0  0  1
+#
 """
     lcontrast(obj, i::Int)
 
@@ -290,6 +348,11 @@ function lcontrast(obj, i::Int)
 end
 
 """
+    typeiii(obj)
+
+Calculate F-statistics for Tests of Between-Subjects Effects. 
+Sum of squares and MS not calculated.
+
 """
 function typeiii(obj)
     V           = vcov(obj) 
@@ -300,12 +363,11 @@ function typeiii(obj)
     fac         = Vector{String}(undef, c)
     F           = Vector{Float64}(undef,c)
     df          = Vector{Float64}(undef, c)
-    ndf         = Vector{Float64}(undef, c)
     pval        = Vector{Float64}(undef, c)
     for i = 1:c
         if typeof(obj.mf.f.rhs.terms[i]) <: InterceptTerm{true}
             fac[i] = "(Intercept)"
-        elseif typeof(obj.mf.f.rhs.terms[i]) <: InterceptTerm{false}
+        elseif typeof(obj.mf.f.rhs.terms[i]) <: InterceptTerm{false} # If zero intercept
             push!(d, i)
             fac[i] = ""
             continue
@@ -318,15 +380,19 @@ function typeiii(obj)
                 L[i, :] .= 0
             end
         end
-
+        RL = rank(L) # Rank of L matrix
+        # F-statistics computed:
+        # F[i]    = (L'*B' * pinv(L * V * L') * L * B) / rank(L)
+        # As V is symmetric we can calc only upper triangle
+        # θ = L * V * L'
         θ = zeros(size(L, 1), size(L, 1))
         mulαβαtinc!(θ, L, V)
-
-
-        F[i]    = (B' * L' * pinv(Symmetric(θ)) * L * B)/rank(L)
-
-        df[i]  = rank(L)
-
+        LB = L * B
+        # Then F can be computed:
+        # F[i]    = (LB' * pinv(Symmetric(θ)) * LB)/rank(L)
+        # I think this is more efficient:
+        F[i]    = mulαtβα(LB, pinv(Symmetric(θ))) / RL
+        df[i]   = RL
         pval[i] = ccdf(FDist(df[i], dof_residual(obj)), F[i])
   
     end
@@ -340,6 +406,8 @@ function typeiii(obj)
 
 end
 
+using PrettyTables
+# use PrettyTables because it is fastest way to make table
 function Base.show(io::IO, obj::GroupEffectsTable)
     mx = hcat(obj.name, obj.df, obj.f, obj.p)
     PrettyTables.pretty_table(io, mx; tf = PrettyTables.tf_compact, header = ["Name", "DF" ,"F" ,"Pval"])

From 283e34cc85733c539b0f670e3f185e767bbc6bb2 Mon Sep 17 00:00:00 2001
From: PharmCat <13901158+PharmCat@users.noreply.github.com>
Date: Sat, 5 Nov 2022 19:11:59 +0300
Subject: [PATCH 03/14] fix rank and update

---
 src/ftest.jl | 44 ++++++++++++++++++++++++++++----------------
 1 file changed, 28 insertions(+), 16 deletions(-)

diff --git a/src/ftest.jl b/src/ftest.jl
index 53e7ac5b..3e05a0cc 100644
--- a/src/ftest.jl
+++ b/src/ftest.jl
@@ -327,26 +327,31 @@ end
 L-contrast matrix for `i` fixed effect.
 """
 function lcontrast(obj, i::Int)
-    n = obj.mf.schema.schema.count
+    n = length(obj.mf.f.rhs.terms)
     if i > n || n < 1 error("Factor number out of range 1-$(n)") end
+    p    = size(obj.mm.m, 2) # number of coefs
     inds = findall(x -> x==i, obj.mm.assign)
     if typeof(obj.mf.f.rhs.terms[i]) <: CategoricalTerm
-        mxc   = zeros(size(obj.mf.f.rhs.terms[i].contrasts.matrix, 1), size(obj.mm.m, 2))
+        mxc   = zeros(size(obj.mf.f.rhs.terms[i].contrasts.matrix, 1), p)
         mxcv  = view(mxc, :, inds)
         mxcv .= obj.mf.f.rhs.terms[i].contrasts.matrix
-        mx    = zeros(size(obj.mf.f.rhs.terms[i].contrasts.matrix, 1) - 1, size(obj.mm.m, 2))
-        for i = 2:size(obj.mf.f.rhs.terms[i].contrasts.matrix, 1)
+        mx    = zeros(size(obj.mf.f.rhs.terms[i].contrasts.matrix, 1) - 1, p)
+        for i = 2:size(obj.mf.f.rhs.terms[i].contrasts.matrix, 1) # correct for zero-intercept model
             mx[i-1, :] .= mxc[i, :] - mxc[1, :]
         end
     else
-        mx = zeros(length(inds), size(obj.mm.m, 2))
-        for i = 1:length(inds)
-            mx[i, inds[i]] = 1
+        mx = zeros(length(inds), p) # unknown correctness for zero-intercept model
+        for j = 1:length(inds)
+            mx[j, inds[j]] = 1
         end
     end
     mx
 end
 
+tname(t::AbstractTerm) = "$(t.sym)"
+tname(t::InteractionTerm) = join(tname.(t.terms), " & ")
+tname(t::InterceptTerm) = "(Intercept)"
+
 """
     typeiii(obj)
 
@@ -356,28 +361,36 @@ Sum of squares and MS not calculated.
 """
 function typeiii(obj)
     V           = vcov(obj) 
-    replace!(V, NaN => 0)
+    replace!(V, NaN => 0) # Some values can be NaN - replace it to zero
     B           = coef(obj)   
-    c           = obj.mf.schema.schema.count
+    c           = length(obj.mf.f.rhs.terms)
     d           = Vector{Int}(undef, 0)
     fac         = Vector{String}(undef, c)
     F           = Vector{Float64}(undef,c)
     df          = Vector{Float64}(undef, c)
     pval        = Vector{Float64}(undef, c)
     for i = 1:c
+        # Make L matrix
+        L       = lcontrast(obj, i)
+        #=
         if typeof(obj.mf.f.rhs.terms[i]) <: InterceptTerm{true}
-            fac[i] = "(Intercept)"
-        elseif typeof(obj.mf.f.rhs.terms[i]) <: InterceptTerm{false} # If zero intercept
+            # I think L matrix for Intercept should represent general mean
+            # For this - L element corresponding to factor should be devided on total number of levels (not coefs)
+            # But because I can't get efficient number of levels for InteractionTerm
+            # I cant't calc it correctly
+        end
+        =#
+        if typeof(obj.mf.f.rhs.terms[i]) <: InterceptTerm{false} # If zero intercept (drop)
             push!(d, i)
             fac[i] = ""
             continue
         else
-            fac[i] = string(obj.mf.f.rhs.terms[i].sym)
+            fac[i] = tname(obj.mf.f.rhs.terms[i])
         end
-        L       = lcontrast(obj, i)
+        # For case when cofs is zero (or NaN) we reduce rank of L-matrix
         for c = 1:length(B)
-            if isnan(B[i]) || iszero(B[i])
-                L[i, :] .= 0
+            if isnan(B[c]) || iszero(B[c])
+                L[:, c] .= 0
             end
         end
         RL = rank(L) # Rank of L matrix
@@ -394,7 +407,6 @@ function typeiii(obj)
         F[i]    = mulαtβα(LB, pinv(Symmetric(θ))) / RL
         df[i]   = RL
         pval[i] = ccdf(FDist(df[i], dof_residual(obj)), F[i])
-  
     end
     if length(d) > 0
         deleteat!(fac, d)

From a6c4ddce0e90fa5704d3c9fc6b0035aa6be7c6fd Mon Sep 17 00:00:00 2001
From: PharmCat <13901158+PharmCat@users.noreply.github.com>
Date: Sun, 6 Nov 2022 16:10:46 +0300
Subject: [PATCH 04/14] test

---
 data/rds1.csv    | 37 +++++++++++++++++++++++++++++++++++++
 test/runtests.jl | 27 +++++++++++++++++++++++++++
 2 files changed, 64 insertions(+)
 create mode 100644 data/rds1.csv

diff --git a/data/rds1.csv b/data/rds1.csv
new file mode 100644
index 00000000..082716be
--- /dev/null
+++ b/data/rds1.csv
@@ -0,0 +1,37 @@
+Subject	Sequence	Period	Formulation	Var
+3	TR	1	T	225.95
+1	RT	1	R	181.09
+2	RT	1	R	114.48
+4	RT	1	R	176.91
+5	TR	1	T	147.01
+6	TR	1	T	97.53
+7	RT	1	R	146.60
+8	TR	1	T	45.58
+9	RT	1	R	109.20
+10	RT	1	R	125.61
+11	TR	1	T	92.26
+12	RT	1	R	237.95
+13	TR	1	T	145.46
+14	TR	1	T	179.96
+15	TR	1	T	173.86
+16	RT	1	R	144.00
+17	RT	1	R	185.10
+18	TR	1	T	117.99
+1	RT	2	T	210.14
+2	RT	2	T	98.72
+3	TR	2	R	241.09
+4	RT	2	T	186.65
+5	TR	2	R	139.56
+6	TR	2	R	124.77
+7	RT	2	T	137.62
+8	TR	2	R	57.71
+9	RT	2	T	139.36
+10	RT	2	T	120.43
+11	TR	2	R	116.10
+12	RT	2	T	228.63
+13	TR	2	R	165.09
+14	TR	2	R	181.09
+15	TR	2	R	206.66
+16	RT	2	T	143.25
+17	RT	2	T	192.22
+18	TR	2	R	125.50
diff --git a/test/runtests.jl b/test/runtests.jl
index e89e3466..99fa0070 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1603,3 +1603,30 @@ end
     # 3. 44 / wt == y
     @test GLM.loglik_obs(Binomial(), y, μ, wt, ϕ) ≈ GLM.logpdf(Binomial(Int(wt), μ), 44)
 end
+
+berds1 = CSV.read(joinpath(glm_datadir, "rds1.csv"), DataFrame)
+berds1.Period = categorical(berds1.Period)
+berds1.Subject = categorical(berds1.Subject)
+
+@testset "Tests of Between-Subjects Effects" begin
+    # This is not BE model - no subject 
+    # Against SPSS 28
+    #=
+    GLM Var BY Sequence Period Formulation 
+  /METHOD=SSTYPE(3) 
+  /INTERCEPT=INCLUDE 
+  /PRINT PARAMETER 
+  /CRITERIA=ALPHA(.05) 
+  /DESIGN=Sequence Period Formulation.
+    =#
+    # Intercept not included in test
+    ols  = lm(@formula(Var ~ Sequence+Period+Formulation), berds1)
+    tbl = GLM.typeiii(ols)
+    @test tbl.f[2] ≈ 1.011001 atol = 1.0E-6
+    @test tbl.f[3] ≈ 0.328551 atol = 1.0E-6
+    @test tbl.f[4] ≈ 0.106973 atol = 1.0E-6
+    @test tbl.p[2] ≈ 0.322206 atol = 1.0E-6
+    @test tbl.p[3] ≈ 0.570520 atol = 1.0E-6
+    @test tbl.p[4] ≈ 0.745747 atol = 1.0E-6
+
+end

From c9b2da5f1e0580ce14aa0884611026e190384c8a Mon Sep 17 00:00:00 2001
From: PharmCat <13901158+PharmCat@users.noreply.github.com>
Date: Mon, 7 Nov 2022 16:08:28 +0300
Subject: [PATCH 05/14] fix zero df

---
 src/ftest.jl | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/src/ftest.jl b/src/ftest.jl
index 3e05a0cc..488ab44d 100644
--- a/src/ftest.jl
+++ b/src/ftest.jl
@@ -352,6 +352,8 @@ tname(t::AbstractTerm) = "$(t.sym)"
 tname(t::InteractionTerm) = join(tname.(t.terms), " & ")
 tname(t::InterceptTerm) = "(Intercept)"
 
+
+
 """
     typeiii(obj)
 
@@ -406,7 +408,11 @@ function typeiii(obj)
         # I think this is more efficient:
         F[i]    = mulαtβα(LB, pinv(Symmetric(θ))) / RL
         df[i]   = RL
-        pval[i] = ccdf(FDist(df[i], dof_residual(obj)), F[i])
+        if iszero(df[i])
+            pval[i] = NaN
+        else
+            pval[i] = ccdf(FDist(df[i], dof_residual(obj)), F[i])
+        end
     end
     if length(d) > 0
         deleteat!(fac, d)

From 0fd5f6462c40ea8848665c01ca72c27b2e4b5a5e Mon Sep 17 00:00:00 2001
From: PharmCat <13901158+PharmCat@users.noreply.github.com>
Date: Tue, 8 Nov 2022 18:35:24 +0300
Subject: [PATCH 06/14] use dot

---
 src/ftest.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/ftest.jl b/src/ftest.jl
index 488ab44d..e646abc2 100644
--- a/src/ftest.jl
+++ b/src/ftest.jl
@@ -406,7 +406,7 @@ function typeiii(obj)
         # Then F can be computed:
         # F[i]    = (LB' * pinv(Symmetric(θ)) * LB)/rank(L)
         # I think this is more efficient:
-        F[i]    = mulαtβα(LB, pinv(Symmetric(θ))) / RL
+        F[i]    = dot(LB, pinv(Symmetric(θ)), LB) / RL
         df[i]   = RL
         if iszero(df[i])
             pval[i] = NaN

From a8044d2260fcbdafff859d7aa8644f2cd51ad847 Mon Sep 17 00:00:00 2001
From: PharmCat <13901158+PharmCat@users.noreply.github.com>
Date: Tue, 8 Nov 2022 18:57:30 +0300
Subject: [PATCH 07/14] coeftables out

---
 src/ftest.jl     | 23 +++++++++++++----------
 test/runtests.jl | 12 ++++++------
 2 files changed, 19 insertions(+), 16 deletions(-)

diff --git a/src/ftest.jl b/src/ftest.jl
index e646abc2..70550f23 100644
--- a/src/ftest.jl
+++ b/src/ftest.jl
@@ -235,14 +235,14 @@ end
 # Tests of Between-Subjects Effects
 # Baset on F-statistics
 # L: The s×p full row rank matrix. The rows are estimable functions. s≥1 where p number of coefs
-
+#=
 struct GroupEffectsTable
     name::Vector{String}
     df::Vector{Float64}
     f::Vector{Float64}
     p::Vector{Float64}
 end
-
+=#
 """
 θ + A * B * A'
 
@@ -262,6 +262,7 @@ function mulαβαtinc!(θ::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix
     end
     θ
 end
+#=
 """
 a' * B * a
 
@@ -280,7 +281,7 @@ function mulαtβα(a::AbstractVector, B::AbstractMatrix{T}) where T
     end
     c
 end
-
+=#
 # See SPSS (GLM/UNIANOVA) and SAS (PROC GLM) documentation 
 # https://www.ibm.com/docs/en/spss-statistics/29.0.0?topic=effects-tests-between-subjects
 # L is a s×p matrix corresponding to plan-matrix of Factor
@@ -369,7 +370,7 @@ function typeiii(obj)
     d           = Vector{Int}(undef, 0)
     fac         = Vector{String}(undef, c)
     F           = Vector{Float64}(undef,c)
-    df          = Vector{Float64}(undef, c)
+    df          = Vector{Tuple{Float64, Float64}}(undef, c)
     pval        = Vector{Float64}(undef, c)
     for i = 1:c
         # Make L matrix
@@ -407,11 +408,11 @@ function typeiii(obj)
         # F[i]    = (LB' * pinv(Symmetric(θ)) * LB)/rank(L)
         # I think this is more efficient:
         F[i]    = dot(LB, pinv(Symmetric(θ)), LB) / RL
-        df[i]   = RL
-        if iszero(df[i])
+        df[i]   = (RL, dof_residual(obj))
+        if iszero(df[i][1])
             pval[i] = NaN
         else
-            pval[i] = ccdf(FDist(df[i], dof_residual(obj)), F[i])
+            pval[i] = ccdf(FDist(df[i][1], df[i][2]), F[i])
         end
     end
     if length(d) > 0
@@ -420,13 +421,15 @@ function typeiii(obj)
         deleteat!(df, d)
         deleteat!(pval, d)
     end
-    GroupEffectsTable(fac, df, F, pval)
-
+    #GroupEffectsTable(fac, df, F, pval)
+    CoefTable([df, F, pval], ["DF/DDF", "F", "Pr(>F)"], fac, 3, 2)
 end
 
+#=
 using PrettyTables
 # use PrettyTables because it is fastest way to make table
 function Base.show(io::IO, obj::GroupEffectsTable)
     mx = hcat(obj.name, obj.df, obj.f, obj.p)
     PrettyTables.pretty_table(io, mx; tf = PrettyTables.tf_compact, header = ["Name", "DF" ,"F" ,"Pval"])
-end
\ No newline at end of file
+end
+=#
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 99fa0070..7f2cfeba 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1622,11 +1622,11 @@ berds1.Subject = categorical(berds1.Subject)
     # Intercept not included in test
     ols  = lm(@formula(Var ~ Sequence+Period+Formulation), berds1)
     tbl = GLM.typeiii(ols)
-    @test tbl.f[2] ≈ 1.011001 atol = 1.0E-6
-    @test tbl.f[3] ≈ 0.328551 atol = 1.0E-6
-    @test tbl.f[4] ≈ 0.106973 atol = 1.0E-6
-    @test tbl.p[2] ≈ 0.322206 atol = 1.0E-6
-    @test tbl.p[3] ≈ 0.570520 atol = 1.0E-6
-    @test tbl.p[4] ≈ 0.745747 atol = 1.0E-6
+    @test tbl.cols[2][2] ≈ 1.011001 atol = 1.0E-6
+    @test tbl.cols[2][3] ≈ 0.328551 atol = 1.0E-6
+    @test tbl.cols[2][4] ≈ 0.106973 atol = 1.0E-6
+    @test tbl.cols[3][2] ≈ 0.322206 atol = 1.0E-6
+    @test tbl.cols[3][3] ≈ 0.570520 atol = 1.0E-6
+    @test tbl.cols[3][4] ≈ 0.745747 atol = 1.0E-6
 
 end

From 0aa9ac7536a164af7d18332b7839f51c77ffb1e2 Mon Sep 17 00:00:00 2001
From: PharmCat <13901158+PharmCat@users.noreply.github.com>
Date: Thu, 10 Nov 2022 12:24:20 +0300
Subject: [PATCH 08/14] remove PrettyTables

---
 Project.toml | 2 --
 1 file changed, 2 deletions(-)

diff --git a/Project.toml b/Project.toml
index 0bc32a67..83b8d844 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,7 +6,6 @@ version = "1.8.1"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
-PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
@@ -21,7 +20,6 @@ CSV = "0.7, 0.8"
 CategoricalArrays = "0.8, 0.9, 0.10"
 DataFrames = "0.22, 1"
 Distributions = "0.16, 0.17, 0.18, 0.19, 0.20, 0.21, 0.22, 0.23, 0.24, 0.25"
-PrettyTables = "1, 2"
 RDatasets = "0.5, 0.6, 0.7"
 Reexport = "0.1, 0.2, 1.0"
 SpecialFunctions = "0.6, 0.7, 0.8, 0.9, 0.10, 1, 2.0"

From 80b880881d91566cba4646474432a0a3f62f3322 Mon Sep 17 00:00:00 2001
From: PharmCat <13901158+PharmCat@users.noreply.github.com>
Date: Thu, 10 Nov 2022 12:39:59 +0300
Subject: [PATCH 09/14] tests

---
 test/runtests.jl | 15 +++++++++++++++
 1 file changed, 15 insertions(+)

diff --git a/test/runtests.jl b/test/runtests.jl
index 7f2cfeba..da5e1740 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1620,6 +1620,7 @@ berds1.Subject = categorical(berds1.Subject)
   /DESIGN=Sequence Period Formulation.
     =#
     # Intercept not included in test
+    # Basic model
     ols  = lm(@formula(Var ~ Sequence+Period+Formulation), berds1)
     tbl = GLM.typeiii(ols)
     @test tbl.cols[2][2] ≈ 1.011001 atol = 1.0E-6
@@ -1628,5 +1629,19 @@ berds1.Subject = categorical(berds1.Subject)
     @test tbl.cols[3][2] ≈ 0.322206 atol = 1.0E-6
     @test tbl.cols[3][3] ≈ 0.570520 atol = 1.0E-6
     @test tbl.cols[3][4] ≈ 0.745747 atol = 1.0E-6
+    # Zero intercep
+    ols  = lm(@formula(Var ~ 0+Sequence+Period+Formulation), berds1)
+    tbl = GLM.typeiii(ols)
 
+    @test tbl.cols[2][1] ≈ 1.011001 atol = 1.0E-6
+    @test tbl.cols[2][2] ≈ 0.328551 atol = 1.0E-6
+    @test tbl.cols[2][3] ≈ 0.106973 atol = 1.0E-6
+    @test tbl.cols[3][1] ≈ 0.322206 atol = 1.0E-6
+    @test tbl.cols[3][2] ≈ 0.570520 atol = 1.0E-6
+    @test tbl.cols[3][3] ≈ 0.745747 atol = 1.0E-6
+    # Crossed factors
+    ols  = lm(@formula(Var ~ 1+Sequence&Period), berds1)
+    tbl = GLM.typeiii(ols)
+    @test tbl.cols[2][2] ≈ 0.482175 atol = 1.0E-6
+    @test tbl.cols[3][2] ≈ 0.696996 atol = 1.0E-6
 end

From 9844ddf8209639e62f4894aafeda7aca3e2c389a Mon Sep 17 00:00:00 2001
From: PharmCat <13901158+PharmCat@users.noreply.github.com>
Date: Mon, 2 Jan 2023 00:03:44 +0300
Subject: [PATCH 10/14] update

---
 Project.toml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Project.toml b/Project.toml
index 83b8d844..fb5e9819 100644
--- a/Project.toml
+++ b/Project.toml
@@ -13,7 +13,7 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
-StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
+StatsModels = "857c15b0-5950-435b-a230-016a4feb6864"
 
 [compat]
 CSV = "0.7, 0.8"
@@ -26,7 +26,7 @@ SpecialFunctions = "0.6, 0.7, 0.8, 0.9, 0.10, 1, 2.0"
 StatsAPI = "1.4"
 StatsBase = "0.33.5"
 StatsFuns = "0.6, 0.7, 0.8, 0.9, 1.0"
-StatsModels = "0.6.23"
+StatsModels = "0.1.0"
 julia = "1"
 
 [extras]

From d6301183f5cd6da8ba54a296f1c341962789190f Mon Sep 17 00:00:00 2001
From: PharmCat <13901158+PharmCat@users.noreply.github.com>
Date: Thu, 26 Jan 2023 15:06:07 +0300
Subject: [PATCH 11/14] update, clean code

---
 src/ftest.jl | 65 ++++++++--------------------------------------------
 1 file changed, 10 insertions(+), 55 deletions(-)

diff --git a/src/ftest.jl b/src/ftest.jl
index 70550f23..e9975921 100644
--- a/src/ftest.jl
+++ b/src/ftest.jl
@@ -235,14 +235,6 @@ end
 # Tests of Between-Subjects Effects
 # Baset on F-statistics
 # L: The s×p full row rank matrix. The rows are estimable functions. s≥1 where p number of coefs
-#=
-struct GroupEffectsTable
-    name::Vector{String}
-    df::Vector{Float64}
-    f::Vector{Float64}
-    p::Vector{Float64}
-end
-=#
 """
 θ + A * B * A'
 
@@ -251,37 +243,21 @@ Change θ (only upper triangle). B is symmetric.
 function mulαβαtinc!(θ::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix)
     axb  = axes(B, 1)
     sa   = size(A, 1)
-    for m ∈ 1:sa
-        for n ∈ m:sa
-            for j ∈ axb
-                @inbounds for i ∈ axb
-                    θ[m, n] +=  A[m, i] * B[i, j] * A[n, j]
+    for j ∈ axb
+        for i ∈ axb
+            @inbounds Bij = B[i, j]
+            for n ∈ 1:sa
+                @inbounds Anj = A[n, j]
+                BijAnj = Bij * Anj
+                @simd for m ∈ 1:n
+                    @inbounds θ[m, n] +=  A[m, i] * BijAnj
                 end
             end
         end
     end
     θ
 end
-#=
-"""
-a' * B * a
 
-"""
-function mulαtβα(a::AbstractVector, B::AbstractMatrix{T}) where T
-    if length(a) != size(B, 2)::Int  || size(B, 2)::Int  != size(B, 1)::Int  error("Dimention error") end
-    axbm  = axes(B, 1)
-    axbn  = axes(B, 2)
-    c = zero(T)
-    for i ∈ axbm
-        ct = zero(T)
-        for j ∈ axbn
-            @inbounds  ct += B[i, j] * a[j]
-        end
-        @inbounds c += ct * a[i]
-    end
-    c
-end
-=#
 # See SPSS (GLM/UNIANOVA) and SAS (PROC GLM) documentation 
 # https://www.ibm.com/docs/en/spss-statistics/29.0.0?topic=effects-tests-between-subjects
 # L is a s×p matrix corresponding to plan-matrix of Factor
@@ -307,7 +283,7 @@ end
 # 1 0  0  0  0  0
 # 1 0  0  0  0  0
 #
-# Then yoy wil have L matrix fo intercept:
+# Then you wil have L matrix for intercept:
 #
 # 1 0  0  0  0  0 
 #
@@ -353,8 +329,6 @@ tname(t::AbstractTerm) = "$(t.sym)"
 tname(t::InteractionTerm) = join(tname.(t.terms), " & ")
 tname(t::InterceptTerm) = "(Intercept)"
 
-
-
 """
     typeiii(obj)
 
@@ -375,14 +349,6 @@ function typeiii(obj)
     for i = 1:c
         # Make L matrix
         L       = lcontrast(obj, i)
-        #=
-        if typeof(obj.mf.f.rhs.terms[i]) <: InterceptTerm{true}
-            # I think L matrix for Intercept should represent general mean
-            # For this - L element corresponding to factor should be devided on total number of levels (not coefs)
-            # But because I can't get efficient number of levels for InteractionTerm
-            # I cant't calc it correctly
-        end
-        =#
         if typeof(obj.mf.f.rhs.terms[i]) <: InterceptTerm{false} # If zero intercept (drop)
             push!(d, i)
             fac[i] = ""
@@ -406,7 +372,6 @@ function typeiii(obj)
         LB = L * B
         # Then F can be computed:
         # F[i]    = (LB' * pinv(Symmetric(θ)) * LB)/rank(L)
-        # I think this is more efficient:
         F[i]    = dot(LB, pinv(Symmetric(θ)), LB) / RL
         df[i]   = (RL, dof_residual(obj))
         if iszero(df[i][1])
@@ -421,15 +386,5 @@ function typeiii(obj)
         deleteat!(df, d)
         deleteat!(pval, d)
     end
-    #GroupEffectsTable(fac, df, F, pval)
     CoefTable([df, F, pval], ["DF/DDF", "F", "Pr(>F)"], fac, 3, 2)
-end
-
-#=
-using PrettyTables
-# use PrettyTables because it is fastest way to make table
-function Base.show(io::IO, obj::GroupEffectsTable)
-    mx = hcat(obj.name, obj.df, obj.f, obj.p)
-    PrettyTables.pretty_table(io, mx; tf = PrettyTables.tf_compact, header = ["Name", "DF" ,"F" ,"Pval"])
-end
-=#
\ No newline at end of file
+end
\ No newline at end of file

From 6ae9c1436e0d752a41233057f6ce677a94950bf3 Mon Sep 17 00:00:00 2001
From: PharmCat <13901158+PharmCat@users.noreply.github.com>
Date: Thu, 26 Jan 2023 15:11:53 +0300
Subject: [PATCH 12/14] update

---
 src/ftest.jl | 65 ++++++++--------------------------------------------
 1 file changed, 10 insertions(+), 55 deletions(-)

diff --git a/src/ftest.jl b/src/ftest.jl
index 70550f23..e9975921 100644
--- a/src/ftest.jl
+++ b/src/ftest.jl
@@ -235,14 +235,6 @@ end
 # Tests of Between-Subjects Effects
 # Baset on F-statistics
 # L: The s×p full row rank matrix. The rows are estimable functions. s≥1 where p number of coefs
-#=
-struct GroupEffectsTable
-    name::Vector{String}
-    df::Vector{Float64}
-    f::Vector{Float64}
-    p::Vector{Float64}
-end
-=#
 """
 θ + A * B * A'
 
@@ -251,37 +243,21 @@ Change θ (only upper triangle). B is symmetric.
 function mulαβαtinc!(θ::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix)
     axb  = axes(B, 1)
     sa   = size(A, 1)
-    for m ∈ 1:sa
-        for n ∈ m:sa
-            for j ∈ axb
-                @inbounds for i ∈ axb
-                    θ[m, n] +=  A[m, i] * B[i, j] * A[n, j]
+    for j ∈ axb
+        for i ∈ axb
+            @inbounds Bij = B[i, j]
+            for n ∈ 1:sa
+                @inbounds Anj = A[n, j]
+                BijAnj = Bij * Anj
+                @simd for m ∈ 1:n
+                    @inbounds θ[m, n] +=  A[m, i] * BijAnj
                 end
             end
         end
     end
     θ
 end
-#=
-"""
-a' * B * a
 
-"""
-function mulαtβα(a::AbstractVector, B::AbstractMatrix{T}) where T
-    if length(a) != size(B, 2)::Int  || size(B, 2)::Int  != size(B, 1)::Int  error("Dimention error") end
-    axbm  = axes(B, 1)
-    axbn  = axes(B, 2)
-    c = zero(T)
-    for i ∈ axbm
-        ct = zero(T)
-        for j ∈ axbn
-            @inbounds  ct += B[i, j] * a[j]
-        end
-        @inbounds c += ct * a[i]
-    end
-    c
-end
-=#
 # See SPSS (GLM/UNIANOVA) and SAS (PROC GLM) documentation 
 # https://www.ibm.com/docs/en/spss-statistics/29.0.0?topic=effects-tests-between-subjects
 # L is a s×p matrix corresponding to plan-matrix of Factor
@@ -307,7 +283,7 @@ end
 # 1 0  0  0  0  0
 # 1 0  0  0  0  0
 #
-# Then yoy wil have L matrix fo intercept:
+# Then you wil have L matrix for intercept:
 #
 # 1 0  0  0  0  0 
 #
@@ -353,8 +329,6 @@ tname(t::AbstractTerm) = "$(t.sym)"
 tname(t::InteractionTerm) = join(tname.(t.terms), " & ")
 tname(t::InterceptTerm) = "(Intercept)"
 
-
-
 """
     typeiii(obj)
 
@@ -375,14 +349,6 @@ function typeiii(obj)
     for i = 1:c
         # Make L matrix
         L       = lcontrast(obj, i)
-        #=
-        if typeof(obj.mf.f.rhs.terms[i]) <: InterceptTerm{true}
-            # I think L matrix for Intercept should represent general mean
-            # For this - L element corresponding to factor should be devided on total number of levels (not coefs)
-            # But because I can't get efficient number of levels for InteractionTerm
-            # I cant't calc it correctly
-        end
-        =#
         if typeof(obj.mf.f.rhs.terms[i]) <: InterceptTerm{false} # If zero intercept (drop)
             push!(d, i)
             fac[i] = ""
@@ -406,7 +372,6 @@ function typeiii(obj)
         LB = L * B
         # Then F can be computed:
         # F[i]    = (LB' * pinv(Symmetric(θ)) * LB)/rank(L)
-        # I think this is more efficient:
         F[i]    = dot(LB, pinv(Symmetric(θ)), LB) / RL
         df[i]   = (RL, dof_residual(obj))
         if iszero(df[i][1])
@@ -421,15 +386,5 @@ function typeiii(obj)
         deleteat!(df, d)
         deleteat!(pval, d)
     end
-    #GroupEffectsTable(fac, df, F, pval)
     CoefTable([df, F, pval], ["DF/DDF", "F", "Pr(>F)"], fac, 3, 2)
-end
-
-#=
-using PrettyTables
-# use PrettyTables because it is fastest way to make table
-function Base.show(io::IO, obj::GroupEffectsTable)
-    mx = hcat(obj.name, obj.df, obj.f, obj.p)
-    PrettyTables.pretty_table(io, mx; tf = PrettyTables.tf_compact, header = ["Name", "DF" ,"F" ,"Pval"])
-end
-=#
\ No newline at end of file
+end
\ No newline at end of file

From 99c73ba48b572e691cb3285e6dfa031c5a977a78 Mon Sep 17 00:00:00 2001
From: PharmCat <13901158+PharmCat@users.noreply.github.com>
Date: Thu, 26 Jan 2023 15:26:30 +0300
Subject: [PATCH 13/14] fix

---
 Project.toml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Project.toml b/Project.toml
index fb5e9819..83b8d844 100644
--- a/Project.toml
+++ b/Project.toml
@@ -13,7 +13,7 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
-StatsModels = "857c15b0-5950-435b-a230-016a4feb6864"
+StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
 
 [compat]
 CSV = "0.7, 0.8"
@@ -26,7 +26,7 @@ SpecialFunctions = "0.6, 0.7, 0.8, 0.9, 0.10, 1, 2.0"
 StatsAPI = "1.4"
 StatsBase = "0.33.5"
 StatsFuns = "0.6, 0.7, 0.8, 0.9, 1.0"
-StatsModels = "0.1.0"
+StatsModels = "0.6.23"
 julia = "1"
 
 [extras]

From ff4721938884172234d3f1042922a2e2d682f8bd Mon Sep 17 00:00:00 2001
From: PharmCat <13901158+PharmCat@users.noreply.github.com>
Date: Thu, 2 Feb 2023 16:46:14 +0300
Subject: [PATCH 14/14] update

---
 src/ftest.jl     | 49 +++++++++++++++++++++++++++++++++++++-----------
 test/runtests.jl | 17 ++++++++++++++++-
 2 files changed, 54 insertions(+), 12 deletions(-)

diff --git a/src/ftest.jl b/src/ftest.jl
index bb6d460c..94cf6b42 100644
--- a/src/ftest.jl
+++ b/src/ftest.jl
@@ -303,16 +303,42 @@ end
 L-contrast matrix for `i` fixed effect.
 """
 function lcontrast(obj, i::Int)
-    n = length(obj.mf.f.rhs.terms)
+    n = length(obj.formula.rhs.terms)
+    cn = length(coef(obj))
     if i > n || n < 1 error("Factor number out of range 1-$(n)") end
-    p    = size(obj.mm.m, 2) # number of coefs
-    inds = findall(x -> x==i, obj.mm.assign)
-    if typeof(obj.mf.f.rhs.terms[i]) <: CategoricalTerm
-        mxc   = zeros(size(obj.mf.f.rhs.terms[i].contrasts.matrix, 1), p)
+    term = obj.formula.rhs.terms[i]
+    prev = 0
+    if i > 1
+        for j = 1:i-1
+            prev += width(obj.formula.rhs.terms[j])
+        end
+    end
+    #=
+    if isa(term, CategoricalTerm)
+        cm = term.contrasts.matrix
+        mx = zeros(Float64, size(cm, 1), cn)
+        view(mx, :, prev+1:prev+width(term)) .= cm
+    elseif isa(term, InteractionTerm)
+        m = width(term)
+        mx = zeros(Float64, m, cn)
+        for j = 1:m 
+            mx[j, j+prev] = 1
+        end
+    else
+        mx = zeros(Float64, 1, cn)
+        mx[1, prev+1] = 1
+    end
+    mx
+    =#
+    
+    p    = length(coef(obj)) # number of coefs
+    inds = prev+1:prev+width(term)
+    if typeof(term) <: CategoricalTerm
+        mxc   = zeros(size(term.contrasts.matrix, 1), p)
         mxcv  = view(mxc, :, inds)
-        mxcv .= obj.mf.f.rhs.terms[i].contrasts.matrix
-        mx    = zeros(size(obj.mf.f.rhs.terms[i].contrasts.matrix, 1) - 1, p)
-        for i = 2:size(obj.mf.f.rhs.terms[i].contrasts.matrix, 1) # correct for zero-intercept model
+        mxcv .= term.contrasts.matrix
+        mx    = zeros(size(term.contrasts.matrix, 1) - 1, p)
+        for i = 2:size(term.contrasts.matrix, 1) # correct for zero-intercept model
             mx[i-1, :] .= mxc[i, :] - mxc[1, :]
         end
     else
@@ -322,6 +348,7 @@ function lcontrast(obj, i::Int)
         end
     end
     mx
+    
 end
 
 tname(t::AbstractTerm) = "$(t.sym)"
@@ -339,7 +366,7 @@ function typeiii(obj)
     V           = vcov(obj) 
     replace!(V, NaN => 0) # Some values can be NaN - replace it to zero
     B           = coef(obj)   
-    c           = length(obj.mf.f.rhs.terms)
+    c           = length(obj.formula.rhs.terms)
     d           = Vector{Int}(undef, 0)
     fac         = Vector{String}(undef, c)
     F           = Vector{Float64}(undef,c)
@@ -348,12 +375,12 @@ function typeiii(obj)
     for i = 1:c
         # Make L matrix
         L       = lcontrast(obj, i)
-        if typeof(obj.mf.f.rhs.terms[i]) <: InterceptTerm{false} # If zero intercept (drop)
+        if typeof(obj.formula.rhs.terms[i]) <: InterceptTerm{false} # If zero intercept (drop)
             push!(d, i)
             fac[i] = ""
             continue
         else
-            fac[i] = tname(obj.mf.f.rhs.terms[i])
+            fac[i] = tname(obj.formula.rhs.terms[i])
         end
         # For case when cofs is zero (or NaN) we reduce rank of L-matrix
         for c = 1:length(B)
diff --git a/test/runtests.jl b/test/runtests.jl
index 2a41149d..3cd91d90 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1854,6 +1854,7 @@ berds1.Subject = categorical(berds1.Subject)
   /DESIGN=Sequence Period Formulation.
     =#
     # Intercept not included in test
+    
     # Basic model
     ols  = lm(@formula(Var ~ Sequence+Period+Formulation), berds1)
     tbl = GLM.typeiii(ols)
@@ -1863,19 +1864,33 @@ berds1.Subject = categorical(berds1.Subject)
     @test tbl.cols[3][2] ≈ 0.322206 atol = 1.0E-6
     @test tbl.cols[3][3] ≈ 0.570520 atol = 1.0E-6
     @test tbl.cols[3][4] ≈ 0.745747 atol = 1.0E-6
+    #=
+    GLM Var BY Sequence Period Formulation 
+  /METHOD=SSTYPE(3) 
+  /INTERCEPT=EXCLUDE 
+  /PRINT PARAMETER 
+  /CRITERIA=ALPHA(.05) 
+  /DESIGN=Sequence Period Formulation.
+    =#
+
     # Zero intercep
     ols  = lm(@formula(Var ~ 0+Sequence+Period+Formulation), berds1)
     tbl = GLM.typeiii(ols)
-
     @test tbl.cols[2][1] ≈ 1.011001 atol = 1.0E-6
     @test tbl.cols[2][2] ≈ 0.328551 atol = 1.0E-6
     @test tbl.cols[2][3] ≈ 0.106973 atol = 1.0E-6
     @test tbl.cols[3][1] ≈ 0.322206 atol = 1.0E-6
     @test tbl.cols[3][2] ≈ 0.570520 atol = 1.0E-6
     @test tbl.cols[3][3] ≈ 0.745747 atol = 1.0E-6
+
     # Crossed factors
     ols  = lm(@formula(Var ~ 1+Sequence&Period), berds1)
     tbl = GLM.typeiii(ols)
     @test tbl.cols[2][2] ≈ 0.482175 atol = 1.0E-6
     @test tbl.cols[3][2] ≈ 0.696996 atol = 1.0E-6
+
+    # Crossed factors (zero - intercept)
+    ols  = lm(@formula(Var ~ 0+Sequence&Period), berds1)
+    tbl = GLM.typeiii(ols)
+    @test tbl.cols[2][1] ≈ 87.103976 atol = 1.0E-6
 end