diff --git a/src/krig/simple.jl b/src/krig/simple.jl
index 46ad3eb..53fc0df 100644
--- a/src/krig/simple.jl
+++ b/src/krig/simple.jl
@@ -40,10 +40,10 @@ function krigmean(fitted::FittedKriging{<:SimpleKriging}, weights::KrigingWeight
 
   cols = Tables.columns(values(d))
   @inbounds map(1:k) do j
-    sum(1:k) do p
+    μ[j] + sum(1:k) do p
       λₚ = @view λ[p:k:end, j]
       zₚ = Tables.getcolumn(cols, vars[p])
-      μ[p] + sum(i -> λₚ[i] ⦿ (zₚ[i] - μ[p]), eachindex(λₚ, zₚ))
+      sum(i -> λₚ[i] ⦿ (zₚ[i] - μ[p]), eachindex(λₚ, zₚ))
     end
   end
 end
diff --git a/test/krig.jl b/test/krig.jl
index 1dc03da..357e865 100644
--- a/test/krig.jl
+++ b/test/krig.jl
@@ -355,6 +355,13 @@
     @test all(≥(0), var(okdist))
     @test all(≥(0), var(ukdist))
     @test all(≥(0), var(dkdist))
+
+    # simple CoKriging with non-zero mean
+    sk = GeoStatsModels.fit(SK(γ, [0.7, 0.2, 0.1]), data)
+    for i in 1:3
+      pred = GeoStatsModels.predict(sk, (:a, :b, :c), pset[i])
+      @test pred ≈ [j == i for j in 1:3]
+    end
   end
 
   @testset "Transiogram" begin