Comparison with different resolution rasters

I am not clear exactly what you are wanting to test here. Do you have a hypothesis? You can certainly take a sampling approach in evaluating the "accuracy" of the lower resolution binominal classification along with hypothesis testing the distributional equivalence of the binominal classes.

Here is a worked example, in R, of what testing a hypothesis and evaluating accuracy would look like. The example is provided in R because ArcGIS just does not have the tools to perform this type of analysis.

First, let's add our required libraries and create some example data. Please note that the higher resolution raster contains random values on a proportional scale and the lower resolution data is binominal thus, emulating your problem. The results in the example cannot be interpreted literally due to the binominal data being a thresholded function of the higher resolution proportional data but, it gives the analytical framework for real data.

library(raster)
library(sp)
library(exactextractr)
library(rfUtilities)

r.high <- raster(nrows=180, ncols=360, xmn=571823.6, xmx=616763.6, ymn=4423540, 
            ymx=4453690, resolution=100, crs = CRS("+proj=utm +zone=12 +datum=NAD83 
            +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"))
    r.high[] <- runif(ncell(r.high),0,1)
r.low <- aggregate(r.high, fact=6, fun=median)
  r.low <- reclassify(r.low, matrix(c(0,0.5,0,0.5001,1,1), 
                      ncol=3, byrow=TRUE))

par(mfrow=c(1,2))
  plot(r.high)
  plot(r.low)

Now that we have a reproducible example we can create a random sample, across the higher resolution raster (which will contain the associated cell values) and then extract the values for the low-resolution data.

rs <- sampleRandom(r.high, 5000, sp=TRUE)
  rs@data <- data.frame(rs@data, extract(r.low, rs))
    names(rs) <- c("high_res", "low_res")

For hypothesis testing, using the Kruskal-Wallis Test, we can decide whether the population distributions are identical without assuming a normal distribution. The null hypothesis is that the snow proportions in the higher resolution data are identical populations across the binominal [0,1] values.

kruskal.test(rs$high_res ~ rs$low_res, data = rs) 

We can also use the log loss statistic to evaluate the accuracy, of the low-resolution data, based on penalizing the magnitude of the "miss-classification". Make sure that, if you are dealing with percentages, that you rescale to a 0-1.

logLoss(y=rs$low_res, p=rs$high_res)

We can also evaluate the data by splitting out the two classes. The Kolmogorov-Smirnov test evaluates the distributional equivalence and Mann-Whitney is a nonparametric test where the null hypothesis is that for randomly selected values x and y from two populations, the probability of x being greater than y is equal to the probability of y being greater than x.

ks.test(rs@data$high_res[which(rs$low_res == 1)], 
        rs@data$high_res[which(rs$low_res == 0)])
wilcox.test(rs@data$high_res[which(rs$low_res == 1)], 
            rs@data$high_res[which(rs$low_res == 0)], 
            alternative = "g")

Alternately you can just intersect the data by coercing the lower resolution data to polygons and extracting the values from the higher resolution data. You can then summarize the values or relate the data into a dataset and calculate an evaluation statistic.

Here, we coerce to polygons and extract high-resolution raster values to low-resolution polygons (pixels)

r.low <- rasterToPolygons(r.low, dissolve=FALSE)
  v <- exact_extract(r.high, as(r.low, "sf"))

Now, simply calculate a few statistical moments

v.stats <- lapply(v, function(x) summary(x[,1]))
  head( v.stats <- as.data.frame(do.call(rbind, v.stats)) )

Or, we can also relate values and calculate log loss

v.df <- do.call(rbind, v)
 v.df <- data.frame(v.df, r.low$layer)
    v.df <- data.frame(ID=unlist(lapply(strsplit(rownames(v.df), "[.]"), 
                        function(x) x[1])), v.df)   

r.low@data <- data.frame(ID=1:nrow(r.low), r.low@data) 
  names(r.low) <- c("ID", "snow_class")
  
v.df <- merge(v.df,r.low@data, by="ID", all = TRUE)   
  names(v.df) <- c("ID", "snow.pct", "pixel_fraction", "snow_class")
 
logLoss(v.df[,"snow_class"], v.df[,"snow.pct"])

To account for potential misalignment issues you can filter the resulting extracted data by the "coverage_fraction" column, or as I later change it to "pixel_fraction". This is the fraction of the intersection of each given pixel and the polygon. So, a 0.2 would be a 20% overlap in the intersection. You can choose a threshold and remove any pixel value that does not adhere to the fractional overlap.