How to model multiplicative effect of parameters/fit data at individual predictor level -

i having difficulty in fitting model on data. basically, have data evaluation of phenotypic property (i.e. hard) of 65 palm trees 5 judges. evaluation scheme, each judge provides score each sample. 3 judges sample data 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   yij=αi+βiθj+εij = judge, j = product 

here, αi judge main coefficients, i judge coefficients due difference in scoring pattern , θj product coefficients , εi assessor dependent.

i trying fit model using lme function in r, difficulty facing fit interaction term because model here fitted parameters rather co-variates.

this model looks quite accurate kind of data. have seen bayesian version (http://www.r-bloggers.com/extending-the-sensory-profiling-data-model/) of , don't know how using mixed-modelling approach or in frequentist way.

my queries here are:

a) can appropriate method fit kind of model? had referred literature description iterative generalized least squares, multi-level model, separate regression model, weighted least-square model given. still not getting how use , fit estimated value of parameters in interaction terms , separate coefficients both interaction parameters?

b) how can heterogeneous error in form?

c) r package can use?