We need priors for the four regression parameters, and for the error variance. The BUGS/JAGS framework does not let us use improper priors, so we have to pick things that we think are "wide" enough.
For the error variance, we also have to pick something that's always positive.
Finally, BUGS/JAGS thinks in terms of "precision" $\tau =1/{\sigma}^{2}$, which is confusing, but you can get used to it. It's best to start by guessing what you think a "large" standard deviation would be, then taking $1/{\sigma}^{2}$ to get the relevant value of $\tau $. A plausible rule of thumb might be 10 times the expected value of your parameter. For example, if you expected that the difference between masses of two species might be 100g, you could use a prior with a standard deviation of 1000g or a precision of ${10}^{-6}$. The normal priors below have low precision, so they have high variance.
rjags
R2jags
You can use inits=NULL to have JAGS pick the initial values randomly from the priors. For more complex models you might want to pick starting values for the chains yourself (see the ?jags documentation).
Examine the chains and output
The lump.mcmc.list function from the emdbook package can be useful for converting a set of MCMC chains into a single long chain.
Further notes
Categorical variables
Implementing models with categorical variables (say, a t-test or an ANOVA) is a little bit more tedious than the multiple regression analysis shown above. There are two basic strategies:
pass the categorical variable as a vector of numeric codes (i.e. pass as.numeric(f) rather than f to JAGS in your data list), and make your parameters a simple vector of the means corresponding to each level, e.g. you could have a vector of parameters catparam and specify the predicted value as pred[i] = catparam[f[i]]
construct a design matrix using the model.matrix function in R and pass the whole model matrix (X)to JAGS. Then, to get the predicted value for each observation, just add up the relevant columns of the model matrix: e.g. pred[i] = beta[1]*X[i,1]+beta[2]*X[i,2]+beta[3]*X[i,3]. You can also use the built-in inprod (inner product) function: pred[i] = inprod(X[i,],beta)
The second approach is a little bit harder to understand but generalizes to more complicated situations, and gives you answers that will more closely match the analogous analyses in R (e.g. using lm).