first commit

coxnet.m

@@ -0,0 +1,76 @@
+function fit = coxnet(x,is_sparse,irs,pcs,y,weights,offset,parm,nobs,nvars,...
+    jd,vp,cl,ne,nx,nlam,flmin,ulam,thresh,isd,maxit,family)
+
+% Internal glmnet function. See also glmnet, cvglmnet.
+
+%time --- column 1
+%status --- column 2
+
+ty = y(:,1);
+tevent = y(:,2);
+if (any(ty <= 0))
+    error('negative event times encountered;  not permitted for Cox family');
+end
+if isempty(offset)
+    offset = ty * 0;
+    is_offset = false;
+else
+    is_offset = true;
+end
+
+if (is_sparse)
+    error('Cox model not implemented for sparse x in glmnet');
+else
+    task = 41;
+    [lmu,ca,ia,nin,dev,alm,nlp,jerr,dev0,ot] = glmnetMex(task,parm,x,ty,jd,vp,ne,nx,nlam,flmin,ulam,thresh,isd,weights,tevent,cl,maxit,offset);
+
+end
+
+if (jerr ~= 0)
+    errmsg = err(jerr,maxit,nx,family);
+    if (errmsg.fatal)
+        error(errmsg.msg);
+    else
+        warning(errmsg.msg);
+    end
+end
+
+ninmax = max(nin);
+lam = alm;
+if (ulam == 0.0)
+    lam = fix_lam(lam); % first lambda is infinity; changed to entry point
+end
+
+dd=[nvars, lmu];
+if ninmax > 0
+    ca = ca(1:ninmax,:);
+    df = sum(abs(ca) > 0, 1);
+    ja = ia(1:ninmax);
+    [ja1,oja] = sort(ja);
+    beta = zeros(nvars, lmu);
+    beta (ja1,:) = ca(oja,:);
+else
+    beta = zeros(nvars,lmu);
+    df = zeros(1,lmu);
+end
+
+fit.beta = beta;
+fit.dev = dev;
+fit.nulldev = dev0;
+fit.df = df';
+fit.lambda = lam;
+fit.npasses = nlp;
+fit.jerr = jerr;
+fit.dim = dd;
+fit.offset = is_offset;
+fit.class = 'coxnet';
+
+
+function new_lam = fix_lam(lam)
+
+new_lam = lam;
+if (length(lam) > 2)
+    llam=log(lam);
+    new_lam(1)=exp(2*llam(2)-llam(3));
+end
+

cvcompute.m

@@ -0,0 +1,22 @@
+function cvcpt = cvcompute(mat, weights, foldid, nlams)
+
+% Internal glmnet function. See also cvglmnet.
+
+% Compute the weighted mean and SD within folds, and hence the se of the mean
+
+wisum = accumarray(foldid,weights);
+nfolds = max(foldid);
+outmat = NaN(nfolds,size(mat,2));
+good = zeros(nfolds,size(mat,2));
+mat(isinf(mat)) = NaN;
+for i = 1:nfolds
+    mati = mat(foldid==i,:);
+    wi = weights(foldid==i,:);
+    outmat(i,:) = wtmean(mati,wi);
+    good(i,1:nlams(i)) = 1;
+end
+N = sum(good,1);
+cvcpt.cvraw = outmat;
+cvcpt.weights = wisum;
+cvcpt.N = N;
+end

cvcoxnet.m

@@ -0,0 +1,87 @@
+function result = cvcoxnet(object, lambda, x, y, weights, offset, foldid, ...
+    type, grouped, keep)
+
+% Internal glmnet function. See also cvglmnet.
+
+if nargin < 10 || isempty(keep)
+    keep = false;
+end
+
+typenames = struct('deviance','Partial Likelihood Deviance');
+if strcmp(type, 'default')
+    type = 'deviance';
+end
+if ~any(strcmp(type, {'deviance'}))
+    warning('Only ''deviance''  available for Cox models; changed to type=''deviance''');
+    type = 'deviance';
+end
+
+if isempty(offset)
+    offset = zeros(size(x,1),1);
+end
+nfolds = max(foldid);
+
+if (length(weights)/nfolds < 10) && ~grouped
+    warning('Option grouped=true enforced for cv.coxnet, since < 3 observations per fold');
+    grouped = true;
+end
+cvraw = NaN(nfolds,length(lambda));
+
+for i = 1:nfolds
+    which = foldid == i;
+    fitobj = object{i};
+    coefmat = glmnetPredict(fitobj,[],[],'coefficients');
+    if (grouped)
+        plfull = cox_deviance([],y,x,offset,weights,coefmat);
+        plminusk = cox_deviance([],y(~which,:),x(~which,:),offset(~which),...
+            weights(~which),coefmat);
+        cvraw(i,1:length(plfull)) = plfull - plminusk;
+    else
+        plk = cox_deviance([],y(which,:),x(which,:),offset(which),...
+            weights(which),coefmat);
+        cvraw(i,1:length(plk)) = plk;
+    end
+end
+status = y(:,2);
+N = nfolds - sum(isnan(cvraw),1);
+weights = accumarray(reshape(foldid,[],1),weights.*status);
+cvraw = bsxfun(@rdivide,cvraw,weights);  %even some weight = 0 does matter because of adjustment in wtmean!
+cvm = wtmean(cvraw,weights);
+sqccv = (bsxfun(@minus,cvraw,cvm)).^2;
+cvsd = sqrt(wtmean(sqccv,weights)./(N-1));
+result.cvm = cvm; result.cvsd = cvsd; result.name = typenames.(type);
+if (keep)
+    warning('keep=TRUE not implemented for coxnet');
+end
+
+
+function result = cox_deviance(pred, y, x, offset, weights, beta)
+
+ty = y(:,1);
+tevent = y(:,2);
+nobs = length(ty); nvars = size(x,2);
+if isempty(weights)
+    weights = ones(nobs,1);
+end
+if isempty(offset)
+    offset = zeros(nobs,1);
+end
+if isempty(beta)
+    beta = []; nvec = 1; nvars = 0;
+else
+    nvec = size(beta,2);
+end
+
+task = 42;
+[flog, jerr] = glmnetMex(task,x,ty,tevent,offset,weights,nvec,beta);
+
+
+if (jerr ~= 0)
+    errmsg = err(jerr,0,0,'cox');
+    if (errmsg.fatal)
+        error(errmsg.msg);
+    else
+        warning(errmsg.msg);
+    end
+end
+result = -2 * flog;

cvelnet.m

@@ -0,0 +1,70 @@
+function result = cvelnet(object, lambda, x, y, weights, offset, foldid, ...
+    type, grouped, keep)
+
+% Internal glmnet function. See also cvglmnet.
+
+if nargin < 10 || isempty(keep)
+    keep = false;
+end
+
+typenames = struct('deviance','Mean-Squared Error','mse','Mean-Squared Error','mae','Mean Absolute Error');
+if strcmp(type, 'default')
+    type = 'mse';
+end
+if ~any(strcmp(type, {'mse','mae','deviance'}))
+    warning('Only ''mse'', ''deviance'' or ''mae''  available for Gaussian models; ''mse'' used');
+    type = 'mse';
+end
+
+if ~isempty(offset)
+    y = y - offset;
+end
+
+predmat = NaN(length(y),length(lambda));
+nfolds = max(foldid);
+nlams = nfolds;
+
+for i = 1:nfolds
+    which = foldid == i;
+    fitobj = object{i};
+    fitobj.offset = false;
+    preds = glmnetPredict(fitobj,x(which,:));
+    nlami = length(object{i}.lambda);
+    predmat(which,1:nlami) = preds;
+    nlams(i) = nlami;
+end
+
+N = size(y,1) - sum(isnan(predmat),1);
+
+yy = repmat(y, 1, length(lambda));
+switch type
+    case 'mse'
+        cvraw = (yy - predmat).^2;
+    case 'deviance'
+        cvraw = (yy - predmat).^2;
+    case 'mae'
+        cvraw = abs(yy - predmat);
+end
+
+if (length(y)/nfolds < 3) && grouped
+    warning('Option grouped=false enforced in cv.glmnet, since < 3 observations per fold');
+    grouped = false;
+end
+
+if (grouped)
+    cvob = cvcompute(cvraw,weights,foldid,nlams);
+    cvraw = cvob.cvraw; 
+    weights = cvob.weights;
+    N = cvob.N;
+end
+
+cvm = wtmean(cvraw,weights);
+sqccv = (bsxfun(@minus,cvraw,cvm)).^2;
+cvsd = sqrt(wtmean(sqccv,weights)./(N-1));
+
+result.cvm = cvm; result.cvsd = cvsd; result.name = typenames.(type);
+
+if (keep)
+    result.fit_preval = predmat;
+end
+end

cvfishnet.m

@@ -0,0 +1,78 @@
+function result = cvfishnet(object,lambda,x,y,weights,offset,foldid,type,grouped,keep)
+
+% Internal glmnet function. See also cvglmnet.
+
+if nargin < 10 || isempty(keep)
+    keep = false;
+end
+
+typenames = struct('mse','Mean-Squared Error','mae','Mean Absolute Error','deviance','Poisson Deviance');
+if strcmp(type, 'default')
+    type = 'deviance';
+end
+if ~any(strcmp(type, {'mse','mae','deviance'}))
+    warning('Only ''mse'', ''deviance'' or ''mae''  available for Poisson models; ''deviance'' used');
+    type = 'deviance';
+end
+
+is_offset = ~isempty(offset);
+
+predmat = NaN(length(y),length(lambda));
+nfolds = max(foldid);
+nlams = nfolds;
+
+for i = 1:nfolds
+    which = foldid == i;
+    fitobj = object{i};
+    if (is_offset)
+        off_sub = offset(which);
+    else
+        off_sub = [];
+    end
+    preds = glmnetPredict(fitobj,x(which,:),[],[],[],off_sub);
+    nlami = length(object{i}.lambda);
+    predmat(which,1:nlami) = preds;
+    nlams(i) = nlami;
+end
+
+N = size(y,1) - sum(isnan(predmat),1);
+
+yy = repmat(y, 1, length(lambda));
+switch type
+    case 'mse'
+        cvraw = (yy - predmat).^2;
+    case 'mae'
+        cvraw = abs(yy - predmat);
+    case 'deviance'
+        cvraw = devi(yy, predmat);
+end
+
+if (length(y)/nfolds < 3) && grouped
+    warning('Option grouped=false enforced in cv.glmnet, since < 3 observations per fold');
+    grouped = false;
+end
+
+if (grouped)
+    cvob = cvcompute(cvraw,weights,foldid,nlams);
+    cvraw = cvob.cvraw; 
+    weights = cvob.weights;
+    N = cvob.N;
+end
+
+cvm = wtmean(cvraw,weights);
+sqccv = (bsxfun(@minus,cvraw,cvm)).^2;
+cvsd = sqrt(wtmean(sqccv,weights)./(N-1));
+
+result.cvm = cvm; result.cvsd = cvsd; result.name = typenames.(type);
+
+if (keep)
+    result.fit_preval = predmat;
+end
+
+
+function result = devi(yy, eta)
+
+deveta = yy .* eta - exp(eta);
+devy = yy .* log(yy) - yy;
+devy(yy == 0) = 0;
+result = 2 * (devy - deveta);