multivariate polynomial fitting

[Alexia369]

loss_function = 
    """  

    function my_loss_function(tree, dataset::Dataset{T,L}, options)::L where {T,L}

        tree.degree != 2 && return L(Inf)

        P = tree.l
        Q = tree.r

        s1, flag = eval_tree_array(P, dataset.X, options)
        !flag && return L(Inf)

        s2, flag = eval_tree_array(Q, dataset.X, options)
        !flag && return L(Inf)

        J_fit = s1 .^ 2 + s2 .^ 3

        loss = sum((J_fit .- dataset.y) .^ 2)

        return loss
    end
    """

J_fit = s1 .^ 2 + s2 .^ 3, I hope to change s1, s2 to get dataset.y by multivariate polynomial fitting, just like univariate polynomial regression

f = Polynomials.fit(s, dataset.y, 2)
J_fit = f.(s)

How to write Julia code to implement multivariate polynomial regression in Loss Function

[Alexia369]

loss_function = 
    """ 
    using GLM
    using DataFrames

    function my_loss_function(tree, dataset::Dataset{T,L}, options)::L where {T,L}

        tree.degree != 2 && return L(Inf)

        P = tree.l
        Q = tree.r

        s1, flag = eval_tree_array(P, dataset.X, options)
        !flag && return L(Inf)

        s2, flag = eval_tree_array(Q, dataset.X, options)
        !flag && return L(Inf)
        
        data = DataFrame(X1=s1, X2=s2, y=dataset.y)
        model = GLM.lm(@formula(y ~ X1 + X2 + I(X1^2) + I(X1*X2) + I(X2^2)), data)
        J_fit = predict(model, data)

        loss = sum((J_fit .- dataset.y) .^ 2)

        return loss
    end
    """

I used the GLM library, but it didn't work. Is it a code problem?