Subtracting the pixel value from the actual or expected value of the dark object is the idea. I will advice that you do that on the natural DN before converting to TOA reflectance. And since you are using QGIS for the analysis, I will suggest that you use the "Semi Auto Classification Plugin" to do the raster processing.
Fundamentally the question here is "what does 'scientifically valid' mean". If you are looking to do spectral modelling on the data, then the answer is possibly different than if you are looking at doing classification / image segmentation. Pansharpening (depending on the method) is simply going to change the range of the values a fairly small amount and shouldn't put your reflectance values outside the realm of possibility.
All in all, it depends a lot on what application you are going to be using the data for. Furthermore, the impact of pansharpening may also be worth documenting as a partial side result in whatever study you are performing. The result may be that it doesn't add anything, except four times as many pixels, meaning four times as long a processing time, which in some cases is a showstopper.
Edit: My database of articles on this topic is not huge, but I have these two where pansharpend data is used (with reasonable results) for image segmentation:
Shackelford, A. K., & Davis, C. H. (2003). A combined fuzzy pixel-based and object-based approach for classification of high-resolution multispectral data over urban areas. IEEE Transactions on Geoscience and Remote Sensing, 41(10), 2354–2364. http://doi.org/10.1109/TGRS.2003.815972
Fernández, I., Aguilar, F. J., Aguilar, M. A., & Álvarez, M. F. (2014). Influence of Data Source and Training Size on Impervious Surface Areas Classification Using VHR Satellite and Aerial Imagery Through an Object-Based Approach. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 7(12), 4681–4691.
Best Answer
The two products are comparable with some initial considerations. As you can see, the resolutions of the products are more or less the same:
To give a more visual comparison of the spectral bands covered by the different sensors, the handbook provides another figure:
From these information, you can guess, most of the bands are directly comparable based on their resolutions, and spectral bands, you just have to compare the correct bands (e.g. ETM+ B1 = OLI B2).
Furthermore, the two sensors record the observed data in the same units (both of them are spectrometers after all): absolute radiance.
From Landsat 8 handbook p. 61.
From Landsat 7 handbook.
The only difference (from this perspective) is in the quantization algorithm they use. To compress the huge amount of raw data (32bit float), the original values are quantized (16bit int), which is called digital numbers (DNs). From the DN rasters, you can compute the original values with a scale, and an offset.
The real difference between the quantizing methods, is Landsat 7 level 1 data only provide a scale, and an offset to calculate radiance values. With it, based on the Earth-Sun distance in astronomical units in the date of acquisition, and a mean solar exoatmospheric irradiance, which is a pre-calibration constant for every non-thermal band, you can calculate the reflectance value for the given scene (without correction for the solar elevation angle).
On the other hand, Landsat 8 level 1 products offer another scale and offset, with which you can easily calculate the reflectance value directly without the extra hassle. You just have to correct the value for the solar elevation angle, which is basically a division with sin(θ).
Note:
Do not be confused with the cos(θs) provided in the Landsat 7 handbook for solar elevation correction. θs is zenith angle, which is complementary with θ. Thus cos(θs) = sin(θ).