diff --git a/src/matrix/kernelkroneckermat.jl b/src/matrix/kernelkroneckermat.jl
index 113d0f53d..90ca1d2da 100644
--- a/src/matrix/kernelkroneckermat.jl
+++ b/src/matrix/kernelkroneckermat.jl
@@ -82,10 +82,10 @@ efficiently inverted or decomposed. See also [`kernelmatrix`](@ref).
 function kronecker_kernelmatrix(
     k::Union{IndependentMOKernel,IntrinsicCoregionMOKernel}, x::MOI, y::MOI
 ) where {MOI<:IsotopicMOInputsUnion}
-    x.out_dim == y.out_dim ||
-        throw(DimensionMismatch("`x` and `y` must have the same `out_dim`"))
+    x.out_axis == y.out_axis ||
+        throw(DimensionMismatch("`x` and `y` must have the same `out_axis`"))
     Kfeatures = kernelmatrix(k.kernel, x.x, y.x)
-    Koutputs = _mo_output_covariance(k, x.out_dim)
+    Koutputs = _mo_output_covariance(k, length(x.out_axis))
     return _kernelmatrix_kroneckerjl_helper(MOI, Kfeatures, Koutputs)
 end
 
@@ -93,7 +93,7 @@ function kronecker_kernelmatrix(
     k::Union{IndependentMOKernel,IntrinsicCoregionMOKernel}, x::MOI
 ) where {MOI<:IsotopicMOInputsUnion}
     Kfeatures = kernelmatrix(k.kernel, x.x)
-    Koutputs = _mo_output_covariance(k, x.out_dim)
+    Koutputs = _mo_output_covariance(k, length(x.out_axis))
     return _kernelmatrix_kroneckerjl_helper(MOI, Kfeatures, Koutputs)
 end
 
diff --git a/src/mokernels/independent.jl b/src/mokernels/independent.jl
index 1f7811b14..a6e5e2811 100644
--- a/src/mokernels/independent.jl
+++ b/src/mokernels/independent.jl
@@ -27,15 +27,15 @@ function (κ::IndependentMOKernel)((x, px)::Tuple{Any,Int}, (y, py)::Tuple{Any,I
     return κ.kernel(x, y) * (px == py)
 end
 
-_mo_output_covariance(k::IndependentMOKernel, out_dim) = Eye{Bool}(out_dim)
+_mo_output_covariance(k::IndependentMOKernel, out_dim::Integer) = Eye{Bool}(out_dim)
 
 function kernelmatrix(
     k::IndependentMOKernel, x::MOI, y::MOI
 ) where {MOI<:IsotopicMOInputsUnion}
-    x.out_dim == y.out_dim ||
-        throw(DimensionMismatch("`x` and `y` must have the same `out_dim`"))
+    x.out_axis == y.out_axis ||
+        throw(DimensionMismatch("`x` and `y` must have the same `out_axis`"))
     Kfeatures = kernelmatrix(k.kernel, x.x, y.x)
-    Koutputs = _mo_output_covariance(k, x.out_dim)
+    Koutputs = _mo_output_covariance(k, length(x.out_axis))
     return _kernelmatrix_kron_helper(MOI, Kfeatures, Koutputs)
 end
 
@@ -43,10 +43,10 @@ if VERSION >= v"1.6"
     function kernelmatrix!(
         K::AbstractMatrix, k::IndependentMOKernel, x::MOI, y::MOI
     ) where {MOI<:IsotopicMOInputsUnion}
-        x.out_dim == y.out_dim ||
-            throw(DimensionMismatch("`x` and `y` must have the same `out_dim`"))
+        x.out_axis == y.out_axis ||
+            throw(DimensionMismatch("`x` and `y` must have the same `out_axis`"))
         Kfeatures = kernelmatrix(k.kernel, x.x, y.x)
-        Koutputs = _mo_output_covariance(k, x.out_dim)
+        Koutputs = _mo_output_covariance(k, length(x.out_axis))
         return _kernelmatrix_kron_helper!(K, MOI, Kfeatures, Koutputs)
     end
 end
diff --git a/src/mokernels/intrinsiccoregion.jl b/src/mokernels/intrinsiccoregion.jl
index 0a940796b..89b81434e 100644
--- a/src/mokernels/intrinsiccoregion.jl
+++ b/src/mokernels/intrinsiccoregion.jl
@@ -42,7 +42,7 @@ function (k::IntrinsicCoregionMOKernel)((x, px)::Tuple{Any,Int}, (y, py)::Tuple{
     return k.B[px, py] * k.kernel(x, y)
 end
 
-function _mo_output_covariance(k::IntrinsicCoregionMOKernel, out_dim)
+function _mo_output_covariance(k::IntrinsicCoregionMOKernel, out_dim::Integer)
     @assert size(k.B) == (out_dim, out_dim)
     return k.B
 end
@@ -50,10 +50,10 @@ end
 function kernelmatrix(
     k::IntrinsicCoregionMOKernel, x::MOI, y::MOI
 ) where {MOI<:IsotopicMOInputsUnion}
-    x.out_dim == y.out_dim ||
-        throw(DimensionMismatch("`x` and `y` must have the same `out_dim`"))
+    x.out_axis == y.out_axis ||
+        throw(DimensionMismatch("`x` and `y` must have the same `out_axis`"))
     Kfeatures = kernelmatrix(k.kernel, x.x, y.x)
-    Koutputs = _mo_output_covariance(k, x.out_dim)
+    Koutputs = _mo_output_covariance(k, length(x.out_axis))
     return _kernelmatrix_kron_helper(MOI, Kfeatures, Koutputs)
 end
 
@@ -61,10 +61,10 @@ if VERSION >= v"1.6"
     function kernelmatrix!(
         K::AbstractMatrix, k::IntrinsicCoregionMOKernel, x::MOI, y::MOI
     ) where {MOI<:IsotopicMOInputsUnion}
-        x.out_dim == y.out_dim ||
-            throw(DimensionMismatch("`x` and `y` must have the same `out_dim`"))
+        x.out_axis == y.out_axis ||
+            throw(DimensionMismatch("`x` and `y` must have the same `out_axis`"))
         Kfeatures = kernelmatrix(k.kernel, x.x, y.x)
-        Koutputs = _mo_output_covariance(k, x.out_dim)
+        Koutputs = _mo_output_covariance(k, length(x.out_axis))
         return _kernelmatrix_kron_helper!(K, MOI, Kfeatures, Koutputs)
     end
 end
diff --git a/src/mokernels/moinput.jl b/src/mokernels/moinput.jl
index 2c2fd1f62..cf81b9f2f 100644
--- a/src/mokernels/moinput.jl
+++ b/src/mokernels/moinput.jl
@@ -7,7 +7,7 @@
 julia> x = [1, 2, 3];
 
 julia> KernelFunctions.MOInputIsotopicByFeatures(x, 2)
-6-element KernelFunctions.MOInputIsotopicByFeatures{Int64, Vector{Int64}, Int64}:
+6-element KernelFunctions.MOInputIsotopicByFeatures{Int64, Int64, Vector{Int64}, Base.OneTo{Int64}}:
  (1, 1)
  (1, 2)
  (2, 1)
@@ -24,10 +24,14 @@ The first `out_dim` elements represent all outputs for the first input, the seco
 
 See [Inputs for Multiple Outputs](@ref) in the docs for more info.
 """
-struct MOInputIsotopicByFeatures{S,T<:AbstractVector{S},Tout_dim<:Integer} <:
-       AbstractVector{Tuple{S,Int}}
+struct MOInputIsotopicByFeatures{S,I,T<:AbstractVector{S},Tout_axis<:AbstractVector{I}} <:
+       AbstractVector{Tuple{S,I}}
     x::T
-    out_dim::Tout_dim
+    out_axis::Tout_axis
+end
+
+function MOInputIsotopicByFeatures(x::AbstractVector, out_dim::Integer)
+    return MOInputIsotopicByFeatures(x, Base.OneTo(out_dim))
 end
 
 """
@@ -39,7 +43,7 @@ end
 julia> x = [1, 2, 3];
 
 julia> KernelFunctions.MOInputIsotopicByOutputs(x, 2)
-6-element KernelFunctions.MOInputIsotopicByOutputs{Int64, Vector{Int64}, Int64}:
+6-element KernelFunctions.MOInputIsotopicByOutputs{Int64, Int64, Vector{Int64}, Base.OneTo{Int64}}:
  (1, 1)
  (2, 1)
  (3, 1)
@@ -54,10 +58,14 @@ As shown above, an `MOInputIsotopicByOutputs` represents a vector of tuples.
 The first `length(x)` elements represent the inputs for the first output, the second
 `length(x)` elements represent the inputs for the second output, etc.
 """
-struct MOInputIsotopicByOutputs{S,T<:AbstractVector{S},Tout_dim<:Integer} <:
-       AbstractVector{Tuple{S,Int}}
+struct MOInputIsotopicByOutputs{S,I,T<:AbstractVector{S},Tout_axis<:AbstractVector{I}} <:
+       AbstractVector{Tuple{S,I}}
     x::T
-    out_dim::Tout_dim
+    out_axis::Tout_axis
+end
+
+function MOInputIsotopicByOutputs(x::AbstractVector, out_dim::Integer)
+    return MOInputIsotopicByOutputs(x, Base.OneTo(out_dim))
 end
 
 const IsotopicMOInputsUnion = Union{MOInputIsotopicByFeatures,MOInputIsotopicByOutputs}
@@ -66,26 +74,30 @@ function Base.getindex(inp::MOInputIsotopicByOutputs, ind::Integer)
     @boundscheck checkbounds(inp, ind)
     output_index, feature_index = fldmod1(ind, length(inp.x))
     feature = @inbounds inp.x[feature_index]
-    return feature, output_index
+    out_indices = axes(inp.out_axis, 1)
+    out_idx = out_indices[begin + output_index - 1]
+    return feature, @inbounds inp.out_axis[out_idx]
 end
 
 function Base.getindex(inp::MOInputIsotopicByFeatures, ind::Integer)
     @boundscheck checkbounds(inp, ind)
-    feature_index, output_index = fldmod1(ind, inp.out_dim)
+    feature_index, output_index = fldmod1(ind, length(inp.out_axis))
     feature = @inbounds inp.x[feature_index]
-    return feature, output_index
+    out_indices = axes(inp.out_axis, 1)
+    out_idx = out_indices[begin + output_index - 1]
+    return feature, @inbounds inp.out_axis[out_idx]
 end
 
-Base.size(inp::IsotopicMOInputsUnion) = (inp.out_dim * length(inp.x),)
+Base.size(inp::IsotopicMOInputsUnion) = (length(inp.out_axis) * length(inp.x),)
 
 function Base.vcat(x::MOInputIsotopicByFeatures, y::MOInputIsotopicByFeatures)
-    x.out_dim == y.out_dim || throw(DimensionMismatch("out_dim mismatch"))
-    return MOInputIsotopicByFeatures(vcat(x.x, y.x), x.out_dim)
+    x.out_axis == y.out_axis || throw(DimensionMismatch("out_axis mismatch"))
+    return MOInputIsotopicByFeatures(vcat(x.x, y.x), x.out_axis)
 end
 
 function Base.vcat(x::MOInputIsotopicByOutputs, y::MOInputIsotopicByOutputs)
-    x.out_dim == y.out_dim || throw(DimensionMismatch("out_dim mismatch"))
-    return MOInputIsotopicByOutputs(vcat(x.x, y.x), x.out_dim)
+    x.out_axis == y.out_axis || throw(DimensionMismatch("out_axis mismatch"))
+    return MOInputIsotopicByOutputs(vcat(x.x, y.x), x.out_axis)
 end
 
 """
@@ -98,7 +110,7 @@ A data type to accommodate modelling multi-dimensional output data.
 julia> x = [1, 2, 3];
 
 julia> MOInput(x, 2)
-6-element KernelFunctions.MOInputIsotopicByOutputs{Int64, Vector{Int64}, Int64}:
+6-element KernelFunctions.MOInputIsotopicByOutputs{Int64, Int64, Vector{Int64}, Base.OneTo{Int64}}:
  (1, 1)
  (2, 1)
  (3, 1)
@@ -138,7 +150,7 @@ julia> Y = [1.1 2.1 3.1; 1.2 2.2 3.2]
 julia> inputs, outputs = prepare_isotopic_multi_output_data(x, ColVecs(Y));
 
 julia> inputs
-6-element KernelFunctions.MOInputIsotopicByFeatures{Float64, Vector{Float64}, Int64}:
+6-element KernelFunctions.MOInputIsotopicByFeatures{Float64, Int64, Vector{Float64}, Base.OneTo{Int64}}:
  (1.0, 1)
  (1.0, 2)
  (2.0, 1)
@@ -186,7 +198,7 @@ julia> Y = [1.1 1.2; 2.1 2.2; 3.1 3.2]
 julia> inputs, outputs = prepare_isotopic_multi_output_data(x, RowVecs(Y));
 
 julia> inputs
-6-element KernelFunctions.MOInputIsotopicByOutputs{Float64, Vector{Float64}, Int64}:
+6-element KernelFunctions.MOInputIsotopicByOutputs{Float64, Int64, Vector{Float64}, Base.OneTo{Int64}}:
  (1.0, 1)
  (2.0, 1)
  (3.0, 1)
diff --git a/src/mokernels/slfm.jl b/src/mokernels/slfm.jl
index 521d07de7..70ab2007f 100644
--- a/src/mokernels/slfm.jl
+++ b/src/mokernels/slfm.jl
@@ -34,8 +34,8 @@ function (κ::LatentFactorMOKernel)((x, px)::Tuple{Any,Int}, (y, py)::Tuple{Any,
 end
 
 function kernelmatrix(k::LatentFactorMOKernel, x::MOInput, y::MOInput)
-    x.out_dim == y.out_dim || error("`x` and `y` should have the same output dimension")
-    x.out_dim == size(k.A, 1) ||
+    x.out_axis == y.out_axis || error("`x` and `y` should have the same output dimension")
+    length(x.out_axis) == size(k.A, 1) ||
         error("Kernel not compatible with the given multi-output inputs")
 
     # Weights matrix ((out_dim x out_dim) x length(k.g))
diff --git a/test/Project.toml b/test/Project.toml
index e16f39a6c..2a54c4331 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -9,6 +9,7 @@ Functors = "d9f16b24-f501-4c13-a1f2-28368ffc5196"
 Kronecker = "2c470bb0-bcc8-11e8-3dad-c9649493f05e"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
+OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
diff --git a/test/mokernels/moinput.jl b/test/mokernels/moinput.jl
index ce1a4757e..e56ef3393 100644
--- a/test/mokernels/moinput.jl
+++ b/test/mokernels/moinput.jl
@@ -4,6 +4,8 @@
     type_2 = AbstractVector{Tuple{AbstractVector{Vector{Float64}},Int}}
 
     @testset "isotopicbyoutputs" begin
+        @test KernelFunctions.MOInputIsotopicByOutputs(x, []) == []
+
         ibo = KernelFunctions.MOInputIsotopicByOutputs(x, 3)
         ibo2 = KernelFunctions.MOInputIsotopicByOutputs(x, 2)
 
@@ -27,9 +29,20 @@
         @test ibo[7] == (x[3], 2)
         @test all([(x_, i) for i in 1:3 for x_ in x] .== ibo)
         @inferred getindex(ibo, 1)
+
+        # test non-standard array
+        ibo3 = KernelFunctions.MOInputIsotopicByOutputs(x, OA.Origin(0)([:a, :b, :c]))
+        ibo4 = KernelFunctions.MOInputIsotopicByOutputs(x, [:a, :b, :c])
+        @test ibo3 == ibo4
+        @test length(ibo3) == 12
+        @test lastindex(ibo3) == 12
+        @test firstindex(ibo3) == 1
+        @test_throws BoundsError ibo3[0]
+        @test_throws BoundsError ibo3[13]
     end
 
     @testset "isotopicbyfeatures" begin
+        @test KernelFunctions.MOInputIsotopicByFeatures([], [1.0, 2.0]) == []
         ibf = KernelFunctions.MOInputIsotopicByFeatures(x, 3)
 
         @test isa(ibf, type_1) == true
@@ -48,6 +61,15 @@
         @test ibf[7] == (x[3], 1)
         @test all([(x_, i) for x_ in x for i in 1:3] .== ibf)
         @inferred getindex(ibf, 1)
+
+        # test non-standard array
+        ibf2 = KernelFunctions.MOInputIsotopicByFeatures(x, OA.Origin(0)(["a", "b", "c"]))
+        ibf3 = KernelFunctions.MOInputIsotopicByFeatures(x, ["a", "b", "c"])
+        @test ibf2 == ibf3
+        @test lastindex(ibf2) == 12
+        @test firstindex(ibf2) == 1
+        @test_throws BoundsError ibf2[0]
+        @test_throws BoundsError ibf2[13]
     end
 
     @testset "prepare_isotopic_multi_output_data" begin
diff --git a/test/runtests.jl b/test/runtests.jl
index caf43cb91..d6ea88b06 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -18,6 +18,7 @@ using ForwardDiff: ForwardDiff
 using ReverseDiff: ReverseDiff
 using FiniteDifferences: FiniteDifferences
 using Compat: only
+import OffsetArrays as OA
 
 using KernelFunctions: SimpleKernel, metric, kappa, ColVecs, RowVecs, TestUtils