From what I understand about the plm package, those two approaches should be identical.

However, the fixed effects produced from this explicit specification are shown to be "reference dependent" [i.e. relative to the default reference in your factor(index)]

    tfe <- plm(y ~ x1 + x2 + factor(index), data, model = "within", index = c("id", "index"))

In contrast, fixef() returns the fixed effects in levels (by default). For you to get the same fixed effect estimates, by specifying the following:

    fixef(object = tfe, effect = "individual", type = "dfirst")

The equivalent for the individual level fixed effects would be:

    fixef(object = tfe, effect = "time", type = "dfirst")

From what I understand about the plm package, those two approaches should be identical.

However, the fixed effects produced from this explicit specification are shown to be "reference dependent" [i.e. relative to the default reference in your factor(index)]

    tfe <- plm(y ~ x1 + x2 + factor(index), data, model = "within", index = c("id", "index"))

In contrast, fixef() returns the fixed effects in levels (by default). For you to get the same fixed effect estimates, by specifying the following:

    fixef(object = tfe, effect = "individual", type = "dfirst")

The equivalent for the individual level fixed effects would be:

    fixef(object = tfe, effect = "time", type = "dfirst")

Computing R-Squared
Also, please see this post for computing R^2 and Adjusted R^2 manually for the full model (i.e. including both the fixed and specified effects): http://karthur.org/2016/fixed-effects-panel-models-in-r.html