Consider the code below. It fits a multi-predictor model, but plots the fit as a function of one predictor at a time, with the others held fixed at their mean values. In contrast, your code has each predictor increase along with the others.
This has its own issues. For one thing, the superimposed points represent data, and may or may not be comparable to the line that constrains predictors to their means. But it is one of the things you can do to visualize a multi-predictor fit as a function of one variable at a time.
x = mvnrnd([0 0],[1 .9;.9 1],200);
y = 5*x(:,1) + 5*x(:,2) + randn(200,1);
[b,dev,stats] = glmfit(x,y,'normal');
% Create an array with all predictors at their mean
meanx = mean(x);
xpred = repmat(meanx,100,1);
xplot = xpred;
ypred = zeros(100,size(x,2));
for k= 1:size(x,2)
% Allow the kth predictor to vary, others stay fixed
xplot(:,k)= linspace(min(x(:,k)),max(x(:,k)),100);
xpred(:,k) = xplot(:,k);
ypred(:,k) = glmval(b,xpred,'identity',stats);
xpred(:,k) = meanx(k); % put the mean back
end
for j= 1:size(x,2)
figure(j);
plot(x(:,j),y,'k.',xplot(:,j),ypred(:,k),'b-','LineWidth',1)
end
If you run this you'll see that the lines don't match the points very well. That's because I generated the x values with correlation 0.9, meaning that it's unlikely that one predictor will stay fixed while the other increases. If you adjust 0.9 down toward 0.0, you should see better agreement between the fit and the data points.
