Looking for some help fitting a binomial model with a fixed intercept. I've seen this question covered successfully with linear models by subtracting the explicit intercept from the regression and and then fit the intercept-free model:
intercept <- 1.0
fit <- lm(I(x - intercept) ~ 0 + y, lin)
summary(fit)
Or by specifying an offset:
m(x ~ y +0 + offset(rep(1, nrow(lin))),
But I can't figure out how to do this with a binomial regression. Example code below:
library(ggplot2)
library(lme4)
library(dplyr)
library(gather)
# simulate dummy dataframe
x <- data.frame(time = c(1, 1, 1, 1, 1, 1,1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2,3, 3, 3, 3, 3, 3, 3, 3, 3,4, 4, 4, 4, 4, 4, 4, 4, 4),
type = c('a', 'a', 'a', 'b', 'b', 'b','c','c','c','a', 'a', 'a', 'b', 'b', 'b','c','c','c','a', 'a', 'a', 'b', 'b', 'b','c','c','c','a', 'a', 'a', 'b', 'b', 'b','c','c','c'), randef=c('aa','ab','ac','ba','bb','bc','ca','cb','cc','aa','ab','ac','ba','bb','bc','ca','cb','cc','aa','ab','ac','ba','bb','bc','ca','cb','cc','aa','ab','ac','ba','bb','bc','ca','cb','cc'),
surv = sample(x = 1:200, size = 36, replace = TRUE),
nonsurv= sample(x = 1:200, size = 36, replace = TRUE))
# convert to survival and non survival into individuals following
# https://stackoverflow.com/questions/51519690/convert-cbind-format-for- binomial-glm-in-r-to-a-dataframe-with-individual-rows
x_long <- x %>%
gather(code, count, surv, nonsurv)
# function to repeat a data.frame
x_df <- function(x, n){
do.call('rbind', replicate(n, x, simplify = FALSE))
}
# loop through each row and then rbind together
x_full <- do.call('rbind',
lapply(1:nrow(x_long),
FUN = function(i)
x_df(x_long[i,], x_long[i,]$count)))
# create binary_code
x_full$binary <- as.numeric(x_full$code == 'surv')
### binomial glm with interaction between time and type:
summary(fm2<-glm(binary ~ time*type, data = x_full, family = "binomial"))
### plot glm in ggplot2
ggplot(x_full, aes(x = time, y = as.numeric(x_full$binary), fill= x_full$type)) +
geom_smooth(method="glm", aes(color = factor(x_full$type)), method.args = list(family = "binomial"))
What i'd like to do is force the intercept through 1 (i.e. all surviving in the first time point), almost like a survival curve. Specifying an offset doesn't seem to work in this example, as it is binomial probability, so specifying a 1 or a 0 doesn't work, and the number of surviving vs non-surviving individuals in the first time point varies among factors).
Stack Overflow, i've learnt more from you than any of my undergraduate classes or postgraduate supervisors. Any suggestions on this one?
