Step by step procedure on how to run nested logistic regression in R

Question

Please provide me with detailed (as possible) steps on how to do nested logistic regression in R. I'm new to R so it would help me a lot if i can get a detailed answer.

We tested how fisher's decision to exit the fishery is affected by different socioeconomic factors. Dependent variable: 0 - stay; 1 - exit Predictors: Age (continuous); Education (categorical); number of children (continuous), etc.

Our respondents were from different towns. A reviewer of our paper instructed us to account for the town where the respondent comes from by using nested logistic regression.

Your help is very much appreciated. Many thanks.

Answer 1

install.packages("mlogit")    
library(mlogit)

my.data <- YOUR.DATA    

nested.logit <- mlogit(stay.exit~ age + education + children , my.data,
shape='long', alt.var='town.list', nests=list(town.list))

See page 19 of the mlogit manual for examples of nested logit model calls. You'll have to review the documentation yourself to ensure you're getting at what you need in terms of options. http://cran.r-project.org/web/packages/mlogit/mlogit.pdf

Segue: I usually like to have a peek at all models by town.list before I look at nested models:

Note: if your categorical variables are not factored you'll have to surround them with as.factor(variable) in your model formula

# Show a little love for plyr
library(plyr)

## RNG
set.seed(123454321)

## Create a list object to store your models
my.models <- list()

## import your data
my.data <- YOUR.DATA

## Create a loop that runs by the list of towns
for(x in 1:length(mydata$town.list) {
## subset data in each step by the town
dat <- subset(my.data, town == town.list[x])
## Save the model to it's own place in the list, identified by town
my.models[[town.list[x]]] <- glm(formula = stay.exit ~ age + education + children, 
family = binomial(link = "logit"), 
data=dat)
}

## View summaries for all models
llply(my.models, summary)

## Access specific models
my.models$<TOWN NAME>

Answer 2

If I understand correctly, you want a varying-intercept model by town, i.e., a hierarchical model? If that's right, just use lme4 package.

Here is an example. Assuming you have in your data frame a variable (factor) called town, and that your data frame is called "fish", just run:

library(lme4)
library(arm) # to use the function display, much better than summary
nest.reg <- glmer(decision ~ age + education + children + (1|town), family = binomial, data = fish)
coef(nest.reg) # this will give the estimated coeficients by town (in this case, only the intercepts will vary).
fixef(nest.reg) # this will give the model averaging over all towns.
ranef(nest.reg) # the errors (specificity) at the town level. If you sum fixef with ranef you will get coef results

Finnaly, it's important to compare the estimated within and between town variation

display(nest.reg) # this will show you, among other things, the estimated residual variatio at the town and individual level. It's the error terms by town and by individual (Residual). The ratio between them will tell you how much information the is between town and within each town. 

Take a look at the last edition of Gelman and Hill's book for more information about multilevel regression using lme4.

Ps.: It's possible to include varying-slope by town as well. If that's what you need, just ask in the comments.