Построить две кривые в логистической регрессии в R

Question

Я использую логистическую регрессию в R (glm). Затем мне удается построить результат. Мой код выглядит следующим образом:

 temperature.glm = glm(Response~Temperature, data=mydata,family=binomial)

 plot(mydata$Temperature,mydata$Response, ,xlab="Temperature",ylab="Probability of Response")
 curve(predict(temperature.glm,data.frame(Temperature=x),type="resp"),add=TRUE, col="red")
 points(mydata$Temperature,fitted(temperature.glm),pch=20)
 title(main="Response-Temperature with Fitted GLM Logistic Regression Line") 

Мои вопросы:

1. Как можно построить две кривые логистической регрессии на одном графике?
2. Я получил эти два коэффициента из другого программного обеспечения для статистики. Как я могу создать случайные данные, подключить эти два набора коэффициентов (Набор 1 и Набор 2), а затем получить две кривые логистической регрессии?

Модели:

                   SET 1
 (Intercept)     -88.4505
 Temperature       2.9677

                  SET 2
 (Intercept)    -88.585533
 Temperature      2.972168

mydata в 2 столбцах и ~ 700 строк.

Response Temperature 
1 29.33 
1 30.37 
1 29.52 
1 29.66 
1 29.57 
1 30.04 
1 30.58 
1 30.41 
1 29.61 
1 30.51 
1 30.91 
1 30.74 
1 29.91 
1 29.99 
1 29.99 
1 29.99 
1 29.99 
1 29.99 
1 29.99 
1 30.71 
0 29.56 
0 29.56 
0 29.56 
0 29.56 
0 29.56 
0 29.57 
0 29.51

Answer 1

1. Чтобы построить кривую, вам просто нужно определить связь между откликом и предиктором и указать диапазон значений предиктора, для которых вы хотите построить эту кривую.Например:

   dat <- structure(list(Response = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
     1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 
     0L, 0L), Temperature = c(29.33, 30.37, 29.52, 29.66, 29.57, 30.04, 
     30.58, 30.41, 29.61, 30.51, 30.91, 30.74, 29.91, 29.99, 29.99, 
     29.99, 29.99, 29.99, 29.99, 30.71, 29.56, 29.56, 29.56, 29.56, 
     29.56, 29.57, 29.51)), .Names = c("Response", "Temperature"), 
     class = "data.frame", row.names = c(NA, -27L))
   
   temperature.glm <- glm(Response ~ Temperature, data=dat, family=binomial)
   
   plot(dat$Temperature, dat$Response, xlab="Temperature", 
        ylab="Probability of Response")
   curve(predict(temperature.glm, data.frame(Temperature=x), type="resp"), 
         add=TRUE, col="red")
   # To add an additional curve, e.g. that which corresponds to 'Set 1':
   curve(plogis(-88.4505 + 2.9677*x), min(dat$Temperature), 
         max(dat$Temperature), add=TRUE, lwd=2, lty=3)
   legend('bottomright', c('temp.glm', 'Set 1'), lty=c(1, 3), 
          col=2:1, lwd=1:2, bty='n', cex=0.8)
   

   Во втором вызове curve выше мы говорим, что логистическая функция определяет отношения между x и y.Результат plogis(z) эквивалентен результату, полученному при оценке 1/(1+exp(-z)).Аргументы min(dat$Temperature) и max(dat$Temperature) определяют диапазон x, для которого следует оценить y.Нам не нужно указывать функции, что x относится к температуре;это неявно, когда мы указываем, что ответ должен оцениваться для этого диапазона значений предикторов.

   

2. Как вы можете видеть, функция curveпозволяет построить кривую без необходимости моделирования данных предиктора (например, температуры).Если вам все еще нужно сделать это, например, чтобы построить некоторые смоделированные результаты испытаний Бернулли, которые соответствуют определенной модели, то вы можете попробовать следующее:

   n <- 100 # size of random sample
   
   # generate random temperature data (n draws, uniform b/w 27 and 33)
   temp <- runif(n, 27, 33)
   
   # Define a function to perform a Bernoulli trial for each value of temp, 
   #   with probability of success for each trial determined by the logistic
   #   model with intercept = alpha and coef for temperature = beta.
   # The function also plots the outcomes of these Bernoulli trials against the 
   #   random temp data, and overlays the curve that corresponds to the model
   #   used to simulate the response data.
   sim.response <- function(alpha, beta) {
     y <- sapply(temp, function(x) rbinom(1, 1, plogis(alpha + beta*x)))  
     plot(y ~ temp, pch=20, xlab='Temperature', ylab='Response')
     curve(plogis(alpha + beta*x), min(temp), max(temp), add=TRUE, lwd=2)    
     return(y)
   }
   

   Примеры:

   # Simulate response data for your model 'Set 1'
   y <- sim.response(-88.4505, 2.9677)
   
   # Simulate response data for your model 'Set 2'
   y <- sim.response(-88.585533, 2.972168)
   
   # Simulate response data for your model temperature.glm
   # Here, coef(temperature.glm)[1] and coef(temperature.glm)[2] refer to
   #   the intercept and slope, respectively
   y <- sim.response(coef(temperature.glm)[1], coef(temperature.glm)[2])
   

   На рисунке ниже показан график, полученный в первом примере выше, то есть результаты одного испытания Бернулли для каждого значения случайного вектора температуры, а также кривая, описывающая модель, из которой были смоделированы данные.