It doesn't seem as they have implemented it directly; you'll have to compute the residuals as the difference between the resulting model prediction and the input data at the desired data points...
There's an example for the Fisher iris data set that doesn't specifically compute the residuals themselves, but calculates the predicted values and plots them in comparison to the input data...just follow its lead...
load fisheriris
t = table(species,meas(:,1),meas(:,2),meas(:,3),meas(:,4), ...
'VariableNames',{'species','meas1','meas2','meas3','meas4'});
Meas = dataset([1 2 3 4]','VarNames',{'Measurements'});
rm = fitrm(t,'meas1-meas4~species','WithinDesign',Meas);
Yhat=predict(rm,t([1 51 101],:))
Y=t([1 51 101],[2:end])
Res=Y-Yhat
Use your data a results similarly...
