diff --git a/NEWS.md b/NEWS.md
index fe401143..33f24ebc 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,61 @@
+# v0.8.0
+
+- Term collections are now `Vector{AbstractTerm}` instead of tuples (#354).
+  This is a breaking change for package developers; user-facing `@formula`
+  syntax and `modelcols`/`modelmatrix` are unchanged.  The motivation is
+  compile latency: with tuples, every distinct number and order of terms in a
+  formula triggered a fresh specialization of the whole `apply_schema`,
+  `modelcols`, `coefnames`, etc. pipeline, so fitting a model with a
+  never-seen formula cost ~0.1-0.8s of compilation.  With vectors, this drops
+  to milliseconds.
+
+  - `+` on terms returns a `Vector{AbstractTerm}`, and `&` distributes over
+    vectors of terms instead of tuples.
+
+  - `InteractionTerm` and `MatrixTerm` are no longer parametric: both store
+    their terms in a `terms::Vector{AbstractTerm}` field.  The constructors
+    still accept any iterable of terms (including tuples), so
+    `InteractionTerm((a, b))` keeps working.  Code that dispatched on the
+    element types, e.g. `InteractionTerm{<:NTuple{N,CategoricalTerm}}`, must
+    switch to a run-time check such as
+    `all(t -> t isa CategoricalTerm, it.terms)`.
+
+  - `TupleTerm` is removed.  Methods that took a `TupleTerm` should take an
+    `AbstractVector{<:AbstractTerm}` instead.  `TermOrTerms` is now
+    `Union{AbstractTerm, AbstractVector{<:AbstractTerm}}`.
+
+  - `collect_matrix_terms` returns a `Vector{AbstractTerm}` (with the
+    `MatrixTerm` first) in the mixed matrix/non-matrix case, and `modelcols`
+    on a vector of terms returns a vector of the per-term columns.
+
+  - `apply_schema` on a vector of terms applies the schema term by term, from
+    left to right, and combines the results with `+`, so duplicates are
+    dropped and terms that expand to several terms are flattened.
+
+  - An empty vector of terms is valid everywhere a vector of terms is: it
+    yields an empty `MatrixTerm` of width 0 whose `modelcols` is a matrix with
+    no columns, like `InterceptTerm{false}`.  Constructing an
+    `InteractionTerm` with no terms throws an `ArgumentError`.
+
+  - `FormulaTerm`, `InteractionTerm`, and `MatrixTerm` now have content-based
+    `==` and `hash`, since the default field-identity comparison no longer
+    holds with vector fields, and `FunctionTerm` gains a `hash` consistent
+    with its existing `==`.  As a consequence, formulas compare equal after
+    `apply_schema` when their terms do, and hash-based deduplication
+    (`unique`, `Set`) works for formulas containing function calls.
+
+- The right-hand side of a formula built with `~` or `@formula` is always a
+  `Vector{AbstractTerm}`, even when it contains a single term (#354).
+  Previously `@formula(y ~ x).rhs` was the bare term `x`; it is now `[x]`.
+  This removes the lone-term special case, so code that handled both a term
+  and a tuple on the right-hand side can handle a vector only.  The left-hand
+  side is still a bare term, and after `apply_schema` the right-hand side is
+  still collected into a single `MatrixTerm` when every term is a matrix
+  term, so the `(y, X)` returned by `modelcols(f, data)` is unchanged.  The
+  same applies to `+` (`a + a` is now `[a]` rather than `a`) and to
+  `apply_schema` on a vector of terms, which returns a vector even when only
+  one term remains.
+
 # v0.7.0
 
 - `FunctionTerm` rework (#183)
diff --git a/Project.toml b/Project.toml
index bc723cb4..be60736f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "StatsModels"
 uuid = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
-version = "0.7.10"
+version = "0.8.0"
 
 [deps]
 DataAPI = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
diff --git a/docs/src/formula.md b/docs/src/formula.md
index 16c2f462..533962d1 100644
--- a/docs/src/formula.md
+++ b/docs/src/formula.md
@@ -299,7 +299,7 @@ symbols or strings (`Term`) and numbers (`ConstantTerm`), which makes it easy to
 work with collections of mixed type:
 
 ```jldoctest 1
-julia> ts = term.((1, :a, "b"))
+julia> ts = term.([1, :a, "b"])
 1
 a(unknown)
 b(unknown)
diff --git a/docs/src/internals.md b/docs/src/internals.md
index b5035569..0e06002d 100644
--- a/docs/src/internals.md
+++ b/docs/src/internals.md
@@ -57,26 +57,27 @@ expression returned by the `@formula` macro is evaluated.  At this point, the
 julia> using StatsModels;
 
 julia> dump(Term(:a) & Term(:b))
-InteractionTerm{Tuple{Term, Term}}
-  terms: Tuple{Term, Term}
+InteractionTerm
+  terms: Array{AbstractTerm}((2,))
     1: Term
       sym: Symbol a
     2: Term
       sym: Symbol b
 
 julia> dump(Term(:a) + Term(:b))
-Tuple{Term, Term}
+Array{AbstractTerm}((2,))
   1: Term
     sym: Symbol a
   2: Term
     sym: Symbol b
 
 julia> dump(Term(:y) ~ Term(:a))
-FormulaTerm{Term, Term}
+FormulaTerm{Term, Vector{AbstractTerm}}
   lhs: Term
     sym: Symbol y
-  rhs: Term
-    sym: Symbol a
+  rhs: Array{AbstractTerm}((1,))
+    1: Term
+      sym: Symbol a
 ```
 
 !!! note
@@ -254,7 +255,7 @@ terms:
 
 * `Term`s become `ContinuousTerm`s or `CategoricalTerm`s
 * `ConstantTerm`s become `InterceptTerm`s
-* Tuples of terms become [`MatrixTerm`](@ref)s where appropriate to explicitly indicate
+* Vectors of terms become [`MatrixTerm`](@ref)s where appropriate to explicitly indicate
   they should be concatenated into a single model matrix
 * Any model-specific (context-specific) interpretation of the terms is made, including
   transforming calls to functions that have special meaning in particular
@@ -274,7 +275,7 @@ Predictors:
   b(unknown) & c(unknown)
 
 julia> typeof(f)
-FormulaTerm{Term, Tuple{ConstantTerm{Int64}, Term, Term, Term, InteractionTerm{Tuple{Term, Term}}}}
+FormulaTerm{Term, Vector{AbstractTerm}}
 
 julia> f = apply_schema(f, schema(f, df))
 FormulaTerm
@@ -288,7 +289,7 @@ Predictors:
   b(continuous) & c(DummyCoding:3→2)
 
 julia> typeof(f)
-FormulaTerm{ContinuousTerm{Float64}, MatrixTerm{Tuple{InterceptTerm{true}, ContinuousTerm{Float64}, ContinuousTerm{Float64}, CategoricalTerm{DummyCoding, Matrix{Float64}, 2}, InteractionTerm{Tuple{ContinuousTerm{Float64}, CategoricalTerm{DummyCoding, Matrix{Float64}, 2}}}}}}
+FormulaTerm{ContinuousTerm{Float64}, MatrixTerm}
 ```
 
 This transformation is done by calling `apply_schema(term, schema, modeltype)`
@@ -575,9 +576,8 @@ julia> poly(my_col, my_degree)
 poly(a, 3)
 
 julia> poly.([:a, :b], my_degree)
-2-element Vector{PolyTerm{Term, ConstantTerm{Int64}}}:
- poly(a, 3)
- poly(b, 3)
+poly(a, 3)
+poly(b, 3)
 ```
 
 These run-time `PolyTerm`s are "schema-less" though, and to be able to construct
@@ -610,7 +610,7 @@ regression as above (which used `@formula(y ~ 1 + poly(a, 2) + poly(b, 2)`), but
 with the predictor names and the polynomial degree stored in variables:
 
 ```jldoctest 1
-julia> poly_vars = (:a, :b); poly_deg = 2;
+julia> poly_vars = [:a, :b]; poly_deg = 2;
 
 julia> poly_formula = term(:y) ~ term(1) + poly.(poly_vars, poly_deg)
 FormulaTerm
diff --git a/docs/src/temporal_terms.md b/docs/src/temporal_terms.md
index 6a3cdc9a..f75e78a4 100644
--- a/docs/src/temporal_terms.md
+++ b/docs/src/temporal_terms.md
@@ -50,7 +50,7 @@ Predictors:
   lead(x, 2)
 
 julia> modelmatrix(f, df)
-5×3 reshape(::Matrix{Union{Missing, Int64}}, 5, 3) with eltype Union{Missing, Int64}:
+5×3 Matrix{Union{Missing, Int64}}:
   2   missing   6
   4   missing   8
   6  2         10
@@ -88,7 +88,7 @@ Predictors:
   lead(x, 2)
 
 julia> modelmatrix(f2, df)
-5×3 reshape(::Matrix{Union{Missing, Int64}}, 5, 3) with eltype Union{Missing, Int64}:
+5×3 Matrix{Union{Missing, Int64}}:
   2   missing   6
   4   missing   8
   6  2         10
diff --git a/src/modelframe.jl b/src/modelframe.jl
index c9fed27f..6151a7be 100644
--- a/src/modelframe.jl
+++ b/src/modelframe.jl
@@ -113,7 +113,7 @@ keyword arguments are passed to [`apply_schema`](@ref).
     [`modelcols`](@ref) pipeline directly
 
 """
-function StatsAPI.modelmatrix(t::Union{AbstractTerm, TupleTerm}, data;
+function StatsAPI.modelmatrix(t::TermOrTerms, data;
                                hints=Dict{Symbol,Any}(), mod::Type{M}=StatisticalModel) where M
     Tables.istable(data) ||
         throw(ArgumentError("expected data in a Table, got $(typeof(data))"))
diff --git a/src/schema.jl b/src/schema.jl
index d05b0ad7..73175c85 100644
--- a/src/schema.jl
+++ b/src/schema.jl
@@ -11,9 +11,9 @@
 terms(t::FormulaTerm) = union(terms(t.lhs), terms(t.rhs))
 terms(t::InteractionTerm) = terms(t.terms)
 terms(t::FunctionTerm) = mapreduce(terms, union, t.args)
-terms(t::AbstractTerm) = [t]
+terms(t::AbstractTerm) = AbstractTerm[t]
 terms(t::MatrixTerm) = terms(t.terms)
-terms(t::TupleTerm) = mapreduce(terms, union, t)
+terms(ts::AbstractVector{<:AbstractTerm}) = mapreduce(terms, union, ts, init=AbstractTerm[])
 
 needs_schema(::AbstractTerm) = true
 needs_schema(::ConstantTerm) = false
@@ -82,13 +82,13 @@ julia> ts = [Term(:x), Term(:y)];
 
 julia> schema(ts, d)
 StatsModels.Schema with 2 entries:
-  x => x
   y => y
+  x => x
 
 julia> schema(ts, d, Dict(:x => HelmertCoding()))
 StatsModels.Schema with 2 entries:
-  x => x
   y => y
+  x => x
 
 julia> schema(term(:y), d, Dict(:y => CategoricalTerm))
 StatsModels.Schema with 1 entry:
@@ -101,8 +101,8 @@ same in a container, but when printed alone are different:
 ```jldoctest 1
 julia> sch = schema(ts, d)
 StatsModels.Schema with 2 entries:
-  x => x
   y => y
+  x => x
 
 julia> term(:x)
 x(unknown)
@@ -117,16 +117,13 @@ y(continuous)
 schema(data, hints=Dict{Symbol,Any}()) = schema(columntable(data), hints)
 schema(dt::D, hints=Dict{Symbol,Any}()) where {D<:ColumnTable} =
     schema(Term.(collect(fieldnames(D))), dt, hints)
-schema(ts::AbstractVector{<:AbstractTerm}, data, hints::Dict{Symbol}) =
-    schema(ts, columntable(data), hints)
-
-# handle hints:
-schema(ts::AbstractVector{<:AbstractTerm}, dt::ColumnTable,
-       hints::Dict{Symbol}=Dict{Symbol,Any}()) =
-    sch = Schema(t=>concrete_term(t, dt, hints) for t in ts)
 
 schema(f::TermOrTerms, data, hints::Dict{Symbol}) =
-    schema(filter(needs_schema, terms(f)), data, hints)
+    schema(f, columntable(data), hints)
+
+# handle hints:
+schema(f::TermOrTerms, dt::ColumnTable, hints::Dict{Symbol}=Dict{Symbol,Any}()) =
+    Schema(t => concrete_term(t, dt, hints) for t in filter(needs_schema, terms(f)))
 
 schema(f::TermOrTerms, data) = schema(f, data, Dict{Symbol,Any}())
 
@@ -237,7 +234,16 @@ in _most_ cases, but cause method ambiguity in some.
 """
 apply_schema(t, schema) = apply_schema(t, schema, Nothing)
 apply_schema(t, schema, Mod::Type) = t
-apply_schema(terms::TupleTerm, schema, Mod::Type) = reduce(+, apply_schema.(terms, Ref(schema), Mod))
+# sequential (left-to-right) application over a vector of terms, combined with
+# `+` so that duplicates are dropped and nested vectors are flattened; always
+# returns a Vector{AbstractTerm}
+function apply_schema(terms::AbstractVector{<:AbstractTerm}, schema, Mod::Type)
+    out = AbstractTerm[]
+    for t in terms
+        out = out + apply_schema(t, schema, Mod)
+    end
+    return out
+end
 
 apply_schema(t::Term, schema::Schema, Mod::Type) = schema[t]
 apply_schema(ft::FormulaTerm, schema::Schema, Mod::Type) =
@@ -250,7 +256,7 @@ apply_schema(it::InteractionTerm, schema::Schema, Mod::Type) =
 apply_schema(t::Union{ContinuousTerm, CategoricalTerm}, schema::Schema, Mod::Type) =
     get(schema, term(t.sym), t)
 apply_schema(t::MatrixTerm, sch::Schema, Mod::Type) =
-    MatrixTerm(apply_schema.(t.terms, Ref(sch), Mod))
+    MatrixTerm(apply_schema(t.terms, sch, Mod))
 
 # TODO: special case this for <:RegressionModel ?
 function apply_schema(t::ConstantTerm, schema::Schema, Mod::Type)
@@ -390,7 +396,7 @@ has_schema(t::ConstantTerm) = false
 has_schema(t::Term) = false
 has_schema(t::Union{ContinuousTerm,CategoricalTerm}) = true
 has_schema(t::InteractionTerm) = all(has_schema(tt) for tt in t.terms)
-has_schema(t::TupleTerm) = all(has_schema(tt) for tt in t)
+has_schema(t::AbstractVector{<:AbstractTerm}) = all(has_schema(tt) for tt in t)
 has_schema(t::MatrixTerm) = has_schema(t.terms)
 has_schema(t::FormulaTerm) = has_schema(t.lhs) && has_schema(t.rhs)
 # FunctionTerms may always be transformed by apply_schema
@@ -503,8 +509,8 @@ end
 drop_term(from, to) = symequal(from, to) ? ConstantTerm(1) : from
 drop_term(from::FormulaTerm, to) = FormulaTerm(from.lhs, drop_term(from.rhs, to))
 drop_term(from::MatrixTerm, to) = MatrixTerm(drop_term(from.terms, to))
-drop_term(from::TupleTerm, to) =
-    tuple((t for t = from if !symequal(t, to))...)
+drop_term(from::AbstractVector{<:AbstractTerm}, to) =
+    AbstractTerm[t for t in from if !symequal(t, to)]
 function drop_term(from::InteractionTerm, t)
     terms = drop_term(from.terms, t)
     length(terms) > 1 ? InteractionTerm(terms) : terms[1]
@@ -537,7 +543,7 @@ The data variables that this term refers to.
 termvars(::AbstractTerm) = Symbol[]
 termvars(t::Union{Term, CategoricalTerm, ContinuousTerm}) = [t.sym]
 termvars(t::InteractionTerm) = mapreduce(termvars, union, t.terms)
-termvars(t::TupleTerm) = mapreduce(termvars, union, t, init=Symbol[])
+termvars(ts::AbstractVector{<:AbstractTerm}) = mapreduce(termvars, union, ts, init=Symbol[])
 termvars(t::MatrixTerm) = termvars(t.terms)
 termvars(t::FormulaTerm) = union(termvars(t.lhs), termvars(t.rhs))
 termvars(t::FunctionTerm) = mapreduce(termvars, union, t.args, init=Symbol[])
diff --git a/src/temporal_terms.jl b/src/temporal_terms.jl
index 3a8f25a3..dd15d6d0 100644
--- a/src/temporal_terms.jl
+++ b/src/temporal_terms.jl
@@ -26,7 +26,7 @@ struct LeadLagTerm{T<:AbstractTerm, F<:Union{typeof(lead), typeof(lag)}} <: Abst
     nsteps::Int
 end
 
-terms(t::LeadLagTerm) = (t.term, )
+terms(t::LeadLagTerm) = AbstractTerm[t.term]
 
 function apply_schema(t::FunctionTerm{F}, sch::Schema, ctx::Type) where F<:Union{typeof(lead), typeof(lag)}
     opname = string(nameof(F.instance))
diff --git a/src/terms.jl b/src/terms.jl
index 3617e28b..94fb2a62 100644
--- a/src/terms.jl
+++ b/src/terms.jl
@@ -1,6 +1,5 @@
 abstract type AbstractTerm end
-const TermOrTerms = Union{AbstractTerm, Tuple{AbstractTerm, Vararg{AbstractTerm}}}
-const TupleTerm = Tuple{TermOrTerms, Vararg{TermOrTerms}}
+const TermOrTerms = Union{AbstractTerm, AbstractVector{<:AbstractTerm}}
 
 Base.broadcastable(x::AbstractTerm) = Ref(x)
 
@@ -50,11 +49,19 @@ any type (captured by the type parameters).
 
 * `lhs::L`: The left-hand side (e.g., response)
 * `rhs::R`: The right-hand side (e.g., predictors)
+
+When a formula is built with `~` (or the [`@formula`](@ref) macro), the
+right-hand side is always a `Vector{AbstractTerm}`, even when it holds a single
+term.  After [`apply_schema`](@ref), the right-hand side is collected into a
+[`MatrixTerm`](@ref) (or a vector of terms when some terms are not matrix
+terms, see [`collect_matrix_terms`](@ref)).
 """
 struct FormulaTerm{L,R} <: AbstractTerm
     lhs::L
     rhs::R
 end
+Base.:(==)(a::FormulaTerm, b::FormulaTerm) = a.lhs == b.lhs && a.rhs == b.rhs
+Base.hash(t::FormulaTerm, h::UInt) = hash(t.rhs, hash(t.lhs, hash(:FormulaTerm, h)))
 
 """
     FunctionTerm{F,Args} <: AbstractTerm
@@ -87,25 +94,27 @@ Response:
 Predictors:
   (a,b)->log(1 + a + b)
 
-julia> typeof(f.rhs)
+julia> ft = only(f.rhs)
+(a,b)->log(1 + a + b)
+
+julia> typeof(ft)
 FunctionTerm{typeof(log), Vector{FunctionTerm{typeof(+), Vector{AbstractTerm}}}}
 
-julia> typeof(only(f.rhs.args))
+julia> typeof(only(ft.args))
 FunctionTerm{typeof(+), Vector{AbstractTerm}}
 
-julia> only(f.rhs.args).args
-3-element Vector{AbstractTerm}:
- 1
- a(unknown)
- b(unknown)
+julia> only(ft.args).args
+1
+a(unknown)
+b(unknown)
 
-julia> f.rhs.f(1 + 3 + 4)
+julia> ft.f(1 + 3 + 4)
 2.0794415416798357
 
-julia> modelcols(f.rhs, (a=3, b=4))
+julia> modelcols(ft, (a=3, b=4))
 2.0794415416798357
 
-julia> modelcols(f.rhs, (a=[3, 4], b=[4, 5]))
+julia> modelcols(ft, (a=[3, 4], b=[4, 5]))
 2-element Vector{Float64}:
  2.0794415416798357
  2.302585092994046
@@ -119,9 +128,10 @@ end
 width(::FunctionTerm) = 1
 
 Base.:(==)(a::FunctionTerm, b::FunctionTerm) = a.f == b.f && a.args == b.args && a.exorig == b.exorig
+Base.hash(t::FunctionTerm, h::UInt) = hash(t.exorig, hash(t.args, hash(t.f, h)))
 
 """
-    InteractionTerm{Ts} <: AbstractTerm
+    InteractionTerm <: AbstractTerm
 
 Represents an _interaction_ between two or more individual terms.
 
@@ -130,7 +140,7 @@ Generated by combining multiple `AbstractTerm`s with `&` (which is what calls to
 
 # Fields
 
-* `terms::Ts`: the terms that participate in the interaction.
+* `terms::Vector{AbstractTerm}`: the terms that participate in the interaction.
 
 # Example
 
@@ -161,13 +171,23 @@ julia> modelcols(t, d)
  0.0       4.33378
 
 julia> modelcols(t.terms, d)
-([1, 2, 3, 4, 5, 6, 7, 8, 9], [0.236781883208121, 0.9437409715735081, 0.4456708824294644, 0.7636794266904741, 0.14507148958283067, 0.021124039581375875, 0.15254507694061115, 0.617492416565387, 0.48153065407402607], [0.0 0.0; 1.0 0.0; … ; 1.0 0.0; 0.0 1.0])
+3-element Vector{Array}:
+ [1, 2, 3, 4, 5, 6, 7, 8, 9]
+ [0.236781883208121, 0.9437409715735081, 0.4456708824294644, 0.7636794266904741, 0.14507148958283067, 0.021124039581375875, 0.15254507694061115, 0.617492416565387, 0.48153065407402607]
+ [0.0 0.0; 1.0 0.0; … ; 1.0 0.0; 0.0 1.0]
 ```
 """
-struct InteractionTerm{Ts} <: AbstractTerm
-    terms::Ts
+struct InteractionTerm <: AbstractTerm
+    terms::Vector{AbstractTerm}
+    function InteractionTerm(terms)
+        ts = collect(AbstractTerm, terms)
+        isempty(ts) && throw(ArgumentError("an InteractionTerm needs at least one term"))
+        return new(ts)
+    end
 end
 width(ts::InteractionTerm) = prod(width(t) for t in ts.terms)
+Base.:(==)(a::InteractionTerm, b::InteractionTerm) = a.terms == b.terms
+Base.hash(t::InteractionTerm, h::UInt) = hash(t.terms, hash(:InteractionTerm, h))
 
 """
     InterceptTerm{HasIntercept} <: AbstractTerm
@@ -229,34 +249,41 @@ CategoricalTerm(sym::Symbol, contrasts::ContrastsMatrix{C,T}) where {C,T} =
     CategoricalTerm{C,T,length(contrasts.coefnames)}(sym, contrasts)
 
 """
-    MatrixTerm{Ts} <: AbstractTerm
+    MatrixTerm <: AbstractTerm
 
 A collection of terms that should be combined to produce a single numeric matrix.
 
-A matrix term is created by [`apply_schema`](@ref) from a tuple of terms using
+A matrix term is created by [`apply_schema`](@ref) from a vector of terms using
 [`collect_matrix_terms`](@ref), which pulls out all the terms that are matrix
 terms as determined by the trait function [`is_matrix_term`](@ref), which is
 true by default for all `AbstractTerm`s.
+
+# Fields
+
+* `terms::Vector{AbstractTerm}`: the terms that make up the matrix.
 """
-struct MatrixTerm{Ts<:TupleTerm} <: AbstractTerm
-    terms::Ts
+struct MatrixTerm <: AbstractTerm
+    terms::Vector{AbstractTerm}
+    MatrixTerm(terms) = new(collect(AbstractTerm, terms))
 end
-# wrap single terms in a tuple
-MatrixTerm(t::AbstractTerm) = MatrixTerm((t, ))
-width(t::MatrixTerm) = sum(width(tt) for tt in t.terms)
+# wrap single terms in a vector
+MatrixTerm(t::AbstractTerm) = MatrixTerm(AbstractTerm[t])
+width(t::MatrixTerm) = sum(width(tt) for tt in t.terms; init=0)
+Base.:(==)(a::MatrixTerm, b::MatrixTerm) = a.terms == b.terms
+Base.hash(t::MatrixTerm, h::UInt) = hash(t.terms, hash(:MatrixTerm, h))
 
 """
-    collect_matrix_terms(ts::TupleTerm)
-    collect_matrix_terms(t::AbstractTerm) = collect_matrix_term((t, ))
+    collect_matrix_terms(ts::AbstractVector{<:AbstractTerm})
+    collect_matrix_terms(t::AbstractTerm) = collect_matrix_terms([t])
 
 Depending on whether the component terms are matrix terms (meaning they have
 [`is_matrix_term(T) == true`](@ref is_matrix_term)), `collect_matrix_terms` will
 return
 
 1.  A single `MatrixTerm` (if all components are matrix terms)
-2.  A tuple of the components (if none of them are matrix terms)
-3.  A tuple of terms, with all matrix terms collected into a single `MatrixTerm`
-    in the first element of the tuple, and the remaining non-matrix terms passed
+2.  A vector of the components (if none of them are matrix terms)
+3.  A vector of terms, with all matrix terms collected into a single `MatrixTerm`
+    in the first element of the vector, and the remaining non-matrix terms passed
     through unchanged.
 
 By default all terms are matrix terms (that is,
@@ -267,19 +294,18 @@ random effects terms in
 [MixedModels.jl](https://github.com/dmbates/MixedModels.jl).
 
 """
-function collect_matrix_terms(ts::TupleTerm)
-    ismat = collect(is_matrix_term.(ts))
+function collect_matrix_terms(ts::AbstractVector{<:AbstractTerm})
+    ismat = [is_matrix_term(t) for t in ts]
     if all(ismat)
         MatrixTerm(ts)
     elseif any(ismat)
-        matterms = ts[ismat]
-        (MatrixTerm(ts[ismat]), ts[.!ismat]...)
+        AbstractTerm[MatrixTerm(ts[ismat]); ts[.!ismat]]
     else
         ts
     end
 end
 collect_matrix_terms(t::T) where {T<:AbstractTerm} =
-    is_matrix_term(T) ? MatrixTerm((t, )) : t
+    is_matrix_term(T) ? MatrixTerm(AbstractTerm[t]) : t
 collect_matrix_terms(t::MatrixTerm) = t
 
 
@@ -308,7 +334,8 @@ function Base.show(io::IO, mime::MIME"text/plain", term::AbstractTerm; prefix=""
     print(io, prefix, term)
 end
 
-function Base.show(io::IO, mime::MIME"text/plain", terms::TupleTerm; prefix=nothing)
+function Base.show(io::IO, mime::MIME"text/plain", terms::AbstractVector{<:AbstractTerm};
+                   prefix=nothing)
     for t in terms
         show(io, mime, t; prefix=something(prefix, ""))
         # ensure that there are newlines in between each term after the first
@@ -316,7 +343,7 @@ function Base.show(io::IO, mime::MIME"text/plain", terms::TupleTerm; prefix=noth
         prefix = something(prefix, '\n')
     end
 end
-Base.show(io::IO, terms::TupleTerm) = join(io, terms, " + ")
+Base.show(io::IO, terms::AbstractVector{<:AbstractTerm}) = join(io, terms, " + ")
 
 Base.show(io::IO, ::MIME"text/plain", t::Term; prefix="") =
     print(io, prefix, t.sym, "(unknown)")
@@ -388,31 +415,43 @@ Base.:&(a::ConstantTerm, b::ConstantTerm) = (validate_interaction(a); validate_i
 
 # associative rule
 Base.:&(it::InteractionTerm, term::AbstractTerm) =
-    term in it.terms ? it : InteractionTerm((it.terms..., term))
+    term in it.terms ? it : InteractionTerm(AbstractTerm[it.terms; term])
 Base.:&(term::AbstractTerm, it::InteractionTerm) =
-    term in it.terms ? it : InteractionTerm((term, it.terms...))
+    term in it.terms ? it : InteractionTerm(AbstractTerm[term; it.terms])
 Base.:&(a::InteractionTerm, b::InteractionTerm) =
-    InteractionTerm((union(a.terms, b.terms)..., ))
+    InteractionTerm(union(a.terms, b.terms))
 
 # distributive rule
-Base.:&(term::AbstractTerm, terms::TupleTerm) = term .& terms
-Base.:&(terms::TupleTerm, term::AbstractTerm) = terms .& term
-Base.:&(as::TupleTerm, bs::TupleTerm) = ((a & b for a in as for b in bs)..., )
+Base.:&(term::AbstractTerm, terms::AbstractVector{<:AbstractTerm}) =
+    AbstractTerm[term & t for t in terms]
+Base.:&(terms::AbstractVector{<:AbstractTerm}, term::AbstractTerm) =
+    AbstractTerm[t & term for t in terms]
+Base.:&(as::AbstractVector{<:AbstractTerm}, bs::AbstractVector{<:AbstractTerm}) =
+    AbstractTerm[a & b for a in as for b in bs]
 
 # + concatenates terms
 Base.:+(a::AbstractTerm) = a
-Base.:+(a::AbstractTerm, b::AbstractTerm) = a==b ? a : (a, b)
+Base.:+(a::AbstractTerm, b::AbstractTerm) = a==b ? AbstractTerm[a] : AbstractTerm[a, b]
 
 # associative rule for +
-Base.:+(as::TupleTerm, b::AbstractTerm) = b in as ? as : (as..., b)
-Base.:+(a::AbstractTerm, bs::TupleTerm) = a in bs ? bs : (a, bs...)
-Base.:+(as::TupleTerm, bs::TupleTerm) = (union(as, bs)..., )
+Base.:+(as::AbstractVector{<:AbstractTerm}, b::AbstractTerm) =
+    b in as ? as : AbstractTerm[as; b]
+Base.:+(a::AbstractTerm, bs::AbstractVector{<:AbstractTerm}) =
+    a in bs ? bs : AbstractTerm[a; bs]
+Base.:+(as::AbstractVector{<:AbstractTerm}, bs::AbstractVector{<:AbstractTerm}) =
+    union!(AbstractTerm[], as, bs)
+# needed to beat the specificity of elementwise +(::Array, ::Array) from Base
+Base.:+(as::Vector{<:AbstractTerm}, bs::Vector{<:AbstractTerm}) =
+    union!(AbstractTerm[], as, bs)
 
 # * expansion
 Base.:*(a::TermOrTerms, b::TermOrTerms) = a + b + a&b
 
-cleanup(terms::TupleTerm) = Tuple(sort!(unique!(collect(terms)), by=degree))
-cleanup(x) = x
+# the rhs of a formula built with `~` is always a Vector{AbstractTerm}, sorted
+# by degree; a lone term is wrapped in a vector
+cleanup(terms::AbstractVector{<:AbstractTerm}) =
+    sort!(unique!(AbstractTerm[terms;]), by=degree)
+cleanup(t::AbstractTerm) = AbstractTerm[t]
 
 degree(::ConstantTerm) = 0
 degree(::InterceptTerm) = 0
@@ -447,10 +486,11 @@ function modelcols(t, d::D) where D
 end
 
 """
-    modelcols(ts::NTuple{N, AbstractTerm}, data) where N
+    modelcols(ts::AbstractVector{<:AbstractTerm}, data)
 
-When a tuple of terms is provided, `modelcols` broadcasts over the individual
-terms.  To create a single matrix, wrap the tuple in a [`MatrixTerm`](@ref).
+When a vector of terms is provided, `modelcols` returns a vector of the columns
+generated by the individual terms.  To create a single matrix, wrap the vector
+in a [`MatrixTerm`](@ref).
 
 # Example
 
@@ -459,13 +499,16 @@ julia> using StableRNGs; rng = StableRNG(1);
 
 julia> d = (a = [1:9;], b = rand(rng, 9), c = repeat(["d","e","f"], 3));
 
-julia> ts = apply_schema(term.((:a, :b, :c)), schema(d))
+julia> ts = apply_schema(term.([:a, :b, :c]), schema(d))
 a(continuous)
 b(continuous)
 c(DummyCoding:3→2)
 
 julia> cols = modelcols(ts, d)
-([1, 2, 3, 4, 5, 6, 7, 8, 9], [0.5851946422124186, 0.07733793456911231, 0.7166282400543453, 0.3203570514066232, 0.6530930076222579, 0.2366391513734556, 0.7096838914472361, 0.5577872440804086, 0.05079002172175784], [0.0 0.0; 1.0 0.0; … ; 1.0 0.0; 0.0 1.0])
+3-element Vector{Array}:
+ [1, 2, 3, 4, 5, 6, 7, 8, 9]
+ [0.5851946422124186, 0.07733793456911231, 0.7166282400543453, 0.3203570514066232, 0.6530930076222579, 0.2366391513734556, 0.7096838914472361, 0.5577872440804086, 0.05079002172175784]
+ [0.0 0.0; 1.0 0.0; … ; 1.0 0.0; 0.0 1.0]
 
 julia> reduce(hcat, cols)
 9×4 Matrix{Float64}:
@@ -492,7 +535,7 @@ julia> modelcols(MatrixTerm(ts), d)
  9.0  0.05079    0.0  1.0
 ```
 """
-modelcols(ts::TupleTerm, d::NamedTuple) = modelcols.(ts, Ref(d))
+modelcols(ts::AbstractVector{<:AbstractTerm}, d::NamedTuple) = [modelcols(t, d) for t in ts]
 
 modelcols(t::Term, d::NamedTuple) = getproperty(d, t.sym)
 modelcols(t::ConstantTerm, d::NamedTuple) = t.n
@@ -552,12 +595,14 @@ modelcols(t::InterceptTerm{false}, d) = Matrix{Float64}(undef, size(first(d),1),
 modelcols(t::FormulaTerm, d::NamedTuple) = (modelcols(t.lhs,d), modelcols(t.rhs, d))
 
 function modelcols(t::MatrixTerm, d::ColumnTable)
+    # an empty MatrixTerm (e.g. from an empty rhs) has no columns, like InterceptTerm{false}
+    isempty(t.terms) && return Matrix{Float64}(undef, size(first(d), 1), 0)
     mat = reduce(hcat, [modelcols(tt, d) for tt in t.terms])
     reshape(mat, size(mat, 1), :)
 end
 
 modelcols(t::MatrixTerm, d::NamedTuple) =
-    reduce(vcat, [modelcols(tt, d) for tt in t.terms])
+    isempty(t.terms) ? Float64[] : reduce(vcat, [modelcols(tt, d) for tt in t.terms])
 
 vectorize(x::Tuple) = collect(x)
 vectorize(x::AbstractVector) = x
@@ -577,7 +622,7 @@ StatsAPI.coefnames(t::ContinuousTerm) = string(t.sym)
 StatsAPI.coefnames(t::CategoricalTerm) =
     ["$(t.sym): $name" for name in t.contrasts.coefnames]
 StatsAPI.coefnames(t::FunctionTerm) = string(t.exorig)
-StatsAPI.coefnames(ts::TupleTerm) = reduce(vcat, coefnames.(ts))
+StatsAPI.coefnames(ts::AbstractVector{<:AbstractTerm}) = mapreduce(coefnames, vcat, ts; init=String[])
 StatsAPI.coefnames(t::MatrixTerm) = mapreduce(coefnames, vcat, t.terms; init = String[])
 StatsAPI.coefnames(t::InteractionTerm) =
     kron_insideout((args...) -> join(args, " & "), vectorize.(coefnames.(t.terms))...)
@@ -642,9 +687,9 @@ termnames(t::CategoricalTerm) = string(t.sym)
 termnames(t::Term) = string(t.sym)
 termnames(t::ConstantTerm) = string(t.n)
 termnames(t::FunctionTerm) = string(t.exorig)
-# termnames(TupleTerm)) always returns a vector, even if it's just one element, e.g.,
-# termnames((term(:a),))
-termnames(ts::TupleTerm) = mapreduce(termnames, vcat, ts; init=String[])
+# termnames(::AbstractVector{<:AbstractTerm}) always returns a vector, even if
+# it's just one element, e.g., termnames([term(:a)])
+termnames(ts::AbstractVector{<:AbstractTerm}) = mapreduce(termnames, vcat, ts; init=String[])
 # termnames(MatrixTerm)) always returns a vector, even if it's just one element, e.g.,
 # termnames(MatrixTerm(term(:a)))
 termnames(t::MatrixTerm) = mapreduce(termnames, vcat, t.terms; init=String[])
@@ -656,7 +701,7 @@ termnames(t::InteractionTerm) =
 
 hasintercept(f::FormulaTerm) = hasintercept(f.rhs)
 hasintercept(t::AbstractTerm) = t == InterceptTerm{true}() || t == ConstantTerm(1)
-hasintercept(t::TupleTerm) = any(hasintercept, t)
+hasintercept(t::AbstractVector{<:AbstractTerm}) = any(hasintercept, t)
 hasintercept(t::MatrixTerm) = hasintercept(t.terms)
 
 omitsintercept(f::FormulaTerm) = omitsintercept(f.rhs)
@@ -664,7 +709,7 @@ omitsintercept(t::AbstractTerm) =
     t == InterceptTerm{false}() ||
     t == ConstantTerm(0) ||
     t == ConstantTerm(-1)
-omitsintercept(t::TupleTerm) = any(omitsintercept, t)
+omitsintercept(t::AbstractVector{<:AbstractTerm}) = any(omitsintercept, t)
 omitsintercept(t::MatrixTerm) = omitsintercept(t.terms)
 
 hasresponse(t) = false
@@ -684,13 +729,13 @@ are converted to symbols before wrapping.
 # Example
 
 ```jldoctest
-julia> ts = term.((1, :a, "b"))
+julia> ts = term.([1, :a, "b"])
 1
 a(unknown)
 b(unknown)
 
 julia> typeof(ts)
-Tuple{ConstantTerm{Int64}, Term, Term}
+Vector{AbstractTerm} (alias for Array{AbstractTerm, 1})
 ```
 """
 term(n::Number) = ConstantTerm(n)
diff --git a/src/vif.jl b/src/vif.jl
index 4fb699ec..8c701d77 100644
--- a/src/vif.jl
+++ b/src/vif.jl
@@ -17,8 +17,8 @@ _find_intercept(form::FormulaTerm) = _find_intercept(form.rhs)
 _find_intercept(::AbstractTerm) = nothing
 _find_intercept(::InterceptTerm{true}) = 1
 _find_intercept(m::MatrixTerm) = _find_intercept(m.terms)
-function _find_intercept(t::TupleTerm) 
-    return findfirst(!isnothing ∘ _find_intercept, t)
+function _find_intercept(ts::AbstractVector{<:AbstractTerm})
+    return findfirst(!isnothing ∘ _find_intercept, ts)
 end
 
 # borrowed from Effects.jl
@@ -97,7 +97,7 @@ function StatsAPI.gvif(model::RegressionModel; scale=false)
     tn = last(termnames(model))
     tn = view(tn, axes(tn, 1) .!= intercept)
     trms = get_matrix_term(form.rhs).terms
-    # MatrixTerms.terms is a tuple or vector so always 1-based indexing
+    # MatrixTerm.terms is a vector so always 1-based indexing
     trms = trms[1:length(trms) .!= intercept]
 
     df = width.(trms)
diff --git a/test/extension.jl b/test/extension.jl
index 288938df..ec82c004 100644
--- a/test/extension.jl
+++ b/test/extension.jl
@@ -37,8 +37,7 @@ end
 
         f_plain = apply_schema(f, sch)
         @test f_plain.rhs.terms[1] isa FunctionTerm
-        # this works but == is not defined correctly and apply_schema creates a new instance
-        @test_broken f_plain == apply_schema(f, sch, Nothing)
+        @test f_plain == apply_schema(f, sch, Nothing)
         @test last(modelcols(f_plain, d)) == hcat(d[:x].^3)
         
         f_special = apply_schema(f, sch, PolyModel)
@@ -55,12 +54,12 @@ end
 
         @test collect_matrix_terms(f.rhs) == MatrixTerm((term(:x) + term(:y)))
         @test collect_matrix_terms(f2.rhs) ==
-            (MatrixTerm((term(:x), )), NonMatrixTerm(term(:y)))
+            [MatrixTerm(term(:x)), NonMatrixTerm(term(:y))]
         @test collect_matrix_terms(f3.rhs) ==
-            (MatrixTerm((term(:y), )), NonMatrixTerm(term(:x)))
+            [MatrixTerm(term(:y)), NonMatrixTerm(term(:x))]
         @test collect_matrix_terms(f4.rhs) == f4.rhs
         @test collect_matrix_terms(f5.rhs) ==
-            (MatrixTerm((term(:x), term(:y))), NonMatrixTerm(term(:y)))
+            [MatrixTerm([term(:x), term(:y)]), NonMatrixTerm(term(:y))]
 
         f = apply_schema(f, sch)
         @test f.rhs isa MatrixTerm
@@ -68,27 +67,27 @@ end
         @test modelcols(f.rhs, d) == hcat(d.x, d.y)
 
         f2 = apply_schema(f2, sch)
-        @test f2.rhs isa Tuple{MatrixTerm, NonMatrixTerm}
-        @test f2.rhs == apply_schema((MatrixTerm(term(:x)), NonMatrixTerm(term(:y))), sch)
-        @test modelcols(f2.rhs, d) == (hcat(d.x), d.y)
+        @test f2.rhs isa Vector{AbstractTerm}
+        @test f2.rhs == apply_schema(AbstractTerm[MatrixTerm(term(:x)), NonMatrixTerm(term(:y))], sch)
+        @test modelcols(f2.rhs, d) == [hcat(d.x), d.y]
 
         # matrix term goes first
         f3 = apply_schema(f3, sch)
-        @test f3.rhs isa Tuple{MatrixTerm, NonMatrixTerm}
-        @test f3.rhs == apply_schema((MatrixTerm(term(:y)), NonMatrixTerm(term(:x))), sch)
-        @test modelcols(f3.rhs, d) == (hcat(d.y), d.x)
+        @test f3.rhs isa Vector{AbstractTerm}
+        @test f3.rhs == apply_schema(AbstractTerm[MatrixTerm(term(:y)), NonMatrixTerm(term(:x))], sch)
+        @test modelcols(f3.rhs, d) == [hcat(d.y), d.x]
 
         f4 = apply_schema(f4, sch)
-        @test f4.rhs isa Tuple{NonMatrixTerm, NonMatrixTerm}
-        @test f4.rhs == apply_schema((NonMatrixTerm(term(:x)), NonMatrixTerm(term(:y))), sch)
-        @test modelcols(f4.rhs, d) == (d.x, d.y)
+        @test all(t -> t isa NonMatrixTerm, f4.rhs)
+        @test f4.rhs == apply_schema([NonMatrixTerm(term(:x)), NonMatrixTerm(term(:y))], sch)
+        @test modelcols(f4.rhs, d) == [d.x, d.y]
 
         # matrix terms are gathered
         f5 = apply_schema(f5, sch)
-        @test f5.rhs isa Tuple{MatrixTerm, NonMatrixTerm}
+        @test f5.rhs isa Vector{AbstractTerm}
         @test f5.rhs ==
-            apply_schema((MatrixTerm((term.((:x, :y)))), NonMatrixTerm(term(:y))), sch)
-        @test modelcols(f5.rhs, d) == (hcat(d.x, d.y), d.y)
+            apply_schema(AbstractTerm[MatrixTerm(term.([:x, :y])), NonMatrixTerm(term(:y))], sch)
+        @test modelcols(f5.rhs, d) == [hcat(d.x, d.y), d.y]
         
     end
 
diff --git a/test/formula.jl b/test/formula.jl
index 4e1a4b93..f543f03d 100644
--- a/test/formula.jl
+++ b/test/formula.jl
@@ -7,7 +7,7 @@
     @test !hasresponse(t)
     @test !hasintercept(t)
     @test omitsintercept(t)
-    @test t.rhs == ConstantTerm(0)
+    @test t.rhs == [ConstantTerm(0)]
     @test issetequal(terms(t), [ConstantTerm(0)])
 
     ## empty lhs, intercept on rhs
@@ -21,7 +21,7 @@
     @test hasintercept(t) == false
     @test omitsintercept(t) == true
     @test hasresponse(t)
-    @test t.rhs == ConstantTerm(0)
+    @test t.rhs == [ConstantTerm(0)]
     @test issetequal(terms(t), term.((:y, 0)))
 
     t = @formula(y ~ -1)
@@ -32,59 +32,59 @@
     t = @formula(y ~ 1)
     @test hasresponse(t) == true
     @test hasintercept(t) == true
-    @test t.rhs == onet
+    @test t.rhs == [onet]
     @test issetequal(terms(t), (onet, y))
 
     ## terms add
     t = @formula(y ~ 1 + x1 + x2)
     @test hasintercept(t) == true
-    @test t.rhs == (onet, x1, x2)
+    @test t.rhs == [onet, x1, x2]
     @test issetequal(terms(t), [y, onet, x1, x2])
 
     ## implicit intercept behavior: NO intercept after @formula
     t = @formula(y ~ x1 + x2)
     @test hasintercept(t) == false
     @test omitsintercept(t) == false
-    @test t.rhs == (x1, x2)
+    @test t.rhs == [x1, x2]
     @test issetequal(terms(t), [y, x1, x2])
 
     ## no intercept
     t = @formula(y ~ 0 + x1 + x2)
     @test hasintercept(t) == false
     @test omitsintercept(t) == true
-    @test t.rhs == term.((0, :x1, :x2))
+    @test t.rhs == term.([0, :x1, :x2])
 
     t = @formula(y ~ -1 + x1 + x2)
     @test hasintercept(t) == false
     @test omitsintercept(t) == true
-    @test t.rhs == term.((-1, :x1, :x2))
+    @test t.rhs == term.([-1, :x1, :x2])
 
     t = @formula(y ~ x1 & x2)
     @test hasintercept(t) == false
     @test omitsintercept(t) == false
-    @test t.rhs == x1&x2
+    @test t.rhs == [x1&x2]
     @test issetequal(terms(t), [y, x1, x2])
 
     ## `*` expansion
     t = @formula(y ~ x1 * x2)
     @test hasintercept(t) == false
     @test omitsintercept(t) == false
-    @test t.rhs == (x1, x2, x1&x2)
+    @test t.rhs == [x1, x2, x1&x2]
     @test issetequal(terms(t), [y, x1, x2])
 
     ## associative rule:
     ## +
     t = @formula(y ~ x1 + x2 + x3)
-    @test t.rhs == (x1, x2, x3)
+    @test t.rhs == [x1, x2, x3]
 
     ## &
     t = @formula(y ~ x1 & x2 & x3)
-    @test t.rhs == x1&x2&x3
+    @test t.rhs == [x1&x2&x3]
     @test issetequal(terms(t), [y, x1, x2, x3])
 
     ## distributive property of + and &
     t = @formula(y ~ x1 & (x2 + x3))
-    @test t.rhs == (x1&x2, x1&x3)
+    @test t.rhs == [x1&x2, x1&x3]
     @test issetequal(terms(t), [y, x1, x2, x3])
     
     ## ordering of interaction terms is preserved across distributive
@@ -93,19 +93,19 @@
 
     ## distributive with *
     t = @formula(y ~ (a + b) * c)
-    @test t.rhs == (a, b, c, a&c, b&c)
+    @test t.rhs == [a, b, c, a&c, b&c]
 
     ## three-way *
     t = @formula(y ~ a * b * c)
-    @test t.rhs == (a, b, c, a&b, a&c, b&c, a&b&c)
+    @test t.rhs == [a, b, c, a&b, a&c, b&c, a&b&c]
     @test issetequal(terms(t), (y, a, b, c))
 
     ## Interactions with `1` reduce to main effect.
     t = @formula(y ~ 1 & x1)
-    @test t.rhs == x1
+    @test t.rhs == [x1]
 
     t = @formula(y ~ (1 + x1) & x2)
-    @test t.rhs == (x2, x1&x2)
+    @test t.rhs == [x2, x1&x2]
 
     ## PR #54 breaks formula-level equality because original (un-lowered)
     ## expression is kept on Formula struct.  but functional (RHS) equality
diff --git a/test/modelmatrix.jl b/test/modelmatrix.jl
index cd92bf9d..91d44c9b 100644
--- a/test/modelmatrix.jl
+++ b/test/modelmatrix.jl
@@ -393,7 +393,9 @@
     @testset "#136" begin
         t = (x = rand(100), y = randn(100));
         f = @formula(y ~ x)
-        @test modelcols(f, t) === (t.y, t.x)
+        y, x = modelcols(f, t)
+        @test y === t.y
+        @test only(x) === t.x
     end
 
     @testset "#185 - interactions of scalar terms for row tables" begin
diff --git a/test/protect.jl b/test/protect.jl
index 2fc9c487..ffe616eb 100644
--- a/test/protect.jl
+++ b/test/protect.jl
@@ -16,7 +16,7 @@
     end
 
     @testset "unprotect" begin
-        using StatsModels: TupleTerm, FullRank
+        using StatsModels: FullRank
         # unprotect reverts to treating calls to +, &, and * as term union,
         # interaction, and combined
 
@@ -26,9 +26,9 @@
         sch = schema(d)
         a, b, c = apply_schema.(term.((:a, :b, :c)), Ref(sch))
 
-        ops_types = ((+) => TupleTerm,
+        ops_types = ((+) => Vector{AbstractTerm},
                      (&) => InteractionTerm,
-                     (*) => TupleTerm,
+                     (*) => Vector{AbstractTerm},
                      (~) => FormulaTerm)
 
         for (op, typ) in ops_types, sch in (sch, FullRank(sch))
@@ -45,7 +45,7 @@
 
         # stops once it hits an non-special call still
         f = ft(+, a, ft(log, ft(+, term(1), b)))
-        @test apply_schema(f, sch) == (a, f.args[2])
+        @test apply_schema(f, sch) == [a, f.args[2]]
 
         # testing behavior of modelcols
         f = @formula(0 ~ 1 - unprotect(a&b))
@@ -54,21 +54,21 @@
 
         # ideally you'd also be able to do these but it's hard to make it work...
         # this fails because - doesn't auto-broadcast, and the broadcasting that
-        # happens during FunctionTerm evaluation gets used up by the (a,b) tuple.
+        # happens during FunctionTerm evaluation gets used up by the [a,b] vector.
         f = @formula(0 ~ 1 - unprotect(a+b))
-        @test f.rhs.args[end] isa FunctionTerm{typeof(unprotect)}
+        @test only(f.rhs).args[end] isa FunctionTerm{typeof(unprotect)}
         ff = apply_schema(f, schema(d))
-        @test ff.rhs.terms[1].args[end] isa StatsModels.TupleTerm
+        @test ff.rhs.terms[1].args[end] isa Vector{AbstractTerm}
         @test_broken modelcols(ff.rhs, d) == 1 .- [d.a d.b]
 
         # and even if we define a broadcasting version, still fails because it
-        # gives a tuple of arrays instead of a matrix
+        # gives a vector of arrays instead of a matrix
         my_sub = (x,y) -> x .- y
         ff = apply_schema(@formula(0 ~ my_sub(1, unprotect(a+b))), schema(d))
         @test_broken modelcols(ff.rhs, d) == 1 .- [d.a d.b]
         
         # both of these could be fixed by always returning a matrix when you call
-        # modelcols on a tuple of terms but that would break other things
+        # modelcols on a vector of terms but that would break other things
     end
 
 end
diff --git a/test/schema.jl b/test/schema.jl
index 9f6cf7f3..8d371f46 100644
--- a/test/schema.jl
+++ b/test/schema.jl
@@ -9,10 +9,10 @@
         @test f == apply_schema(f, schema(f, df))
     end
 
-    @testset "lonely term in a tuple" begin
+    @testset "lonely term in a vector" begin
         d = (a = [1,1],)
-        @test apply_schema(ConstantTerm(1), schema(d)) == apply_schema((ConstantTerm(1),), schema(d))
-        @test apply_schema(Term(:a), schema(d)) == apply_schema((Term(:a),), schema(d))
+        @test [apply_schema(ConstantTerm(1), schema(d))] == apply_schema([ConstantTerm(1)], schema(d))
+        @test [apply_schema(Term(:a), schema(d))] == apply_schema([Term(:a)], schema(d))
     end
 
     @testset "hints" begin
diff --git a/test/statsmodel.jl b/test/statsmodel.jl
index e7f751c3..3f8c8220 100644
--- a/test/statsmodel.jl
+++ b/test/statsmodel.jl
@@ -271,9 +271,9 @@ end
     @test termnames(FunctionTerm(log, [Term(:x)], :(log(x)))) == "log(x)"
     @test termnames(InteractionTerm(term.((:a, :b, :c)))) == "a & b & c"
     @test termnames(MatrixTerm(term(:a))) == ["a"]
-    @test termnames(MatrixTerm((term(:a), term(:b)))) == ["a", "b"]
-    @test termnames((term(:a), term(:b))) == ["a", "b"]
-    @test termnames((term(:a),)) == ["a"]
+    @test termnames(MatrixTerm([term(:a), term(:b)])) == ["a", "b"]
+    @test termnames([term(:a), term(:b)]) == ["a", "b"]
+    @test termnames([term(:a)]) == ["a"]
 end
 
 @testset "lrtest" begin
diff --git a/test/terms.jl b/test/terms.jl
index b58900bc..8813e859 100644
--- a/test/terms.jl
+++ b/test/terms.jl
@@ -6,7 +6,7 @@ end
 mimestring(x) = mimestring(MIME"text/plain", x)
 
 struct MultiTerm <: AbstractTerm
-    terms::StatsModels.TupleTerm
+    terms::Vector{AbstractTerm}
 end
 StatsModels.apply_schema(mt::MultiTerm, sch::StatsModels.Schema, Mod::Type) =
     apply_schema.(mt.terms, Ref(sch), Mod)
@@ -55,8 +55,8 @@ StatsModels.apply_schema(mt::MultiTerm, sch::StatsModels.Schema, Mod::Type) =
     @testset "term operators" begin
         a = term(:a)
         b = term(:b)
-        @test a + b == (a, b)
-        @test (a ~ b) == FormulaTerm(a, b)
+        @test a + b == [a, b]
+        @test (a ~ b) == FormulaTerm(a, [b])
         @test string(a~b) == "$a ~ $b"
         @test mimestring(a~b) ==
             """FormulaTerm
@@ -136,19 +136,24 @@ StatsModels.apply_schema(mt::MultiTerm, sch::StatsModels.Schema, Mod::Type) =
 
         ## addition of two identical function terms
         @test f2.rhs + f2.rhs == f2.rhs
+
+        ## hash is consistent with ==, so hash-based deduplication works
+        @test hash(f1.rhs) == hash(f2.rhs)
+        @test hash(f1) == hash(f2)
+        @test length(unique([only(f2.rhs), only(f2.rhs)])) == 1
+        @test length(Set([f1, f2])) == 1
+        @test length(Set([term(:z) & only(f2.rhs), term(:z) & only(f2.rhs)])) == 1
     end
 
-    @testset "expand nested tuples of terms during apply_schema" begin
+    @testset "flatten term vectors during apply_schema" begin
         sch = schema((a=rand(10), b=rand(10), c=rand(10)))
 
-        # nested tuples of terms are expanded by apply_schema
-        terms = (term(:a), (term(:b), term(:c)))
-        terms2 = apply_schema(terms, sch, Nothing)
-        @test terms2 isa NTuple{3, ContinuousTerm}
-        @test terms2 == apply_schema(term.((:a, :b, :c)), sch, Nothing)
+        terms2 = apply_schema(term.([:a, :b, :c]), sch, Nothing)
+        @test terms2 isa Vector{AbstractTerm}
+        @test all(t -> t isa ContinuousTerm, terms2)
 
-        # a term that generates multiple terms after apply_schema
-        mterms = (terms[1], MultiTerm(terms[2]))
+        # a term that generates multiple terms after apply_schema is flattened
+        mterms = AbstractTerm[term(:a), MultiTerm([term(:b), term(:c)])]
         terms3 = apply_schema(mterms, sch, Nothing)
 
         @test terms2 == terms3
@@ -231,41 +236,208 @@ StatsModels.apply_schema(mt::MultiTerm, sch::StatsModels.Schema, Mod::Type) =
 
     end
 
-    @testset "Tuple terms" begin
-        using StatsModels: TermOrTerms, TupleTerm, Term
+    @testset "Term containers" begin
+        using StatsModels: TermOrTerms, Term
         a, b, c = Term.((:a, :b, :c))
 
-        # TermOrTerms - one or more AbstractTerms (if more, a tuple)
-        # empty tuples are never terms
-        @test !(() isa TermOrTerms)
-        @test (a, ) isa TermOrTerms
-        @test (a, b) isa TermOrTerms
-        @test (a, b, a&b) isa TermOrTerms
-        @test !(((), a) isa TermOrTerms)
-        # can't contain further tuples
-        @test !((a, (a,), b) isa TermOrTerms)
-
-        # a tuple of AbstractTerms OR Tuples of one or more terms
-        # empty tuples are never terms
-        @test !(() isa TupleTerm)
-        @test (a, ) isa TupleTerm
-        @test (a, b) isa TupleTerm
-        @test (a, b, a&b) isa TupleTerm
-        @test !(((), a) isa TupleTerm)
-        @test (((a,), a) isa TupleTerm)
-
-        # no methods for operators on term and empty tuple (=no type piracy)
+        # TermOrTerms - a term or a vector of terms
+        @test a isa TermOrTerms
+        @test [a] isa TermOrTerms
+        @test [a, b] isa TermOrTerms
+        @test AbstractTerm[a, b, a&b] isa TermOrTerms
+        @test !([1, 2] isa TermOrTerms)
+
+        # no methods for operators on term and tuples (=no type piracy)
         @test_throws MethodError a + ()
         @test_throws MethodError () + a
         @test_throws MethodError a & ()
         @test_throws MethodError () & a
         @test_throws MethodError a ~ ()
         @test_throws MethodError () ~ a
+        @test_throws MethodError a + (a, b)
+        @test_throws MethodError (a, b) + a
 
-        # show methods of empty tuples preserved
+        # show methods
         @test "$(())" == "()"
-        @test "$((a,b))" == "a + b"
-        @test "$((a, ()))" == "(a, ())"
+        @test "$([a, b])" == "a + b"
+        @test "$(AbstractTerm[a, a & b])" == "a + a & b"
+    end
+
+    @testset "vector of terms" begin
+        using StatsModels: terms, termvars, has_schema, drop_term, cleanup,
+            collect_matrix_terms, hasintercept, omitsintercept, termnames,
+            InterceptTerm, ContinuousTerm, CategoricalTerm, MatrixTerm,
+            InteractionTerm
+        a, b, c = term(:a), term(:b), term(:c)
+        one = term(1)
+
+        @testset "+ and & on vectors" begin
+            # all four (term, vector) combinations, with and without duplicates
+            @test [a, b] + c == [a, b, c]
+            @test [a, b] + a == [a, b]
+            @test c + [a, b] == [c, a, b]
+            @test a + [a, b] == [a, b]
+            @test [a, b] + [b, c] == [a, b, c]
+            @test [a, b] + [a, b] == [a, b]
+            # concretely typed vectors hit the Vector{<:AbstractTerm} method,
+            # not elementwise + from Base
+            @test Term[a, b] + Term[b, c] == [a, b, c]
+            @test Term[a, b] + Term[b, c] isa Vector{AbstractTerm}
+            @test view([a, b, c], 1:2) + c == [a, b, c]
+            # a + a is a one-element vector, not a bare term
+            @test a + a == [a]
+            @test a + a isa Vector{AbstractTerm}
+            @test [a] + a == [a]
+
+            # distributive rule
+            @test a & [b, c] == [a & b, a & c]
+            @test [a, b] & c == [a & c, b & c]
+            @test [a, b] & [c, one] == [a & c, a, b & c, b]
+            @test (a & [b, c]) isa Vector{AbstractTerm}
+
+            # * expands with vectors on either side
+            @test a * [b, c] == [a, b, c, a & b, a & c]
+            @test [a, b] * c == [a, b, c, a & c, b & c]
+        end
+
+        @testset "~ always makes a vector rhs" begin
+            f = a ~ b
+            @test f.rhs == [b]
+            @test f.rhs isa Vector{AbstractTerm}
+            @test f.lhs == a
+            # duplicates dropped and sorted by degree, without touching the input
+            rhs = AbstractTerm[a & b, c, a, c, one]
+            f = a ~ rhs
+            @test f.rhs == [one, c, a, a & b]
+            @test rhs == AbstractTerm[a & b, c, a, c, one]
+            @test cleanup(a) == [a]
+            @test ([a, b] ~ c) == FormulaTerm([a, b], [c])
+        end
+
+        @testset "traversal" begin
+            v = AbstractTerm[one, a, b & c]
+            @test terms(v) == [one, a, b, c]
+            @test terms(v) isa Vector{AbstractTerm}
+            @test terms(AbstractTerm[]) == AbstractTerm[]
+            @test terms(MatrixTerm(v)) == terms(v)
+            @test termvars(v) == [:a, :b, :c]
+            @test termvars(AbstractTerm[]) == Symbol[]
+            @test termvars(MatrixTerm(v)) == [:a, :b, :c]
+            @test hasintercept(v) && !omitsintercept(v)
+            @test !hasintercept([a, b]) && omitsintercept([term(0), a])
+            @test termnames(v) == ["1", "a", "b & c"]
+            @test termnames([a]) == ["a"]
+            @test drop_term(v, b & c) == [one, a]
+            @test drop_term(v, b & c) isa Vector{AbstractTerm}
+            @test drop_term(v, term(:z)) == v
+        end
+
+        @testset "schema" begin
+            d = (a=rand(10), b=rand(10), c=repeat(["u", "v"], 5))
+            sch = schema(d)
+            @test !has_schema([a, b])
+            @test !has_schema(AbstractTerm[apply_schema(a, sch), b])
+            @test has_schema(apply_schema([a, b], sch))
+
+            # schema and apply_schema accept a vector, a lone term, and a formula
+            @test schema([a, b], d).schema == schema(a + b, d).schema
+            @test keys(schema(a, d).schema) == Set([a])
+            @test keys(schema(AbstractTerm[one, a, b & c], d).schema) == Set([a, b, c])
+
+            ts = apply_schema(AbstractTerm[a, b, c, a & c], sch)
+            @test ts isa Vector{AbstractTerm}
+            @test ts[1] isa ContinuousTerm && ts[3] isa CategoricalTerm
+            @test ts[4] isa InteractionTerm
+            @test ts[4].terms == [ts[1], ts[3]]
+            # a lone term in a vector stays a vector; duplicates are dropped
+            @test apply_schema([a], sch) == [apply_schema(a, sch)]
+            @test apply_schema([a, a], sch) == [apply_schema(a, sch)]
+            @test apply_schema(AbstractTerm[], sch) == AbstractTerm[]
+            # the formula rhs collapses to a MatrixTerm
+            f = apply_schema(term(:a) ~ b + c, sch)
+            @test f.rhs isa MatrixTerm
+            @test f.rhs == MatrixTerm(apply_schema([b, c], sch))
+            @test width(f.rhs) == 2
+        end
+
+        @testset "modelcols, coefnames, modelmatrix" begin
+            d = (a=collect(1.0:5.0), b=collect(6.0:10.0), c=repeat(["u", "v", "u", "v", "u"]))
+            sch = schema(d)
+            ts = apply_schema([a, b, c], sch)
+            cols = modelcols(ts, d)
+            @test cols isa Vector
+            @test cols[1] == d.a && cols[2] == d.b
+            @test size(cols[3]) == (5, 1)
+            @test coefnames(ts) == ["a", "b", "c: v"]
+            @test modelmatrix(ts, d) == modelmatrix(MatrixTerm(ts), d)
+            @test modelmatrix(a + b, d) == [d.a d.b]
+            @test modelmatrix(a, d) == reshape(d.a, :, 1)
+            @test width(ts[3] & ts[1]) == 1
+        end
+
+        @testset "MatrixTerm and collect_matrix_terms" begin
+            @test MatrixTerm(a) == MatrixTerm([a])
+            @test MatrixTerm((a, b)) == MatrixTerm([a, b])
+            @test MatrixTerm(a).terms isa Vector{AbstractTerm}
+            @test collect_matrix_terms(a) == MatrixTerm(a)
+            @test collect_matrix_terms(MatrixTerm(a)) == MatrixTerm(a)
+            @test collect_matrix_terms([a, b]) == MatrixTerm([a, b])
+            @test collect_matrix_terms(Term[a, b]) == MatrixTerm([a, b])
+            @test InteractionTerm((a, b)) == InteractionTerm([a, b])
+            @test InteractionTerm([a, b]).terms isa Vector{AbstractTerm}
+        end
+
+        @testset "== and hash" begin
+            @test hash(MatrixTerm([a, b])) == hash(MatrixTerm((a, b)))
+            @test hash(a & b) == hash(InteractionTerm([a, b]))
+            @test hash(a ~ b) == hash(a ~ [b])
+            @test MatrixTerm([a, b]) != MatrixTerm([b, a])
+            @test a & b != b & a
+            @test (a ~ b) != (b ~ a)
+            @test length(Set([a ~ b + c, a ~ b + c])) == 1
+            @test length(Set([MatrixTerm([a, b]), MatrixTerm([a, b])])) == 1
+            @test length(unique([a & b, a & b, b & a])) == 2
+        end
+
+        @testset "empty vector of terms" begin
+            e = AbstractTerm[]
+            d = (y=rand(5), a=rand(5))
+            @test e + a == [a] && a + e == [a] && e + e == e
+            @test a & e == e && e & a == e && e & e == e
+            @test a * e == [a]
+            @test terms(e) == e && termvars(e) == Symbol[]
+            @test has_schema(e) && !hasintercept(e) && !omitsintercept(e)
+            @test termnames(e) == String[] && coefnames(e) == String[]
+            @test drop_term(e, a) == e
+            @test isempty(schema(e, d).schema)
+            @test apply_schema(e, schema(d)) == e
+            @test collect_matrix_terms(e) == MatrixTerm(e)
+            @test width(MatrixTerm(e)) == 0
+            @test coefnames(MatrixTerm(e)) == String[]
+            @test termnames(MatrixTerm(e)) == String[]
+            @test modelcols(e, d) == []
+            # an empty rhs behaves like `y ~ 0`: a matrix with n rows and no columns
+            f = apply_schema(term(:y) ~ e, schema(d))
+            @test f.rhs == MatrixTerm(e)
+            y, X = modelcols(f, d)
+            @test y == d.y && size(X) == (5, 0)
+            @test size(modelmatrix(e, d)) == (5, 0)
+            @test modelcols(MatrixTerm(e), (y=1.0, a=2.0)) == Float64[]
+            @test coefnames(f) == ("y", String[])
+            @test mimestring(term(:y) ~ e) == "FormulaTerm\nResponse:\n  y(unknown)\nPredictors:"
+            @test hash(e) == hash(AbstractTerm[])
+            # an interaction needs at least one term
+            @test_throws ArgumentError InteractionTerm(e)
+            @test_throws ArgumentError InteractionTerm(())
+        end
+
+        @testset "show" begin
+            @test string(AbstractTerm[]) == ""
+            @test string([a]) == "a"
+            @test string(MatrixTerm([a, b])) == "a + b"
+            @test mimestring([a, b]) == "a(unknown)\nb(unknown)"
+            @test mimestring(Term[a]) == "a(unknown)"
+        end
     end
 
     @testset "concrete_term error messages" begin
@@ -277,10 +449,10 @@ StatsModels.apply_schema(mt::MultiTerm, sch::StatsModels.Schema, Mod::Type) =
     @testset "sort by degree in ~" begin
         one, a, b = term.([1, :a, :b])
         for zero_deg in [one, InterceptTerm{true}(), InterceptTerm{false}()]
-            @test a + zero_deg == (a, zero_deg)
+            @test a + zero_deg == [a, zero_deg]
             @test (a ~ a + zero_deg) == (a ~ zero_deg + a)
 
-            @test a & b + zero_deg + a == (a & b, zero_deg, a)
+            @test a & b + zero_deg + a == [a & b, zero_deg, a]
             @test (a ~ a & b + zero_deg + a) == (a ~ zero_deg + a + a & b)
         end
     end