Logistic Model Generalization

Logistic Model Generalization

Consider the epilepsy data given in HSAUR2 package of R.

For this epilepsy data, define a response Y as follows: Y = 1 ifthe seizure rate > 4 and 0 otherwise. Now, fit the generalized linearmixed effect model to Y using the logistic model under the binomialfamily conditional on the random intercept term. The conditionallogistic model should include baseline seizure count, treatment,log(age), and period (as a numeric variable) as the fixed effectexplanatory variables.
(a) Interpret the regression coefficient for the treatment variable.
(b) Predict the probability of Y = 1 for subject 29 at period= 4. Notethat you have to predict conditional on the random effect.
(c) Estimate the marginal probability of Y = 1 at period= 4 for thesubject that would have the same value for the explanatory variablesas of subject 29? 

Solution 

mixedlogistic.R

library(HSAUR2)

x <- epilepsy

x$seizure.rate<-ifelse(x$seizure.rate>4,1,0)

x$age<-log(x$age)

colnames(x) <-c(“treatment”,”base”,”logage”,”response” ,”period”,”subject”)

#we’ll be taking subject as the random effect to use mixed model logistic regression

library(lme4)

model<- glmer(x$response~x$treatment+x$base+x$logage+x$period+(1 | subject),data=x , family=”binomial”)

#this is our reqd. model

summary(model)

#look at p-value of x$treatmentProgabide , the p-value is 0.0757 , implying that the effect is not too significant.

pred.R 

library(HSAUR2)

x <- epilepsy

x$seizure.rate<-ifelse(x$seizure.rate>4,1,0)

x$age<-log(x$age)

x$treatment<-as.factor(x$treatment)

colnames(x) <-c(“treatment”,”base”,”logage”,”response” ,”period”,”subject”)

#we’ll be taking subject as the random effect to use mixed model logistic regression

library(lme4)

model<- glmer(x$response~x$treatment+x$base+x$logage+x$period+(1 | subject),data=x , family=”binomial”)

#this is our reqd. model

summary(model)

#look at p-value of x$treatmentProgabide , the p-value is 0.0757 , implying that the effect is not too significant.

pred<- predict(model, newdata=epilepsy)

prob<- exp(pred)/(1+exp(pred)) #as P(Y=1 given predictors) =e^z/(1+e^z)

#look at the id 293 at the leftmost coloumn in epilepsy , that’s the 29th subject at 4th period

#so we look at the 293 marked element of prob , that’s our answer for part 2 , i.e. answer is 0.99592817

 

#now the 3rd part requires the marginal probability , that’s mathematically equivalent to running an ordinary logistic regression

#and finding the probability P(Y=1) for subject 29 at period 4

model2 <- glm(x$response ~.,family=binomial(link=’logit’),data=x) #ignore the warnings , since some values are actually very close to 0 or 1

fitted.results<- predict(model2,newdata=x,type=’response’)#ignore the warning , the warning appears since this fit is rank defficient , but it’s computationally correct

print(fitted.results)

#again look at id 293 of fiited.results. Our answer is 1.000000e+00 , which actually means that the marginal probability is very close to 1