Usually raster dowscaling methods exploit interpolation techniques. However, I am not sure how to proceed (in R, QGIS, GEE, or Python) if a raster needs to be downscaled (=i.e. its resolution increased, e.g. from 1 km to 100 m) such that the sum of the underlying pixels yields the value of the original lower resolution pixel. Of course I would like to to this in a way that not all pixels underlying the coarser resolution pixel have the same value, but based on the distribution of pixels on another raster layer (i.e. with a statistical regression approach).
the easiest approach is to increase the resolution of your raster (i think downscaling is an awkward term myself) and then to perform the regression or redistribution or further analysis on your new raster with the secondary data.
When restricting the higher resolution cells to be the sum of the original cell, in R it's very straightforward;
require(raster) # a dummy raster of 10 x 10 cells r <- raster(xmn=0, xmx=100, ymn=0, ymx=100, ncol=10, nrow=10) r <- sample(1:4, 100, replace=T) # disaggregation factor, in this case we're going to increase the res by 10 fact <- 10 # disaggregate and increase resolution # by not setting an interpolation method, we force the new cells to hold the value of the 'parent' cell r.inc <- disaggregate(r, fact) # so you need to divide the higher res cells by the disaggregation factor squared r.inc <- r.inc/fact^2 # check > cellStats(r,sum) == cellStats(r.inc,sum)  TRUE
then go on to perform your regression, glm etc using a stack