function LUR_Analysis(onefield,randomtest,perremove,correlmatrix)
    %Be sure to update the filepaths below to your current directories.
    %Data files should be tab-delimited .txt files. Each independent
    %variable type (e.g., major roads, % impervious surfaces) should have
    %its own file, with each row indicating buffer length. All independent
    %variables should be kept in their own folder.

    % onefield = 1 -> one one buffer per variable can be selected
    % onefield = 0 -> multiple buffers per variable can be selected
    %onefield = 0;


    % randomtest = 0 -> use the entire dataset
    % randomtest = 1 -> remove a %(1 - perremove) of the dateset to test the model
    %randomtest = 0;
    %perremove = 0.90;


    % correlmatrix = 1 -> output the correlation matrix for all independant
    % variables
    % correlmatrix = 0 -> do not output the correlation matrix
    %correlmatrix = 0;

    %folder with all of the text files containing the independant variables
    %(update the filepath to folder containing LUR data)
    folder = uigetdir('C:\...','Select the folder where the indenpentent variables are located');


    % get the names of all files. dirListing is a struct array.
    dirListing = dir(folder);

    % loop through the files, open the files and load the data
    count = 0;
    textnames = struct('name','textdata');
    data1 = struct('name','data');
    for k = 1:length(dirListing)
        if ~dirListing(k).isdir

            %initialize the data arrays
            fileName = fullfile(folder,dirListing(k).name);
            d = importdata(fileName);
            count = count + 1;
            fieldname = strcat('field',int2str(count));
            p = d.textdata;
            textnames = setfield(textnames,fieldname,p);
            q = d.data;
            data1 = setfield(data1,fieldname,q);
            [m,n] = size(d.data);
            len(count) = n; 
        end
    end

    data1 = rmfield(data1, 'name');
    textnames = rmfield(textnames, 'name');

    startdata = data1;
    starttextnames = textnames;
    startlen = len;


    % select the file containing all of the dependant variables
    [FileName,PathName,FilterIndex] = uigetfile('C:\...','select the file with containing the dependent variables');
    Yfile = strcat(PathName,FileName);
    %Yfile = 'C:\Users\Matt\Dropbox\Research\LUR\Matlab\Datasets\A_Annual_Hourly_Only.txt';
    
    %start loop for each Y value (depedant variables).
    % load Y data (measured NO2 values)
    Ytemp = importdata(Yfile);
    AllYtemp = Ytemp;
    %Allstartdata=startdata;
    [t,ystop]=size(Ytemp.data);
    
    if randomtest
        %run random test 100 times
        for i = 1:500
            randomt(Ytemp,AllYtemp,data1,t,ystop,startdata,starttextnames,startlen,perremove,correlmatrix,onefield,randomtest);  
        end
    else
       main(startdata,starttextnames,startlen,ystop,Ytemp,correlmatrix,onefield,randomtest); 
    end
end



function randomt(Ytemp,AllYtemp,data1,t,ystop,startdata,starttextnames,startlen,perremove,correlmatrix,onefield,randomtest)
%%%%%%%%%%%%%%%% test a random sample of data%%%%%%%%%%%%%%%%%%%
% remove a random set of stations from both the indepedant dataset (data)
% and the depedant varaible dataset (Ytemp).  Creates four new variables.
% Ytemp contains the independant variable values used to develop the model.
% Ytempob contains the independant variable values used to test the model.
% startdata contains the dependant variable values used to develop the model
% startdataobs contains the dependant variable values used to test the
% model.

    %randomize the array indicies
    randorder = randperm(t);
    randdatanum = int16(perremove*t);

     
    %%%%% Gather all of the station data used to develope the model%%%%%%%
    fn = fieldnames(data1);
    for str = fn'
        [e,v]=size(data1.(char(str)));
        temp = zeros(randdatanum,v);
        %temp = data1.(char(str));
        for q = 1:randdatanum
            if randorder(q) > length(data1.(char(str)))
                temp(q,:) = NaN(1,1);
            else
                temp(q,:) = data1.(char(str))(randorder(q),:);       
            end
        end
        startdata.(char(str))=temp;
    end
    temp = zeros(randdatanum,ystop);
    for q = 1:randdatanum
       temp(q,:) = AllYtemp.data(randorder(q),:);       
    end
    Ytemp.data=temp(1:randdatanum,:);

    %Ytempob = Ytemp;
    %startdataobs=startdata;

    %%%%% Gather all of the station data to validate the model %%%%%
    fn = fieldnames(startdata);  
    for str = fn'
        [e,v]=size(data1.(char(str)));
        temp = zeros(t-randdatanum,v);
        q5 = 1;
        %temp = startdata.(char(str));
        for q = randdatanum+1:t
            if randorder(q) > length(data1.(char(str)))
                temp(q5,:) = NaN(1,1);
            else
                temp(q5,:) = data1.(char(str))(randorder(q),:);       
            end
            q5 = q5 + 1;
        end
        startdataobs.(char(str))=temp;
    end
    q5 = 1;
    temp = zeros(t-randdatanum,ystop);
    for q = randdatanum+1:t
       temp(q5,:) = AllYtemp.data(randorder(q),:); 
       q5 = q5 + 1;
    end
    Ytempob.data=temp(1:(t-randdatanum),:);
    ObsModel=zeros(length(Ytempob.data),20);
    clear temp AllYtmp Allstartdata data1 d p q
    main(startdata,starttextnames,startlen,ystop,Ytemp,correlmatrix,onefield,randomtest,ObsModel,Ytempob,startdataobs)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end


function main(startdata,starttextnames,startlen,ystop,Ytemp,correlmatrix,onefield,randomtest,ObsModel,Ytempob,startdataobs)
    % loop through each dependant variable
    for h = 1:ystop
        clear data2 origdata Y Model
        data2 = startdata;
        textnames = starttextnames;
        len = startlen;

        origdata=Ytemp.data(:,h);
        Y = origdata;
        modelnum = 0;
        Model = [];

        % start of loop for each regression iteration
        modeldone=0;
        while modeldone==0
            % combine all independant variables together for analysis

            X = [];
            X(:,1) = Y;
            fn = fieldnames(data2);
             for str = fn'
                 temp = data2.(char(str));

%                 if length(temp) < length(X)
%                     temp(length(temp)+1:length(X))=NaN(1,length(X)-length(temp));
%                     data2.(char(str))=temp;
%                 end
%                 X = [X,temp(1:length(X),:)];

                tempsize = size(temp);
                Xsize = size(X);
%                 if tempsize(1) < Xsize(1)
%                     temp(:,tempsize(2)+1:Xsize(2))=NaN(tempsize(1),Xsize(2)-tempsize(2));
%                     data2.(char(str))=temp;
%                 end               
                X = [X,temp(1:Xsize(1),:)];

             end

            %find the correlation coefficient matrix
            R = corrcoef(X,'rows','complete');

            %%%%%%%%%%%%%%%%%%%%%% output the correlation matrix for all
            %%%%%%%%%%%%%%%%%%%%%% paramters and output them to a file
            if correlmatrix
                se=length(R);
                fid = fopen('C:\...\Correlation_Matrix.txt', 'w');
                for s = 2:se
                    for p = 1:s-1
                        %only output value above 0.95
                        if (R(s,p)>.95) && (s~=p)

                            fieldnum = 0;
                            tot = 0;
                            while tot < s-1
                                fieldnum = fieldnum + 1;
                                tot = tot + len(fieldnum);
                            end
                            pos1 = s - 1 - (tot-len(fieldnum));
                            fieldstr1 = char(fn(fieldnum));
                            text1 = textnames.(fieldstr1)(pos1);

                            fieldnum = 0;
                            tot = 0;
                            while tot < p-1
                                fieldnum = fieldnum + 1;
                                tot = tot + len(fieldnum);
                            end
                            pos2 = p - 1 - (tot-len(fieldnum));
                            fieldstr2 = char(fn(fieldnum));
                            text2 = textnames.(fieldstr2)(pos2);

                            fprintf(fid, '\r\n%s\t\t %s\t\t %5.4f\t',textnames.(fieldstr1){pos1},textnames.(fieldstr2){pos2},R(s,p));
                        end
                    end
                end
                fclose(fid);
                correlmatrix = 0;
            end
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            %only look at the correlation coef between dependant var and
            %independant vars.  Find the highest value and input the data into
            %the model array.
            R1 = R(2:length(R),1);
            [C,I]=max(abs(R1));



            modelnum = modelnum + 1;
            Model.data(:,modelnum)=X(:,I+1);


            %find position of  the variable in the structures textnames and data
            %model is used to store information about each model parameter
            fieldnum = 0;
            tot = 0;
            while tot < I
                fieldnum = fieldnum + 1;
                tot = tot + len(fieldnum);
            end
            pos = I - (tot-len(fieldnum));
            fieldstr = char(fn(fieldnum));
            Model.name(modelnum)=textnames.(fieldstr)(pos);
            Model.corr(modelnum)=C;

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            if randomtest
                ObsModel(:,modelnum)=startdataobs.(fieldstr)(:,pos);
            end
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


            %remove the variable from the existing dataset for next iteration
            if onefield || (len(fieldnum)==1)
                data2 = rmfield(data2, fieldstr);
                textnames = rmfield(textnames, fieldstr);
                len(fieldnum) = [];
            else
                data2.(fieldstr)(:,pos)=[];
                textnames.(fieldstr)(pos)=[];
                len(fieldnum) = len(fieldnum)-1;
            end

            % run the regression analysis on the existing model.
            stats = regstats(origdata,Model.data,'linear',{'yhat','tstat','fstat','r','rsquare','adjrsquare','covb'});


            if modelnum > 1
                statsVIF = regstats(Model.data(:,modelnum),Model.data(:,1:(modelnum-1)),'linear',{'rsquare'});
                VIF = 1/(1-statsVIF.rsquare);
            else
                VIF = 0;
            end
                      
            % test to see if the last variable added to the model is significant and that it 
            % passes the multicolinearity test with the other parameters.  
            % If it is not the model analysis is complete and the last independant
            % variable is removed.  If it is significant statistics are
            % calculated
            if (max(stats.tstat.pval)>0.05) || (VIF > 3.0)
                modeldone=1;
                Model.data(:,modelnum)=[];
                Model.name(modelnum)=[];
                Model.corr(modelnum)=[];
                modelnum=modelnum-1;
            else
%                  varname = char(Model.name(modelnum));
%                  if strcmp(varname(1:4),'Elev') && (stats.tstat.beta(modelnum+1)>0)
%                      Model.data(:,modelnum)=[];
%                      Model.name(modelnum)=[];
%                      Model.corr(modelnum)=[];
%                      modelnum=modelnum-1;
%                  else
                    Model.rsquare(modelnum)=stats.rsquare;
                    Model.adjrsquare(modelnum)=stats.adjrsquare;
                    Y = stats.r;
                    Model.iqr(modelnum) = iqr(Model.data(:,modelnum));
                    Model.bxiqr(modelnum) = Model.iqr(modelnum) * stats.tstat.beta(modelnum+1);
%                  end
            end
        end

        % run the final regression on the final Model parameters. Output data
        % to textfiles.

        stats = regstats(origdata,Model.data,'linear',{'yhat','tstat','fstat','r','rsquare','adjrsquare','covb','mse'});
        stats.meana = mean(origdata(~isnan(origdata)));
        stats.meanb = mean(stats.yhat(~isnan(stats.yhat)));
        stats.stdeva = std(origdata(~isnan(origdata)));
        stats.stdevb = std(stats.yhat(~isnan(stats.yhat)));
        
        error = stats.yhat - origdata;
        averror = mean(error(~isnan(error)));
        avaberror = mean(abs(error(~isnan(error))));
        bias = (stats.yhat-origdata)./origdata;
        avbias = mean(bias(~isnan(bias)));
        avabbias = mean(abs(bias(~isnan(bias))));

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % calculate the Model test values and compare to the observed data for
        % model evaluation.
        if randomtest
           Yfit = zeros(length(Ytempob.data),1); 
           for prand = 1:length(Ytempob.data)
                Yfit(prand) = stats.tstat.beta(1);
                for qrand = 1:modelnum
                     Yfit(prand) = Yfit(prand) + stats.tstat.beta(qrand+1)*ObsModel(prand,qrand);   
                end
           end
           tempYtest = Ytempob.data(:,h);
           statstest = regstats(Ytempob.data(:,h),Yfit,'linear',{'rsquare','adjrsquare','mse'});
           statstest.meana = mean(tempYtest(~isnan(tempYtest)));
           statstest.meanb = mean(Yfit(~isnan(Yfit)));
           statstest.stdeva = std(tempYtest(~isnan(tempYtest)));
           statstest.stdevb = std(Yfit(~isnan(Yfit)));
           rerror = Yfit - tempYtest;
           raverror = mean(rerror(~isnan(rerror)));
           ravaberror = mean(abs(rerror(~isnan(rerror))));
           rbias = (Yfit - tempYtest)./tempYtest;
           ravbias = mean(rbias(~isnan(rbias)));
           ravabbias = mean(abs(rbias(~isnan(rbias))));
           %figure(h)
           %hold on
           %scatter(Ytempob.data(:,h),Yfit,'+')
           %scatter(Ytemp.data(:,h),stats.yhat,'o')
           %hold off



        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


        fid = fopen('C:\...\output.txt', 'a');
        %Yname = cell2str(Ytemp.textdata(1));
        fprintf(fid, '\r\n%s\t %5.4f\t %5.4f\t %5.6f\t %5.6f\t', Ytemp.textdata{h},stats.rsquare,stats.adjrsquare,stats.fstat.sse,stats.mse);
        if randomtest
            %stats. structure contains model building observations,
            %statstest. structure contains model testing observations
            fprintf(fid, '%5.4f\t %5.4f\t %5.4f\t %5.4f\t %5.4f\t %5.4f\t %5.4f\t %5.4f\t %5.4f\t %5.4f\t %5.4f\t %5.4f\t %5.4f\t %5.4f\t %5.4f\t %5.4f\t %5.4f\t %5.4f\t %5.4f\t', stats.meana, stats.meanb, stats.stdeva, stats.stdevb, statstest.rsquare, statstest.adjrsquare, statstest.mse, statstest.meana, statstest.meanb, statstest.stdeva, statstest.stdevb,averror,avaberror,avbias,avabbias,raverror,ravaberror,ravbias,ravabbias);
        end
        for k = 1:modelnum
            fprintf(fid, '%s\t', Model.name{k});
        end

        fclose(fid);
        fid = fopen('C:\...\outputparameters.txt', 'a');
        fprintf(fid, '\r\n\r\n\r\n%s\t %s','Model Name',Ytemp.textdata{h});
        fprintf(fid,'\r\n%s\t %s\t %s\t %s\t %s\t %s\t %s\t %s\t', 'Rsqr', 'Adj_Rsqr', 'N', 'sse', 'dfr', 'ssr', 'f', 'pval');
        fprintf(fid,'\r\n%5.4f\t %5.4f\t %i\t %8.2f\t %i\t %8.2f\t %7.4f\t %5.4f\t', stats.rsquare, stats.adjrsquare, stats.fstat.dfe, stats.fstat.sse,stats.fstat.dfr,stats.fstat.ssr,stats.fstat.f,stats.fstat.pval);
        fprintf(fid,'\r\n%s\t %s\t %s\t %s\t', 'error', 'Ab_error', 'bias', 'ab_bias');
        fprintf(fid,'\r\n%8.4f\t %8.4f\t %8.4f\t %8.4f\t', averror, avaberror, avbias, avabbias);
 
        fprintf(fid,'\r\n\r\n%s\t %s\t %s\t %s\t %s\t %s\t %s\t %s\t %s\t %s\t', 'Parameter','B Coef','Stand Error','T stat','pval','rsqr','adjrsqr','Correl','iqr','Beta x iqr','VIF');
        fprintf(fid,'\r\n%s\t %12.11f\t %12.11f\t %5.4f\t %5.4f\t', 'Intercept',stats.tstat.beta(1),stats.tstat.se(1),stats.tstat.t(1),stats.tstat.pval(1));
        for k = 1:modelnum
            fprintf(fid, '\r\n%s\t %12.11f\t %12.11f\t %3.2f\t %3.2f\t %5.4f\t %5.4f\t %5.4f\t %12.11f\t %5.4f\t', Model.name{k},stats.tstat.beta(k+1),stats.tstat.se(k+1),stats.tstat.t(k+1),stats.tstat.pval(k+1),Model.rsquare(k),Model.adjrsquare(k),Model.corr(k),Model.iqr(k),Model.bxiqr(k));
            para = Model.data(:,k);
            otherpara = Model.data;
            otherpara(:,k)=[];
            statspara = regstats(para,otherpara,'linear',{'rsquare'});
            VIF = 1/(1-statspara.rsquare);
            fprintf(fid, '%5.4f\t',VIF);
        end
        fclose(fid);



    end
    
    
end
%test to see if the model is finished tstat..
%end loop 

% data.field1(:,3)=[] removes column of data from field
% s = rmfield(s, 'fieldname') remove the entire field



