plm: using fixef() to manually calculate fitted values for a fixed effects twoways model

This works for an unbalanced data with effect="individual" and time dummies y ~ x +factor(year):

fitted <- pmodel.response(plm.model)-residuals(plm.model)