Getting (maybe manually) confidence interval of fits after using multi-way clustering package (multiwayvcov)

I am interested in plotting fits with confidence intervals after using two-way clustering package (multiwayvcov).

Here is my reproducible data.

rm(list=ls(all=TRUE))
library(lmtest)
library(multiwayvcov)

dv<-c(1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0)
int1<-c(0.0123, 0.3428, 0.2091, 0.8325, 0.7113, 0.7401, 0.6009, 0.5062, 0.4841, 0.8912, 0.3850, 0.2463, 0.0625, 0.5374, 0.1984)
int2<-c(0, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0)
cont<-c(3, 1, 2, 4, 6, 7, 1, 4, 3, 2, 4, 3, 6, 1, 3)
cluster1<-c(1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3)
cluster2<-c(1, 2, 3, 1, 2, 3, 1, 2, 1, 2, 1, 2, 3, 1, 2)

mydata<-as.data.frame(cbind(dv, int1, int2, cont, cluster1, cluster2))

This is my non-clustered model:

result_lm <- lm(dv~int1+int2+cont,data=mydata)

To get clustered results using "cluster1" and "cluster2", I use functions in the package of "lmtest" and "multiwayvcov" as follows.

cluster_vcov<-cluster.vcov(result_lm, ~cluster1+cluster2)
result_2c<-coeftest(result_lm, cluster_vcov) 

Here, "cluster_vcov" is just a variance-covariance matrix and "result_2c" is just an atomic vector. Thus, I am not able to use "predict" function to plot fits on a new dataset ("datagrid") such as

grid <- seq(0,1,.2)
datagrid <- data.frame(int1=rep(grid,2),
                   int2=c(rep(0,length(grid)),
                          rep(1,length(grid))))
datagrid$cont<-mean(mydata$cont, na.rm=T)

Before moving to what I have done, here is something similar what I would like to have eventually.

fits <- predict(result_lm,newdata=datagrid,interval="confidence")
plotdata <- data.frame(fits,datagrid)
plotdata$int2 <- plotdata$int2==1
ggplot(plotdata,aes(x=int1,y=fit,ymin=lwr,ymax=upr,color=int2)) + geom_line(aes(linetype = int2)) + geom_ribbon(alpha=.2) + theme(legend.position="none") + scale_color_manual(values=c("red", "darkgreen")) + scale_linetype_manual(values=c("dashed", "solid"))

The result is

To address the problem that "result_2c" does not give a dataframe that can be directly used with "predict", I decided to construct a data by myself as follows.

d_twc_result<-data.frame(matrix(0, nrow =4, ncol = 4) )  
colnames(d_twc_result) <- c("Estimate","Std. Error","t value", "Pr(>|t|)") 
rownames(d_twc_result) <-c("(Intercept)", "int1","int2", "cont")

for (j in 1:4){
  for (i in 1:4){
    d_twc_result[i, j]<-result_2c[i,j]
  }
}

Then, using "d_twc_result$Estimate", I generate a vector that corresponds to "fits" that one could get after running "predict".

fits<-c(1:12) 
for (i in 1:12){
  fits[i]<-d_twc_result$Estimate[1]+
    d_twc_result$Estimate[2]*datagrid$int1[i]+
    d_twc_result$Estimate[3]*datagrid$int2[i]+
    d_twc_result$Estimate[4]*datagrid$cont[i]
 }

Yet, I was still not able to construct vectors for "lwr" and "upr", which requires 'residuals' or 'standard error'. What I was actually stuck is that it seems impossible to get 'residuals' or 'standard error' because there is no observation on 'dv' in the dataset "datagrid".

Nevertheless, "predict" works with the dataset "datagrid", so I guess that I am poorly understanding how "predict" works or the concept of fit.

It will be highly appreciated if you could help me to get "lwr" and "upr" (if my understanding of the concept of fit is incorrect). Thank for any comment in advance.