How to calculate the robust standard error of predicted y from a linear regression model in R? How can I calculate the robust standard error of the predicted y from a linear regression model in R? Any suggestion is appreciated.
 A: I assume that you mean heteroskedasticity-consistent standard errors.  
From a fitted regression model, a predicted value is 
$$
\tilde y = \tilde X'\hat\beta
$$
Its variance is
$$
V(\tilde y) = V(\tilde X'\hat\beta)\\
V(\tilde y) = \tilde X' \hat V_\hat\beta \tilde X
$$
where $\hat V_\hat\beta$ is the estimated parameter covariance matrix.  The formula for a heteroskedasticity-consistent parameter covariance matrix is on wikipedia.  The standard errors of the fitted values are the square root of the diagonals of this matrix.
Here is an example:
#Fake data
x1 = rnorm(100)
x2 = rnorm(100)
e = x1*rnorm(100)
y = 10+x1-x2+e
X = cbind(1,x1,x2)

#Linear model
m = lm(y~X-1)
summary(m)
betahat = as.matrix(coef(m))

#Non-HC standard errors of fitted values
se.yhat = sqrt(diag(X%*%vcov(m)%*%t(X)))
se = predict(m,se.fit=TRUE)$se.fit
all.equal(se,se.yhat)#Matrix formula gives same result as se.fit option of predict method for lm's

#Now getting se's of fitted values based on HC-robust paramter covariance matrix
ehat = residuals(m)
vcov.HC = solve(t(X)%*%X) %*% t(X)%*%diag(ehat^2)%*%X %*% solve(t(X)%*%X)

se.yhat.HC = sqrt(diag(X%*%vcov.HC%*%t(X)))

#Showing that this method replicates regression parameter standard errors as given by the coeftest function
library(lmtest)
library(sandwich)
coeftest(m,vcov=vcovHC(m,type="HC0"))
sqrt(diag(vcov.HC))

The example here gives SE's for observed values, but you can put whatever you want into the X matrix, including a 1xp vector for a single observation.
You can't modify a lm object's $V_\beta$ matrix, because it isn't stored, rather it's called on the fly as requested.  But you can do so if you use gam from the mgcv package (you don't need to use semiparametrics necessarily)
#Respecifying as a gam in order to create an object from which you can use the predict method
library(mgcv)
mg = gam(y~X-1)
mg$Vp <- vcov.HC
    summary(mg)
    all.equal(as.numeric(predict(mg,se.fit=TRUE)$se.fit),se.yhat.HC)

This is convenient because you can now use the standard predict method with se.fit=TRUE
