Question

How can I plot the logistic regression? I would like to plot the dependent variable on the y-axis and independent on the x. I called the coefficients and got an output, so no errors on the script.

Here's the data for the independent variable (SupPres):

#Set the range for water supply pressure
SupPres <- c(20:120)
#Create a normal distribution for water supply pressure
SupPres <- rnorm(3000, mean=70, sd=25)   

Logistic regression and creating y-variable:

#Create logistic regression
z=1+2*NozHosUn+3*SupPres+4*PlaceSet+5*Hrs4+6*WatTemp
z <- (z-mean(z))/sd(z)
pr = 1/(1+exp(-z))
y <- rbinom(3000,1,pr)
DishWa=data.frame(y=y, NozHosUn=NozHosUn,SupPres=SupPres,
               PlaceSet=PlaceSet,Hrs4=Hrs4,WatTemp=WatTemp)
glm(y~NozHosUn+SupPres+PlaceSet+Hrs4+WatTemp, data=DishWa, 
    family=binomial)

Please let me know if I can provide more information. Thanks.

Was it helpful?

Solution

A possible solution that takes up a suggestion in Michael J. Crawley's excellent book "Statistics: an introduction using R" is the following code:

attach(mtcars);
model <- glm(formula = am ~ hp + wt, family = binomial);
print(summary(model));
model.h <- glm(formula = am ~ hp, family = binomial);
model.w <- glm(formula = am ~ wt, family = binomial);
op <- par(mfrow = c(1,2));
xv <- seq(0, 350, 1);
yv <- predict(model.h, list(hp = xv), type = "response");
hp.intervals <- cut(hp, 3);
plot(hp, am);
lines(xv, yv);
points(hp,fitted(model.h),pch=20);
am.mean.proportion <- tapply(am, hp.intervals, sum)[[2]] / table(hp.intervals)[[2]];
am.mean.proportion.sd <- sqrt(am.mean.proportion * abs(tapply(am, hp.intervals, sum)[[3]] - tapply(am, hp.intervals, sum)[[1]]) / table(hp.intervals)[[2]]);
points(median(hp), am.mean.proportion, pch = 16);
lines(c(median(hp), median(hp)), c(am.mean.proportion - am.mean.proportion.sd, am.mean.proportion + am.mean.proportion.sd));
xv <- seq(0, 6, 0.01);
yv <- predict(model.w, list(wt = xv), type = "response");
wt.intervals <- cut(wt, 3);
plot(wt, am);
lines(xv, yv);
points(wt,fitted(model.w),pch=20);
am.mean.proportion <- tapply(am, wt.intervals, sum)[[2]] / table(wt.intervals)[[2]];
am.mean.proportion.sd <- sqrt(am.mean.proportion * abs(tapply(am, wt.intervals, sum)[[3]] - tapply(am, wt.intervals, sum)[[1]]) / table(wt.intervals)[[2]]);
points(median(wt), am.mean.proportion, pch = 16);
lines(c(median(wt), median(wt)), c(am.mean.proportion - am.mean.proportion.sd, am.mean.proportion + am.mean.proportion.sd));
detach(mtcars);
par(op);

In addition to the standard glm plot it contains an indicator for the fit of the central third to the data. It shows that the fit to hp is poor while that to wt is rather good.

Example output

The example uses R's built-in "cars" dataset.

OTHER TIPS

There are several ways to plot the results of logistic regression. See https://sites.google.com/site/daishizuka/toolkits/plotting-logistic-regression-in-r and Plot two curves in logistic regression in R for explanations and examples in R.

Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top