Hand computation for plm package in R for predicted values

I am trying to fully understand the ins and outs of the plm package in R. I have a dataset where I computed a fixed effects plm:

mydata <- read.spss("G:/data.sav",use.value.labels=TRUE, to.data.frame = TRUE)

attach(mydata)
Y <- cbind(Y)
X <- cbind(x1,x2,x3)

pdata <- plm.data(mydata, index=c("id","YEAR"))
fixed <- plm(Y ~ X, data=pdata, model= "within")


I am trying to calculate both the training and test dataset predicted values. I found some material on CV here that kind of addresses what I am trying to do, but does not completely answer it. When I calculate by hand:

y = beta1*z + beta2*z + theta*id


the predicted values from

fitted <- as.numeric(fixed$model[[1]] - fixed$residuals)


are not the same as when I sum beta*x1 + beta*x2 + beta*x3 + fixef(fixed). Basically I am trying to figure out how to calculate by hand the in sample predicted values so I can run the calculations for the out of sample data.

Can anyone explain what I maybe doing wrong in doing the calculation by hand?

You need to extract the fixed effects by fixef and match them to the individual index. Here is an example for the Grunfeld data:

data(Grunfeld, package = "plm")
fe <- plm(inv ~ value + capital, data=Grunfeld, model = "within")

temp <- merge(Grunfeld, data.frame(fixef_firm = names(fixef(fe)), fixef = as.numeric(fixef(fe))), all.x =T, by.x = c("firm"), by.y=c("fixef_firm"))
fitted_by_hand <- temp$fixef + fe$coefficients[1] * Grunfeld$value + fe$coefficients[2] * Grunfeld$capital fitted <- fe$model[ , 1] - fe\$residuals

# just to remove attributs and specific classes
fitted_by_hand <- as.numeric(fitted_by_hand)
fitted <- as.numeric(fitted)

all.equal(fitted, fitted_by_hand) # TRUE
cbind(fitted, fitted_by_hand) # see yourself