I'm pretty sure that neither of your ideas (using weights or multiple instances of $x_h$) will work correctly. My understanding of the help page for predict.lm is that you need to use pred.var = summary(mod)$sigma^2 / m where mod is your model and m is your value of $m$.

  The reasoning is that pred.var is used to set the variance of future observations; by default it is assumed to be the same as in the data ($\sigma^2$, estimated by $MSE$). By pretending it is $\sigma^2/m$ (estimated as $MSE/m$), you are using the variance of the average of $m$ predictions and will get the correct result.

This idea is reinforced by looking at the code for predict.lm -- here is a fragment:

hwid <- tfrac * switch(interval, confidence = sqrt(ip), 
    prediction = sqrt(ip + pred.var))

Note that ip here is the same for both confidence and prediction intervals, and pred.var works independently of it.

Addendum

Duh! -- after looking at the fact that the default for pred.var is res.var/weights, it does work to simply set weights equal to the value of $m$.

The weights argument will work -- just set it equal to $m$ in your formula. It will not work to put multiple instances of $x_h$ in newdata -- you'll just get several copies of the same interval.

Another approach is to use pred.var = summary(mod)$sigma^2 / m where mod is your model and m is your value of $m$. The reason this works is that pred.var is used to set the variance of future observations; by default it is assumed to be the same as in the data ($\sigma^2$, estimated by $MSE$). By pretending it is $\sigma^2/m$ (estimated as $MSE/m$), you are using the variance of the average of $m$ predictions and will produce the correct result.

I'm pretty sure that neither of your ideas (using weights or multiple instances of $x_h$) will work correctly. My understanding of the help page for predict.lm is that you need to use pred.var and set it equal to $MSE/m$.

The reasoning is that pred.var is used to set the variance of future observations; by default it is assumed to be the same as in the data ($\sigma^2$, estimated by $MSE$). By pretending it is $\sigma^2/m$ (estimated as $MSE/m$), you are using the variance of the average of $m$ predictions and will get the correct result.

I'm pretty sure that neither of your ideas (using weights or multiple instances of $x_h$) will work correctly. My understanding of the help page for predict.lm is that you need to use pred.var = summary(mod)$sigma^2 / m where mod is your model and m is your value of $m$.

The reasoning is that pred.var is used to set the variance of future observations; by default it is assumed to be the same as in the data ($\sigma^2$, estimated by $MSE$). By pretending it is $\sigma^2/m$ (estimated as $MSE/m$), you are using the variance of the average of $m$ predictions and will get the correct result.

This idea is reinforced by looking at the code for predict.lm -- here is a fragment:

hwid <- tfrac * switch(interval, confidence = sqrt(ip), 
    prediction = sqrt(ip + pred.var))

Note that ip here is the same for both confidence and prediction intervals, and pred.var works independently of it.

Addendum

Duh! -- after looking at the fact that the default for pred.var is res.var/weights, it does work to simply set weights equal to the value of $m$.