You wrote:
For regressor R (which represents some condition x trial part) on subject S, I would like to extract the ROI's average value for the TRs during which R is "on" (with some padding to the left and right).

DP>> Do you mean the voxelwise average over the time window? Or the spatial average over the ROI during that window?

Then, I'd like to come up with one averaged timecourse that represents the average of all those individual timecourses. This represents the average value in the ROI during R (as a time course) for S.

DP>> From this part, it sounds like you're effectively trying to get the impulse response for that regressor averaged over the ROI. It might be easier to extract the impulse response directly with -iresp option to 3dDeconvolve then use 3dMaskAve on the iresp output to get the ROI average impulse response...