I am having difficulty in fitting a model on data. Basically, I have data about the evaluation of phenotypic property (i.e. hard) of 65 palm trees by 5 judges. As an evaluation scheme, each judge provides score to each sample. For 3 judges sample data look like this:
Judge Product Hard
aa 1 5
ab 1 6
ac 1 3
aa 1 7
ab 1 5
ac 1 4
aa 2 5
ab 2 8
ac 2 6
aa 2 7
ab 2 4
ac 2 4
Main objective here is to get product coefficients with less judge errors, for which I want to fit this kind of model:
$$Y_{ij} = α_i + β_iθ_j + ε_{ij}$$ i = judge, j = product
Here, $α_i$ is judge main coefficients, $_i$ is judge coefficients due to difference in their scoring pattern and $θ_j$ is product coefficients and $ε_i$ is assessor dependent.
I was trying to fit this model using lme function in R, but difficulty I am facing to fit the interaction term because model here fitted for parameters rather than co-variates.
This model looks quite accurate for my kind of data. I have seen Bayesian version (http://www.r-bloggers.com/extending-the-sensory-profiling-data-model/) of it and I don't know how to do using mixed-modelling approach or in a frequentist way.
My queries here are:
a) What can be an appropriate method to fit this kind of model? I had referred so much literature where description about iterative generalized least squares, multi-level model, separate regression model, weighted least-square model are given. But still I am not getting how to use and fit estimated value of parameters in interaction terms and get separate coefficients for both interaction parameters?
b) How can I get heterogeneous error in this form?
c) which R package can I use?
