%% Construct bootstrap samples for a nonlinear regression and estimate %% bootstrap confidence intervals. %% Fit simulated dose-response data to a Hill function, given by %% y(x) = ymax/(1 + (EC50/dose)^n) %% Load data DoseResponse = load('DoseResponseData.dat', '-ascii'); dose = DoseResponse(:,1); yExpt = DoseResponse(:,2); %% Call nlinfit to perform nonlinear regression. betaGuess = [1, 1, 1]; % Initial guess for parameter values [betaHat, residuals, Jac, Cov, mse] = nlinfit(dose, yExpt, @Hill, betaGuess); %% Resample the residuals using the bootstrp function nboot = 200; % Number of bootstrap replicates [bootresiduals, bootIndices] = bootstrp(nboot, [], residuals); bootIndices(:,1:3) % Print the first 3 columns of bootIndices bootResiduals = residuals(bootIndices); bootResiduals(:, 1:3) % Print the first 3 columns of bootResiduals %% Generate bootstrap datasets by adding the resampled residuals %%to the best fit curve yMod = Hill(betaHat, dose); yBoot = repmat(yMod, 1, nboot) + bootResiduals; %% Fit each of the bootstrap datasets to the Hill function to %% build up the bootstrap distributions of parameter estimates betaBoot = zeros(nboot, 3); for i=1:nboot betaBoot(i,:) = nlinfit(dose, yBoot(:,i), @Hill, betaGuess); end % Estimate 95% CIs bootCI = prctile(betaBoot,[2.5 97.5]); betaHat bootCI %% Plot histograms of bootstrap distributions for i = 1:3 subplot(3,1,i) hist(betaBoot(:,i), 30) h = findobj(gca,'Type','patch'); set(h, 'EdgeColor', 'w') end axes(subplot(311)) xlabel('y_{max}') axes(subplot(312)) xlabel('log_{10}(EC_{50})') axes(subplot(313)) xlabel('Hill coefficient, n')