mymodel <- (lm(Y ~ X1 + X2 + I(X1*X2), data = mydata));

clusteredSE <- function(fm, dfcw, cluster)
{
  library(sandwich)
  library(lmtest)
  library(zoo)
  M <- length(unique(cluster))
  N <- length(cluster)
  dfc <- (M / (M - 1)) * ((N - 1) / (N - fm$rank))
  u <- apply(estfun(fm), 2, function(x) tapply(x, cluster, sum))
  vcovCL <- dfc * sandwich(fm, meat = crossprod(u) / N) * dfcw
  coeftest(fm, vcovCL)
}

clusteredSE(mymodel, 1, mydata$countries)