diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
new file mode 100644
index 0000000..5e89533
--- /dev/null
+++ b/.github/workflows/CI.yml
@@ -0,0 +1,35 @@
+name: CI
+on:
+  push:
+    branches:
+      - main
+      - master
+    paths-ignore:
+      - 'LICENSE.md'
+      - 'README.md'
+  pull_request:
+    branches:
+      - main
+      - master
+    paths-ignore:
+      - 'LICENSE.md'
+      - 'README.md'
+jobs:
+  ci:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        julia-version: [1, 1.6]
+        julia-arch: [x64]
+        os: [ubuntu-18.04, macos-10.15, windows-2019]
+    steps:
+      - uses: actions/checkout@v2
+      - uses: julia-actions/setup-julia@v1
+        with:
+          version: ${{ matrix.julia-version }}
+      - uses: julia-actions/julia-buildpkg@v0.1
+      - uses: julia-actions/julia-runtest@v0.1
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+      - uses: julia-actions/julia-uploadcodecov@v0.1
+        if:  ${{ startsWith(matrix.os, 'Ubuntu') && startsWith(matrix.julia-version, '1.6') }}
diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml
new file mode 100644
index 0000000..2118253
--- /dev/null
+++ b/.github/workflows/CompatHelper.yml
@@ -0,0 +1,16 @@
+name: CompatHelper
+on:
+  schedule:
+    - cron: 43 7 * * *
+  workflow_dispatch:
+jobs:
+  CompatHelper:
+    runs-on: ubuntu-latest
+    steps:
+      - name: Pkg.add("CompatHelper")
+        run: julia -e 'using Pkg; Pkg.add("CompatHelper")'
+      - name: CompatHelper.main()
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+          COMPATHELPER_PRIV: ${{ secrets.DOCUMENTER_KEY }}
+        run: julia -e 'using CompatHelper; CompatHelper.main()'
\ No newline at end of file
diff --git a/.github/workflows/Documenter.yml b/.github/workflows/Documenter.yml
new file mode 100644
index 0000000..fca4809
--- /dev/null
+++ b/.github/workflows/Documenter.yml
@@ -0,0 +1,23 @@
+name: Documenter
+on:
+  push:
+    branches: [main, master]
+    tags: [v*]
+  pull_request:
+    branches:
+      - main
+      - master
+jobs:
+  Documenter:
+    name: Documentation
+    runs-on: ubuntu-20.04
+    steps:
+      - uses: actions/checkout@v2
+      - uses: julia-actions/setup-julia@v1
+        with:
+          version: 1.6
+      - uses: julia-actions/julia-buildpkg@latest
+      - uses: julia-actions/julia-docdeploy@latest
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}
diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml
new file mode 100644
index 0000000..45ec269
--- /dev/null
+++ b/.github/workflows/TagBot.yml
@@ -0,0 +1,15 @@
+name: TagBot
+on:
+  issue_comment:
+    types:
+      - created
+  workflow_dispatch:
+jobs:
+  TagBot:
+    if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot'
+    runs-on: ubuntu-latest
+    steps:
+      - uses: JuliaRegistries/TagBot@v1
+        with:
+          token: ${{ secrets.GITHUB_TOKEN }}
+          ssh: ${{ secrets.DOCUMENTER_KEY }} # see https://juliadocs.github.io/Documenter.jl/stable/man/hosting/#GitHub-Actions
diff --git a/.travis.yml b/.travis.yml
deleted file mode 100644
index 18d05f6..0000000
--- a/.travis.yml
+++ /dev/null
@@ -1,25 +0,0 @@
-# Documentation: http://docs.travis-ci.com/user/languages/julia/
-language: julia
-os:
-  - linux
-  - osx
-julia:
-  - 1.0
-  - nightly
-notifications:
-  email: false
-after_success:
-  - julia -e 'using Pkg; Pkg.add("Coverage"); using Coverage; Codecov.submit(process_folder())'
-jobs:
-  allow_failures:
-    - julia: nightly
-  fast_finish: true
-  include:
-    - stage: Documentation
-      julia: 1.0
-      script: julia --project=docs -e '
-          using Pkg;
-          Pkg.develop(PackageSpec(path=pwd()));
-          Pkg.instantiate();
-          include("docs/make.jl");'
-      after_success: skip
diff --git a/Project.toml b/Project.toml
index ed7270f..80bff94 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "RegressionFormulae"
 uuid = "545c379f-4ec2-4339-9aea-38f2fb6a8ba2"
-authors = ["Dave Kleinschmidt"]
+authors = ["Dave Kleinschmidt", "Phillip Alday"]
 version = "0.1.0"
 
 [deps]
@@ -8,11 +8,14 @@ Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
 
 [compat]
+StatsBase = "0.33"
 StatsModels = "0.6.7"
-julia = "1"
+julia = "1.6"
 
 [extras]
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test"]
+test = ["StatsBase", "StatsModels", "Test"]
diff --git a/README.md b/README.md
index 77a0a5a..5309e17 100644
--- a/README.md
+++ b/README.md
@@ -4,19 +4,22 @@
 
 [![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://kleinschmidt.github.io/RegressionFormulae.jl/stable)
 [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://kleinschmidt.github.io/RegressionFormulae.jl/dev)
-[![Build Status](https://travis-ci.com/kleinschmidt/RegressionFormulae.jl.svg?branch=master)](https://travis-ci.com/kleinschmidt/RegressionFormulae.jl)
 [![Codecov](https://codecov.io/gh/kleinschmidt/RegressionFormulae.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/kleinschmidt/RegressionFormulae.jl)
 
 Extended [StatsModels.jl
 `@formula`](https://www.github.com/JuliaStats/StatsModels.jl) syntax for
 regression modeling.
 
+Note that the functionality in this package is very new: please verify that the resulting schematized formulae and model coefficient (names) are what you were expecting, especially if you are combining multiple "advanced" formula features.
+
+<!--
 ## Examples
 
 ```julia
 using RegressionFormulae, StatsModels, GLM, DataFrames
 
 ```
+-->
 
 ## Supported syntax ##
 
@@ -29,6 +32,8 @@ using RegressionFormulae, StatsModels, GLM, DataFrames
 Generate all main effects and interactions up to the specified order.  For
 instance, `(a+b+c)^2` generates `a + b + c + a&b + a&c + b&c`, but not `a&b&c`.
 
+**NB:** The presence of interaction terms within the base will result in redundant terms and is currently unsupported.
+
 ## Approach
 
 Extended syntax is supported at two levels.  First, RegressionFormulae.jl
diff --git a/docs/Project.toml b/docs/Project.toml
index dfa65cd..85474d9 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,2 +1,6 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+RegressionFormulae = "545c379f-4ec2-4339-9aea-38f2fb6a8ba2"
+
+[compat]
+Documenter = "0.27"
diff --git a/docs/make.jl b/docs/make.jl
index 4467d33..b8be06a 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,17 +1,16 @@
-using Documenter, RegressionFormula
+using Documenter, RegressionFormulae
 
 makedocs(;
-    modules=[RegressionFormula],
-    format=Documenter.HTML(),
     pages=[
         "Home" => "index.md",
     ],
-    repo="https://github.com/kleinschmidt/RegressionFormula.jl/blob/{commit}{path}#L{line}",
-    sitename="RegressionFormula.jl",
-    authors="Dave Kleinschmidt",
-    assets=String[],
+    repo="https://github.com/kleinschmidt/RegressionFormulae.jl/blob/{commit}{path}#L{line}",
+    sitename="RegressionFormulae.jl",
+    authors="Dave Kleinschmidt and Phillip Alday",
 )
 
 deploydocs(;
-    repo="github.com/kleinschmidt/RegressionFormula.jl",
+    repo="github.com/kleinschmidt/RegressionFormulae.jl",
+    devbranch = "main",
+    push_preview = true
 )
diff --git a/docs/src/index.md b/docs/src/index.md
index c37af91..37775d1 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -1,8 +1,8 @@
-# RegressionFormula.jl
+# RegressionFormulae.jl
 
 ```@index
 ```
 
 ```@autodocs
-Modules = [RegressionFormula]
+Modules = [RegressionFormulae]
 ```
diff --git a/src/RegressionFormulae.jl b/src/RegressionFormulae.jl
index c45ca16..03725b2 100644
--- a/src/RegressionFormulae.jl
+++ b/src/RegressionFormulae.jl
@@ -4,7 +4,7 @@ using StatsModels
 using Combinatorics
 using Base.Iterators
 
-import StatsModels: apply_schema
+using StatsModels: apply_schema
 
 const TermTuple = NTuple{N, AbstractTerm} where N
 const Schemas = Union{StatsModels.Schema, StatsModels.FullRank}
diff --git a/src/fulldummy.jl b/src/fulldummy.jl
index 157e06d..6a1d7df 100644
--- a/src/fulldummy.jl
+++ b/src/fulldummy.jl
@@ -18,6 +18,3 @@ function fulldummy(t::CategoricalTerm)
     )
     t = CategoricalTerm(t.sym, new_contrasts)
 end
-
-fulldummy(x) =
-    throw(ArgumentError("fulldummy isn't supported outside of a MixedModel formula"))
diff --git a/src/nesting.jl b/src/nesting.jl
index c74fc34..8d36b67 100644
--- a/src/nesting.jl
+++ b/src/nesting.jl
@@ -1,5 +1,8 @@
-# TODO: handle nested grouping.  This parses as (a / b) / c instead of
-# (/)(a,b,c) so need to do the reduction manually
+
+_isfulldummy(x::AbstractTerm) = false
+function _isfulldummy(x::CategoricalTerm)
+    return isa(x.contrasts, StatsModels.ContrastsMatrix{StatsModels.FullDummyCoding})
+end
 
 """
     group / term
@@ -7,10 +10,26 @@
 Generate predictors for `term` within each level of `group`.  Implemented as
 `group + fulldummy(group) & term`.
 """
-function Base.:(/)(args::AbstractTerm...)
-    groups = (&)(args[1:end-1]...)
-    term = last(args)
-    return groups + fulldummy(groups) & term
+function Base.:(/)(outer::CategoricalTerm, inner::AbstractTerm)
+    return outer + fulldummy(outer) & inner
+end
+
+function Base.:(/)(outer::TermTuple, inner::AbstractTerm)
+    return outer[1:end-1] + last(outer) / inner
+end
+
+function Base.:(/)(outer::InteractionTerm, inner::AbstractTerm)
+    # we should only get here via expansion where the interaction term,
+    # but who knows what devious things users will try
+    all(_isfulldummy, outer.terms[1:end-1]) ||
+        throw(ArgumentError("Outer interactions in a nesting must consist only " *
+                            " of categorical terms with FullDummyCoding, got $outer"))
+    return outer + outer & inner
+end
+
+function Base.:(/)(outer::AbstractTerm, inner::AbstractTerm)
+    throw(ArgumentError("nesting terms requires categorical grouping term, got $outer / $inner " *
+                        "Manually specify $outer as `CategoricalTerm` in hints/contrasts"))
 end
 
 function StatsModels.apply_schema(
@@ -19,15 +38,9 @@ function StatsModels.apply_schema(
     Mod::Type{<:RegressionModel},
 )
     length(t.args_parsed) == 2 ||
-        throw(ArgumentError("malformed nesting term: $t (Exactly two arguments required"))
+        throw(ArgumentError("malformed nesting term: $t (Exactly two arguments required)"))
 
     args = apply_schema.(t.args_parsed, Ref(sch), Mod)
 
-    map(args[1:end-1]) do arg 
-        typeof(arg) <: CategoricalTerm ||
-            throw(ArgumentError("nesting terms requires categorical grouping term, got $arg.  "*
-                                "Manually specify $first as `CategoricalTerm` in hints/contrasts"))
-    end
-    
-    return (/)(args...)
+    return first(args) / last(args)
 end
diff --git a/src/power.jl b/src/power.jl
index 1ea35b6..fb5cd19 100644
--- a/src/power.jl
+++ b/src/power.jl
@@ -4,11 +4,20 @@ combinations_upto(x, n) = Iterators.flatten(combinations(x, i) for i in 1:n)
     (term1, term2, ...) ^ n
 
 Generate all interactions of terms up to order ``n``.
+
+!!! warning
+    If any term is an `InteractionTerm`, then nonsensical interactions may
+    arise, e.g. `a & a & b`.
 """
-function Base.:(^)(args::TermTuple, deg::ConstantTerm)
+function Base.:(^)(args::TermTuple, deg::ConstantTerm{<:Integer})
+    deg.n > 0 || throw(ArgumentError("power should be greater than zero (got $deg)"))
     tuple(((&)(terms...) for terms in combinations_upto(args, deg.n))...)
 end
 
+function Base.:(^)(::TermTuple, deg::AbstractTerm)
+    throw(ArgumentError("power should be an integer constant (got $deg)"))
+end
+
 function StatsModels.apply_schema(
     t::FunctionTerm{typeof(^)},
     sch::StatsModels.FullRank,
@@ -17,7 +26,7 @@ function StatsModels.apply_schema(
     length(t.args_parsed) == 2 ||
         throw(ArgumentError("invalid term $t: should have exactly two arguments"))
     first, second = t.args_parsed
-    second isa ConstantTerm ||
-        throw(ArgumentError("invalid term $t: power should be a number (got $second)"))
+    second isa ConstantTerm{<:Integer} ||
+        throw(ArgumentError("invalid term $t: power should be an integer (got $second)"))
     apply_schema.(first^second, Ref(sch), ctx)
 end
diff --git a/test/dummymod.jl b/test/dummymod.jl
new file mode 100644
index 0000000..ac2fb35
--- /dev/null
+++ b/test/dummymod.jl
@@ -0,0 +1,103 @@
+isdefined(@__MODULE__, :DUMMY_MOD) || (const DUMMY_MOD = true)
+
+# taken from StatsModels
+# https://github.com/JuliaStats/StatsModels.jl/blob/dee41c287033c0e9c18714a8d63bae61302027a6/test/statsmodel.jl
+
+using StatsBase
+
+# A dummy RegressionModel type
+struct DummyMod <: RegressionModel
+    beta::Vector{Float64}
+    x::Matrix
+    y::Vector
+end
+
+## dumb fit method: just copy the x and y input over
+StatsBase.fit(::Type{DummyMod}, x::Matrix, y::Vector) =
+    DummyMod(collect(1:size(x, 2)), x, y)
+StatsBase.response(mod::DummyMod) = mod.y
+## dumb coeftable: just prints the "beta" values
+StatsBase.coeftable(mod::DummyMod) =
+    CoefTable(reshape(mod.beta, (size(mod.beta,1), 1)),
+              ["'beta' value"],
+              ["" for n in 1:size(mod.x,2)],
+              0)
+# dumb predict: return values predicted by "beta" and dummy confidence bounds
+function StatsBase.predict(mod::DummyMod;
+                           interval::Union{Nothing,Symbol}=nothing)
+    pred = mod.x * mod.beta
+    if interval === nothing
+        return pred
+    elseif interval === :prediction
+        return (prediction=pred, lower=pred .- 1, upper=pred .+ 1)
+    else
+        throw(ArgumentError("value not allowed for interval"))
+    end
+end
+function StatsBase.predict(mod::DummyMod, newX::Matrix;
+                           interval::Union{Nothing,Symbol}=nothing)
+    pred = newX * mod.beta
+    if interval === nothing
+        return pred
+    elseif interval === :prediction
+        return (prediction=pred, lower=pred .- 1, upper=pred .+ 1)
+    else
+        throw(ArgumentError("value not allowed for interval"))
+    end
+end
+StatsBase.dof(mod::DummyMod) = length(mod.beta)
+StatsBase.dof_residual(mod::DummyMod) = length(mod.y) - length(mod.beta)
+StatsBase.nobs(mod::DummyMod) = length(mod.y)
+StatsBase.deviance(mod::DummyMod) = sum((response(mod) .- predict(mod)).^2)
+# Incorrect but simple definition
+StatsModels.isnested(mod1::DummyMod, mod2::DummyMod; atol::Real=0.0) =
+    dof(mod1) <= dof(mod2)
+StatsBase.loglikelihood(mod::DummyMod) = -sum((response(mod) .- predict(mod)).^2)
+StatsBase.loglikelihood(mod::DummyMod, ::Colon) = -(response(mod) .- predict(mod)).^2
+
+# A dummy RegressionModel type that does not support intercept
+struct DummyModNoIntercept <: RegressionModel
+    beta::Vector{Float64}
+    x::Matrix
+    y::Vector
+end
+
+StatsModels.drop_intercept(::Type{DummyModNoIntercept}) = true
+
+## dumb fit method: just copy the x and y input over
+StatsBase.fit(::Type{DummyModNoIntercept}, x::Matrix, y::Vector) =
+    DummyModNoIntercept(collect(1:size(x, 2)), x, y)
+StatsBase.response(mod::DummyModNoIntercept) = mod.y
+## dumb coeftable: just prints the "beta" values
+StatsBase.coeftable(mod::DummyModNoIntercept) =
+    CoefTable(reshape(mod.beta, (size(mod.beta,1), 1)),
+              ["'beta' value"],
+              ["" for n in 1:size(mod.x,2)],
+              0)
+# dumb predict: return values predicted by "beta" and dummy confidence bounds
+function StatsBase.predict(mod::DummyModNoIntercept;
+                           interval::Union{Nothing,Symbol}=nothing)
+    pred = mod.x * mod.beta
+    if interval === nothing
+        return pred
+    elseif interval === :prediction
+        return (prediction=pred, lower=pred .- 1, upper=pred .+ 1)
+    else
+        throw(ArgumentError("value not allowed for interval"))
+    end
+end
+function StatsBase.predict(mod::DummyModNoIntercept, newX::Matrix;
+                           interval::Union{Nothing,Symbol}=nothing)
+    pred = newX * mod.beta
+    if interval === nothing
+        return pred
+    elseif interval === :prediction
+        return (prediction=pred, lower=pred .- 1, upper=pred .+ 1)
+    else
+        throw(ArgumentError("value not allowed for interval"))
+    end
+end
+StatsBase.dof(mod::DummyModNoIntercept) = length(mod.beta)
+StatsBase.dof_residual(mod::DummyModNoIntercept) = length(mod.y) - length(mod.beta)
+StatsBase.nobs(mod::DummyModNoIntercept) = length(mod.y)
+StatsBase.deviance(mod::DummyModNoIntercept) = sum((response(mod) .- predict(mod)).^2)
diff --git a/test/fulldummy.jl b/test/fulldummy.jl
new file mode 100644
index 0000000..627d5f3
--- /dev/null
+++ b/test/fulldummy.jl
@@ -0,0 +1,15 @@
+using StatsModels
+using RegressionFormulae
+using RegressionFormulae: fulldummy
+using Test
+
+include("dummymod.jl")
+
+dat = (; y=zeros(10), a=["u","i","o"], b=["q","w","e"], c=["s","d","f"], x=1:10)
+
+@testset "error checking" begin
+    @test_throws ArgumentError fulldummy(term(:a))
+
+    sch = schema(dat)
+    @test_throws ArgumentError apply_schema(@formula(y ~ fulldummy(x)), sch, RegressionModel)
+end
diff --git a/test/nesting.jl b/test/nesting.jl
new file mode 100644
index 0000000..dea2162
--- /dev/null
+++ b/test/nesting.jl
@@ -0,0 +1,71 @@
+using StatsModels
+using RegressionFormulae
+using Test
+
+include("dummymod.jl")
+
+dat = (; y=zeros(3), a=["u","i","o"], b=["q","w","e"], c=["s","d","f"], x=1:3)
+
+@testset "error checking" begin
+    sch = schema(dat)
+    @test_throws ArgumentError apply_schema(@formula(y ~ x / a), sch, RegressionModel)
+    @test_throws ArgumentError apply_schema(@formula(y ~ /(a, b, c)), sch, RegressionModel)
+    @test !RegressionFormulae._isfulldummy(term(:a))
+end
+
+@testset "single nesting level" begin
+    m = fit(DummyMod, @formula(y ~ 0 + a / b), dat)
+    @test coefnames(m) == ["a: i", "a: o", "a: u",
+                           "a: i & b: q", "a: o & b: q", "a: u & b: q",
+                           "a: i & b: w", "a: o & b: w", "a: u & b: w"]
+
+    m = fit(DummyMod, @formula(y ~ 0 + a / x), dat)
+    @test coefnames(m) == ["a: i", "a: o", "a: u",
+                           "a: i & x", "a: o & x", "a: u & x"]
+
+    m = fit(DummyMod, @formula(y ~ 1 + a / b), dat)
+    @test coefnames(m) == ["(Intercept)", "a: o", "a: u",
+                           "a: i & b: q", "a: o & b: q", "a: u & b: q",
+                           "a: i & b: w", "a: o & b: w", "a: u & b: w"]
+
+    m = fit(DummyMod, @formula(y ~ 1 + a / x), dat)
+    @test coefnames(m) == ["(Intercept)", "a: o", "a: u",
+                           "a: i & x", "a: o & x", "a: u & x"]
+end
+
+@testset "multiple nesting levels" begin
+    m = fit(DummyMod, @formula(y ~ 0 + a / b / c), dat)
+    @test coefnames(m) == ["a: i", "a: o", "a: u",
+                           "a: i & b: q", "a: o & b: q", "a: u & b: q",
+                           "a: i & b: w", "a: o & b: w", "a: u & b: w",
+                           "a: i & b: q & c: f", "a: o & b: q & c: f",
+                           "a: u & b: q & c: f", "a: i & b: w & c: f",
+                           "a: o & b: w & c: f", "a: u & b: w & c: f",
+                           "a: i & b: q & c: s", "a: o & b: q & c: s",
+                           "a: u & b: q & c: s", "a: i & b: w & c: s",
+                           "a: o & b: w & c: s", "a: u & b: w & c: s"]
+    m = fit(DummyMod, @formula(y ~ 0 + a / b / x), dat)
+    @test coefnames(m) == ["a: i", "a: o", "a: u",
+                           "a: i & b: q", "a: o & b: q", "a: u & b: q",
+                           "a: i & b: w", "a: o & b: w", "a: u & b: w",
+                           "a: i & b: q & x", "a: o & b: q & x",
+                           "a: u & b: q & x", "a: i & b: w & x",
+                           "a: o & b: w & x", "a: u & b: w & x"]
+    m = fit(DummyMod, @formula(y ~ 1 + a / b / c), dat)
+    @test coefnames(m) == ["(Intercept)", "a: o", "a: u",
+                           "a: i & b: q", "a: o & b: q", "a: u & b: q",
+                           "a: i & b: w", "a: o & b: w", "a: u & b: w",
+                           "a: i & b: q & c: f", "a: o & b: q & c: f",
+                           "a: u & b: q & c: f", "a: i & b: w & c: f",
+                           "a: o & b: w & c: f", "a: u & b: w & c: f",
+                           "a: i & b: q & c: s", "a: o & b: q & c: s",
+                           "a: u & b: q & c: s", "a: i & b: w & c: s",
+                           "a: o & b: w & c: s", "a: u & b: w & c: s"]
+    m = fit(DummyMod, @formula(y ~ 1 + a / b / x), dat)
+    @test coefnames(m) == ["(Intercept)", "a: o", "a: u",
+                           "a: i & b: q", "a: o & b: q", "a: u & b: q",
+                           "a: i & b: w", "a: o & b: w", "a: u & b: w",
+                           "a: i & b: q & x", "a: o & b: q & x",
+                           "a: u & b: q & x", "a: i & b: w & x",
+                           "a: o & b: w & x", "a: u & b: w & x"]
+end
diff --git a/test/power.jl b/test/power.jl
new file mode 100644
index 0000000..e0c3754
--- /dev/null
+++ b/test/power.jl
@@ -0,0 +1,41 @@
+using StatsModels
+using RegressionFormulae
+using Test
+
+include("dummymod.jl")
+
+dat = (; y=zeros(3), a=1:3, b=11:13, c=21:23, d=31:33, e=["u", "i", "o"])
+
+
+@testset "error checking" begin
+    @test_throws ArgumentError (term(:b),) ^ term(:a)
+    @test_throws ArgumentError (term(:b),) ^ term(2.5)
+    @test_throws ArgumentError (term(:b),) ^ term(-2)
+
+    sch = schema(dat)
+    @test_throws ArgumentError apply_schema(@formula(y ~ ^(a, b, c)), sch, RegressionModel)
+    @test_throws ArgumentError apply_schema(@formula(y ~ (a+b)^2.5), sch, RegressionModel)
+end
+
+@testset "powers of sums" begin
+    m = fit(DummyMod, @formula(y ~ (a + b + c + d)^3), dat)
+    @test coefnames(m) == ["(Intercept)", "a", "b", "c", "d",
+                           "a & b", "a & c", "a & d", "b & c", "b & d", "c & d",
+                           "a & b & c", "a & b & d", "a & c & d", "b & c & d"]
+
+    m = fit(DummyMod, @formula(y ~ (1 + a + b + c + d)^3), dat)
+    @test_broken coefnames(m) == ["(Intercept)", "a", "b", "c", "d",
+                                  "a & b", "a & c", "a & d", "b & c", "b & d", "c & d",
+                                  "a & b & c", "a & b & d", "a & c & d", "b & c & d"]
+
+    m = fit(DummyMod, @formula(y ~ (a + b + e)^2), dat)
+    @test coefnames(m) == ["(Intercept)", "a", "b", "e: o", "e: u",
+                           "a & b", "a & e: o", "a & e: u", "b & e: o", "b & e: u"]
+end
+
+@testset "embedded interactions" begin
+    m = fit(DummyMod, @formula(y ~ (a + b + c * d)^3), dat)
+    cn = coefnames(m)
+    @test_broken !("a & c & c & d" in cn)
+    @test_broken cn == unique(cn)
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 10c97a4..0345cdb 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,6 +1,15 @@
 using RegressionFormulae
+using StatsModels
 using Test
 
-@testset "RegressionFormulae.jl" begin
-    # Write your own tests here.
+@testset "fulldummy" begin
+    include("fulldummy.jl")
+end
+
+@testset "nesting" begin
+    include("nesting.jl")
+end
+
+@testset "powers of terms" begin
+    include("power.jl")
 end