Nice code. For someone who started using SAS a few weeks ago, you have made great progress.
Before I discuss the error, let me show you how the SAS log helps you find it. When you run the program, the log shows something like this:
747 train_ind = setdif(s,test_ind);
748 print train_ind;
749 trainset=x[train_ind, ];
750 print trainset;
751 testset =x[test_ind, ];
752 print testset;
753 testset_y = y[test_ind];
754 trainset_y = y[train_ind];
755 run Regress_cv;
756 mse_acum=mse_acum + mse;
757 end;
...
769 end;
ERROR: (execution) Invalid subscript or subscript out of range.
operation : [ at line 749 column 17
operands : x, train_ind,
x 37 rows 2 cols (numeric)
train_ind 1 row 37 cols (numeric)
statement : ASSIGN at line 749 column 7
Now look at the red line near the bottom of the log. It says that the invalid subscript occurred in the '[' operation (subscript) on line 749, column 17 of the log. If you look at the line numbers in the log, you discover that Line 749 is
749 trainset=x[train_ind, ];
If you also print 'j' (the iteration number) in the loop, you discover that the problem occurs when j=2. The question is therefore why the X variable doesn't have enough rows to support the subscript operation for iteration 2?
The answer is that the Regress_cv module does not have any arguments and therefore all variables in the module are GLOBAL variables! (See the doc.) In particular, there is a variable x in the module that overwrites the symbol of the same name in your program.
One solution is to make the module use a local symbol table by passing in the training and test data set. For example, you might define the module as
start Regress_cv(_trainset, _testset);
k = ncol(_trainset);
trainset_y = _trainset[,1]; trainset = _trainset[,2:k];
testset_y = _testset[,1]; testset = _testset[,2:k];
ntrain = nrow(trainset); /* trainset size */
c = j(ntrain, 1, 1); /* column vector of 1s */
x = c||trainset; /* concatenating */
y = trainset_y;
xpxi = inv(x`*x); /* inverse of X'X */
beta = xpxi * (x`*y); /* parameter estimate */
ntest = nrow(testset); /* testset size */
a = j(ntest, 1, 1); /* column vector of 1s */
testset = a||testset;
yhat = testset*beta; /* prediction on testset*/
resid = testset_y-yhat; /* residuals */
sse = ssq(resid); /* SSE */
dfe = nrow(x)-ncol(x); /* error DF */
mse = sse/dfe; /* MSE */
return mse;
finish Regress_cv;
Notice that I've passed in the whole train/test data. You could pass in 4 args if you prefer to subdivide the data in the main program. Using this new module, call it like this:
do i = 2 to 16;
mse_acum=0;
x = pcs_matrix[,2:i];
y = pcs_matrix[,1];
n = nrow(x); /* sample size */
ssquares= ssq(y-sum(y)/n);
aux= 1:n; /* vector from 1:n */
call randseed(1); /* reproducible results*/
s = sample(aux, n, "WOR"); /* sampling without replacement */
do j = 1 to (k-1); /* creating train and test sets and regressing */
test_ind = s[((j-1)#ceil(n/k) + 1):(j#ceil(n/k))];
train_ind = setdif(s,test_ind);
trainset=pcs_matrix[train_ind, 1:i];
testset =pcs_matrix[test_ind, 1:i];
mse = Regress_cv(trainset, testset);
mse_acum=mse_acum + mse;
end;
/* for the last fold since n it's not always divisible by k do:*/
test_ind = s[((k-1)#ceil(n/k) + 1):n];
train_ind = setdif(s,test_ind);
trainset=pcs_matrix[train_ind, 1:i];
testset =pcs_matrix[test_ind, 1:i];
mse = Regress_cv(trainset, testset);
mse_acum=mse_acum + mse;
print mse_acum;
*r2_cross[i-1] = 1 - mse_acum#n/ssquares;
end;
For more about local and global symbols in SAS/IML, see the article "Understanding local and global variables in the SAS/IML language".
