0
$\begingroup$

I would like to generate and plot the 95% confidence intervals for the fitted logistics model code below, with filling between the mean curve and its upper and lower bands. Appreciate suggestions!!

Clear[data, se, a, b]
data = {{1, 1}, {2, 0.75}, {3, 0.611111111}, {4, 0.611111111}, {5, 
    0.611111111}, {6, 0.416666667}, {7, 0.416666667}, {8, 
    0.666666667}, {9, 0.535353535}, {10, 0.45508658}, {11, 
    0.387598676}, {12, 0.340230255}, {13, 0.215230255}, {14, 
    0.188271007}, {15, 0.187838107}, {16, 0.224308695}, {17, 
    0.229369424}, {18, 0.235463951}, {19, 0.2488416}, {20, 
    0.248417357}, {21, 0.206599175}, {22, 0.174664443}, {23, 
    0.171198551}, {24, 0.163201729}, {25, 0.150447012}, {26, 
    0.164363959}, {27, 0.166265958}, {28, 0.140662245}, {29, 
    0.128777833}, {30, 0.129944921}, {31, 0.123349511}, {32, 
    0.111959778}, {33, 0.111616813}, {34, 0.102635032}, {35, 
    0.086485736}, {36, 0.075479421}, {37, 0.07383027}, {38, 
    0.081188587}, {39, 0.090044466}, {40, 0.084028814}, {41, 
    0.096462466}, {42, 0.087077379}};
Length[data]
lm = LogitModelFit[data, x, x];
Normal[lm]
lm["ParameterTable"]
lm["ParameterConfidenceIntervalTable", ConfidenceLevel -> 0.95]
Plot[lm[x], {x, 1, 60}, Epilog :> Point[data], 
 PlotStyle -> {Red, Thin}]
$\endgroup$
4
  • $\begingroup$ You don't have appropriate data to use LogitModelFit. For the response variable one needs either counts or the proportion of successes and the number of trials for each response. In addition, it looks like the responses are highly correlated (i.e., you don't have independent observations). What is the data and how was it collected? (I'm not looking for subject matter details. Just need to determine mainly if the proportions are a ratio of integers or if the proportion is a proportion volume or area or weight or some continuous variable.) $\endgroup$
    – JimB
    Apr 21, 2020 at 2:20
  • $\begingroup$ Ooops! I need to learn to read the documentation better and not assume that Mathematica does it the same way that other software does it. One does NOT enter a count (i.e., the number of successes) but rather a proportion between zero and one. (Although, with the number of trials for each observation being 1, then the observed proportion and the count are identical.) I'll update my answer with a better example. $\endgroup$
    – JimB
    Apr 21, 2020 at 15:10
  • $\begingroup$ How about just use my dataset? $\endgroup$
    – user42700
    Apr 21, 2020 at 16:12
  • 1
    $\begingroup$ You dataset isn't currently complete (or maybe appropriate) yet. If indeed your proportions are based on counts, then your sample size (i.e., number of trials) must be huge for each of the 42 observations. So first of all you have not included the sample size for each observation. Why is that critical? While it doesn't make much difference for the estimates of the parameters, it is critical to figure out appropriate measures of precision (confidence bands, confidence intervals, standard errors, etc.) And those measures of precision require a binomial distribution. And I'm not sure you do. $\endgroup$
    – JimB
    Apr 21, 2020 at 16:29

2 Answers 2

2
$\begingroup$

If the data is not thought to follow a binomial distribution but rather the data consists of proportions associated with a proportion of weight, height, area, volume, etc., then LogitModelFit is not appropriate.

However, one can fit a logistic model using LinearModelFit or NonlinearModelFit and obtain confidence bands directly by transforming the response variable. (For the OP's data this will ignore the assumption of independent observations as serially correlation seems pretty strong. Ignoring this also affects the quality of the confidence bands. A description as to how the data was collected would be necessary to account for the apparent serial correlation.)

data = {{1, 1}, {2, 0.75}, {3, 0.611111111}, {4, 0.611111111}, {5, 0.611111111}, 
  {6, 0.416666667}, {7, 0.416666667}, {8, 0.666666667}, {9, 0.535353535},
  {10, 0.45508658}, {11, 0.387598676}, {12, 0.340230255}, {13, 0.215230255},
  {14, 0.188271007}, {15, 0.187838107}, {16, 0.224308695}, {17, 0.229369424}, 
  {18, 0.235463951}, {19, 0.2488416}, {20, 0.248417357}, {21, 0.206599175},
  {22, 0.174664443}, {23, 0.171198551}, {24, 0.163201729}, {25, 0.150447012},
  {26, 0.164363959}, {27, 0.166265958}, {28, 0.140662245}, {29, 0.128777833}, 
  {30, 0.129944921}, {31, 0.123349511}, {32, 0.111959778}, {33, 0.111616813}, 
  {34, 0.102635032}, {35, 0.086485736}, {36, 0.075479421}, {37, 0.07383027},
  {38, 0.081188587}, {39, 0.090044466}, {40, 0.084028814}, {41, 0.096462466}, 
  {42, 0.087077379}};

(* Drop first point and convert to logits *)
logit = data[[2 ;;]];
logit[[All, 2]] = Log[logit[[All, 2]]] - Log[1 - logit[[All, 2]]];

(* Fit model *)
lm = LinearModelFit[logit, {x, x^2}, x];

(* Plot on logit scale *)
Plot[{lm[x], lm["MeanPredictionBands"], 
  lm["SinglePredictionBands"]}, {x, 1, 42}, Epilog :> Point[logit], 
 PlotStyle -> {Red, Blue, Green},
 Frame -> True, FrameLabel -> (Style[#, Bold, 18] &) /@ {"x", "Log[p/(1-p)]"},
 PlotLegends -> {"Prediction", "95% Mean prediction intervals", 
   "95% Single prediction intervals"}]

Data and predictions on logit scale

(* Plot on original scale *)
Plot[{1 - 1/(1 + Exp[lm[x]]), 
  1 - 1/(1 + Exp[lm["MeanPredictionBands"]]),
  1 - 1/(1 + Exp[lm["SinglePredictionBands"]])}, {x, 1, 42}, 
 PlotRange -> {Automatic, {0, 1}}, Epilog :> Point[data],
 PlotStyle -> {Red, Blue, Green},
 Frame -> True, FrameLabel -> (Style[#, Bold, 18] &) /@ {"x", "p"},
 PlotLegends -> {"Prediction", "95% Mean prediction intervals", 
   "95% Single prediction intervals"}]

Data and predictions on original scale

One could also fit a "logist"-shaped curve with NonlinearModelFit without having to transform the response variable.

$\endgroup$
1
  • $\begingroup$ Hi Jim: That is the route I took - fitting a "logist" shaped exponential to the data as reported ... many thanks to all for the thoughfulness of these responses ... very appreciated ... prg $\endgroup$
    – user42700
    Apr 23, 2020 at 16:50
2
$\begingroup$

The LogitModelFit function expects the response variable to be a proportion between 0 and 1 based on count data from a binomial distribution. No other kind of proportion (proportion of weight, proportion of area, proportion of volume, etc.) is appropriate.

Having said that LogitModelFit will still provide an estimate of the parameters for a logistic fit for any set of proportions. So what's the problem? Because LogitModelFit assumes that the data comes from a binomial distribution, the measures of precision (standard errors, confidence intervals, etc.) all depend on that assumption. That assumption does not hold when the proportions are ratios of non-integer numbers. (For some datasets that assumption doesn't necessarily hold even for ratios of integers.) It's not just the data values that matter. It's also how the data was sampled.

Here is an example:

(* Predictor variable *)
x = Range[42];
(* Number of trials associated with each observation *)
n = ConstantArray[10, Length[x]];
(* Number of successes *)
y = {6, 5, 7, 4, 7, 7, 8, 4, 5, 4, 5, 3, 3, 3, 0, 1, 3, 4, 2, 4, 3, 3,
    0, 0, 3, 2, 1, 2, 2, 1, 1, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0};

(* Construct data for input to LogitModelFit *)
data = Transpose[{x, y/n}]; (* Response variable needs to be a proportion between 0 and 1 *)

(* Perform analysis *)
lmf = LogitModelFit[data, z, z, Weights -> n];
lmf["ParameterTable"]

Parameter table

If Weights -> n was not given, then what happens is that the weights are all set to 1 which gives way too large standard errors of the parameters. (Note that Weights is especially necessary if the sample sizes vary for each observation.)

To find an approximate 95% confidence band one needs to extract the estimated parameters and the associated covariance matrix:

(* Get estimated parameters and the associated covariance matrix *)
coeff = lmf["BestFitParameters"];
cov = lmf["CovarianceMatrix"];

(* Get linear predictor and lower and upper confidence limits *)
xx = Min[data[[All, 1]]] + (Max[data[[All, 1]]] - 
     Min[data[[All, 1]]]) Range[0, 100]/100;
predicted = coeff.{1, #} & /@ xx;
halfWidths = 1.96 ({1, #}.cov.{1, #})^0.5 & /@ xx;
upper = predicted + halfWidths;
lower = predicted - halfWidths;

(* Convert from logit scale to probability *)
predicted = 1 - 1/(1 + Exp[predicted]);
lower = 1 - 1/(1 + Exp[lower]);
upper = 1 - 1/(1 + Exp[upper]);

(* Plot results *)
ListPlot[{predicted, lower, upper}, Joined -> True, 
 PlotStyle -> {Black, Gray, Gray}, Frame -> True]

Predicted results and approximate 95% confidence band

$\endgroup$
4
  • $\begingroup$ JimB: Thank you! ... As long as the data is between 0 and 1 (inclusive) you can fit a Logit Model too ... $\endgroup$
    – user42700
    Apr 21, 2020 at 11:36
  • $\begingroup$ Thanks. I see the error in my ways. I need to learn to read better. Yes, the response variable for LogitModelFit does need to be between 0 and 1. The example that I gave only dealt with a sample from a binomial distribution with the number of trials being just 1. I'll add in a better example. $\endgroup$
    – JimB
    Apr 21, 2020 at 15:07
  • $\begingroup$ Why not use my dataset? $\endgroup$
    – user42700
    Apr 21, 2020 at 16:36
  • $\begingroup$ "Why not use my dataset?" now answered in a comment in the original post. $\endgroup$
    – JimB
    Apr 21, 2020 at 16:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.