I don't know if there's an official way, given that the outputs from legend() beyond the first one are not mentioned in the official documentation. And you're right, when you specify NumColumns and request the 2nd (or subsequent) output, NumColumns seems to be ignored, so the accepted answer to the linked-to question doesn't work in this case.
Here's a workaround (using the approach taken by the other answer to the linked-to question) that may be good enough, which decreases the length of all the lines in the legend.
plot(rand(4))
lgd = legend(["Directional Data","Median Direction",...
"Interquartile Range","Interquartile Slope Data"],...
"NumColumns",4,"Interpreter","latex");
lgd.Position = [0.4,0.0,0.2,0.05];
lgd.ItemTokenSize(1) = 12;
