3.10 In this question we look at the logit model for quantal assay data. At each of a number of doses, xi, groups of ni individuals are given the dose, and ri respond, 1 ≤ i ≤ k. In the logit model, the probability of response at dose xi is written as P(xi)=1/(1 + e−(α+βxi)

)

(i) Study the following Matlab program for evaluating the negative loglikelihood.

function y=logit(t)

%

%LOGIT Calculates the negative log-likelihood for a logistic

% model and quantal response data

% x is the ‘dose’

% n is the number of individuals at each ‘dose’

% r is the number of individuals that respond at each ‘dose’

%_____________________________________________________________ x=[49.06 52.99 56.91 60.84 64.76 68.69 72.61 76.54];

n=[59 60 62 56 63 59 62 60];

r=[6 13 18 28 52 53 61 60];

w=ones(size(x));

alpha=t(1); beta=t(2);

y=w ./(1+exp(-(alpha+beta*x)));

z=w-y;

loglik=r*(log(y))’+(n-r)*(log(z))’;

y=-loglik;

(ii) Obtain the maximum likelihood estimate of the pair of parameters (α, β), using the Matlab fminsearch routine.

(iii) Display the data, by plotting ri/ni vs xi, 1 ≤ i ≤ k.

(iv) Retain the plot on the screen by using the command hold, and then plot the fitted logistic curve.