X1 = [0 0 0 0 25 25 25 25 50 50 50 50]';
X2 = [ 0 10 20 30 0 10 20 30 0 10 20 30]';
Y = [0.51 0.22 0.205 0.234 0.73 0.214 0.212 0.225 0.87 0.234 0.202 0.204]';

scatter3(X1, X2, Y)
n = length(Y);

X = [ones(n,1) X1 X2]; 
szX = size(X);
p = szX(2);   % total number of beta's

betahat = (X' * X)^(-1) * X' * Y;
Yhat = X * betahat;    % predicted values 
resid = Y - Yhat;
SSE = sum(resid.^2)
SST = sum(( Y - mean(Y)).^2);
s2 = SSE/(n-p)
R2 = (SST-SSE)/SST


 hold on
plot3(X1, X2, Yhat, '*')
 hold off
 
 
 Xquad = [ones(n,1) X1 X2 X1.^2  X2.^2  X1.*X2];  
 
szXquad = size(Xquad);
p2 = szXquad(2); 
 
 beta2hat = (Xquad' * Xquad)^(-1) * Xquad' * Y; 
 Yhat2 = Xquad * beta2hat;
 resid2 = Y - Yhat2;
 
SSE = sum(resid2.^2)     % same as  resid2'*resid2
SST = sum(( Y - mean(Y)).^2);
s = sqrt(SSE/(n-p2))
R2 = (SST-SSE)/SST
 
 scatter3(X1, X2, Y)
   hold on
 plot3(X1, X2, Yhat2, '*')
   hold off
   
 
% residual plot 
plot3(X1,X2, resid2,'.')