diff --git a/src/1.JWAS/src/JWAS.jl b/src/1.JWAS/src/JWAS.jl
index 634d5e26..58ee13f6 100644
--- a/src/1.JWAS/src/JWAS.jl
+++ b/src/1.JWAS/src/JWAS.jl
@@ -134,7 +134,7 @@ export dataset
         * If `causal_structure` is provided, e.g., causal_structure = [0.0 0.0 0.0;1.0 0.0 0.0;1.0 0.0 0.0] for
           trait 2 -> trait 1 and trait 3 -> trait 1 (row index affacts column index), phenotypic causal networks will be incorporated using structure equation models.
 * Genomic Prediction
-    * Predicted values for individuals of interest can be obtained based on an user-defined prediction equation `prediction_equation`, e.g., "y1:animal + y1:geno + y1:age".
+    * Predicted values for individuals of interest can be obtained based on an user-defined prediction equation `prediction_equation`, e.g., "y1:animal + y1:age".
     For now, genomic data is always included. Genetic values including effects defined with genotype and pedigre information are returned if `prediction_equation`= false, defaulting to `false`.
     * Individual estimted genetic values and prediction error variances (PEVs) are returned if `outputEBV`=true, defaulting to `true`. Heritability and genetic
     variances are returned if `output_heritability`=`true`, defaulting to `true`. Note that estimation of heritability is computaionally intensive.
diff --git a/src/1.JWAS/src/output.jl b/src/1.JWAS/src/output.jl
index e443421e..73e37f5d 100644
--- a/src/1.JWAS/src/output.jl
+++ b/src/1.JWAS/src/output.jl
@@ -26,12 +26,20 @@ function prediction_setup(model)
         end
     else
         prediction_equation = string.(strip.(split(model.MCMCinfo.prediction_equation,"+")))
-        if mme.MCMCinfo.output_heritability != false
+        if model.MCMCinfo.output_heritability != false
             printstyled("User-defined prediction equation is provided. ","The heritability is the ",
             "proportion of phenotypic variance explained by the value defined by the prediction equation.\n",
             bold=false,color=:green)
         end
+        for i in prediction_equation
+            term_symbol = Symbol(split(i,":")[end])
+            if !(haskey(model.modelTermDict,i) || (isdefined(Main,term_symbol) && typeof(getfield(Main,term_symbol)) == Genotypes))
+                error("Terms $i in the prediction equation is not found.")
+            end
+        end
     end
+    printstyled("Predicted values for individuals of interest will be obtained as the summation of ",
+    prediction_equation, " (Note that genomic data is always included for now).",bold=false,color=:green)
     if length(prediction_equation) == 0 && model.M == false
         println("Default or user-defined prediction equation are not available.")
         model.MCMCinfo.outputEBV = false
@@ -221,10 +229,13 @@ function getEBV(mme,traiti)
         mytrait, effect = split(term,':')
         if mytrait == traiti_name
             sol_term     = map(Float64,location_parameters[(location_parameters[!,:Effect].==effect).&(location_parameters[!,:Trait].==traiti_name),:Estimate])
-            if length(sol_term) == 1 #1-element Array{Float64,1} doesn't work below; convert it to a scalar
+            if VERSION < v"1.6" && length(sol_term) == 1 #1-element Array{Float64,1} doesn't work below; Will be deleted
                 sol_term = sol_term[1]
             end
             EBV_term = mme.output_X[term]*sol_term
+            if VERSION < v"1.6" && length(sol_term) == 1 #1-element Array{Float64,1} doesn't work below; Will be deleted
+                EBV_term = vec(EBV_term)
+            end
             EBV += EBV_term
         end
     end