Given a specific value for a variable, how do you find the predicted value of a fixed effects multivariate regression? I have specified the following model in R:
model <- lm(yield ~ N + N^2 + P + K + S + factor(year) + factor(variety), 
            data=data)

As shown above, I have included fixed effects for year and variety. I want to find the predicted yield from my model, given a specific value (N = 100)for nitrogen for which I have included linear and squared terms.
I have been trying to do this with the predict function but am unsure how to do this with fixed effects and a squared term.

You have to use I(N^2) instead of N^2, because otherwise lm() interprets N^2 as crossing, which is not what you want here.
