Hands-on with GLMs I

Stat 203 Lecture 23

Dr. Janssen

Explorations

Exploring `blocks`

Children were asked to build towers as high as they could out of cubical and cylindrical blocks. The number of blocks used and time taken were recorded and are available in the blocks dataset in the GLMsData package. For now, we’ll consider only the number of blocks \(y\) and the age of the child \(x\).

Explore the data; plot the number of blocks used against the age of the child. Propose a GLM for the data, identifying the EDM, parameters, and systematic component.

Code

library(GLMsData); data(blocks);

par(mfrow=c(1, 2))
plot(jitter(Number)~Age, data=blocks)
plot( Number~cut(Age, 3), data=blocks)

Case Study: Household Size in the Philippines

Setting the Stage

How many other people live with you in your home? The number of people sharing a house differs from country to country and often from region to region. International agencies use household size when determining needs of populations, and the household sizes determine the magnitude of the household needs.

The Philippine Statistics Authority (PSA) spearheads the Family Income and Expenditure Survey (FIES) nationwide. The survey, which is undertaken every three years, is aimed at providing data on family income and expenditure, including levels of consumption by item of expenditure. Our data, from the 2015 FIES, is a subset of 1500 of the 40,000 observations [@PSA]. Our data set focuses on five regions: Central Luzon, Metro Manila, Ilocos, Davao, and Visayas (see Figure @ref(fig:philippinesmap)).

At what age are heads of households in the Philippines most likely to find the largest number of people in their household? Is this association similar for poorer households (measured by the presence of a roof made from predominantly light/salvaged materials)? We begin by explicitly defining our response, \(Y=\) number of household members other than the head of the household. We then define the explanatory variables: age of the head of the household, type of roof (predominantly light/salvaged material or predominantly strong material), and location (Central Luzon, Davao Region, Ilocos Region, Metro Manila, or Visayas). Note that predominantly light/salvaged materials are a combination of light material, mixed but predominantly light material, and mixed but predominantly salvaged material, and salvaged matrial. Our response is a count, so we consider a Poisson regression where the parameter of interest is \(\mu\), the average number of people, other than the head, per household. We will primarily examine the relationship between household size and age of the head of household, controlling for location and income.

Importing Data

fHH1 <- read.csv("https://prof.mkjanssen.org/glm/data/fHH1.csv")
head(fHH1)

  X     location age total numLT5                          roof
1 1 CentralLuzon  65     0      0 Predominantly Strong Material
2 2  MetroManila  75     3      0 Predominantly Strong Material
3 3  DavaoRegion  54     4      0 Predominantly Strong Material
4 4      Visayas  49     3      0 Predominantly Strong Material
5 5  MetroManila  74     3      0 Predominantly Strong Material
6 6      Visayas  59     6      0 Predominantly Strong Material

location = where the house is located (Central Luzon, Davao Region, Ilocos Region, Metro Manila, or Visayas)
age = the age of the head of household
total = the number of people in the household other than the head
numLT5 = the number in the household under 5 years of age
roof = the type of roof in the household (either Predominantly Light/Salvaged Material, or Predominantly Strong Material, where stronger material can sometimes be used as a proxy for greater wealth)

Questions

Exploration

In groups of 2-3, brainstorm some questions about this data.

Exploratory Data Analysis

Code

m <- mean(fHH1$total)
v <- var(fHH1$total)
c(Mean=m,Variance=v)

    Mean Variance 
3.684667 5.534254

Code

#sd(fHH1$total)

ggplot(fHH1, aes(total)) + 
  geom_histogram(binwidth = .25, color = "black", 
                 fill = "white") + 
  xlab("Number in the house excluding head of household") +
  ylab("Count of households")

Code

cuts = cut(fHH1$age,
           breaks=c(15,20,25,30,35,40,45,50,55,60,65,70))
ageGrps <- data.frame(cuts,fHH1)

ggplot(data = ageGrps, aes(x = total)) +
  geom_histogram(binwidth = .25, color = "black", 
                 fill = "white") +
  facet_wrap(cuts) +
  xlab("Household size")

The first reveals a fair amount of variability in the number in each house; responses range from 0 to 16 with many of the respondents reporting between 1 and 5 people in the house. Like many Poisson distributions, this graph is right skewed. It clearly does not suggest that the number of people in a household is a normally distributed response.

The second further shows that responses can be reasonably modeled with a Poisson distribution when grouped by a key explanatory variable: age of the household head. These last two plots together suggest that Assumption 1 (Poisson Response) is satisfactory in this case study.

For Poisson random variables, the variance of \(Y\) (i.e., the square of the standard deviation of \(Y\)), is equal to its mean, where \(Y\) represents the size of an individual household. As the mean increases, the variance increases. So, if the response is a count and the mean and variance are approximately equal for each group of \(X\), a Poisson regression model may be a good choice. In Table @ref(tab:table1chp4) we display age groups by 5-year increments, to check to see if the empirical means and variances of the number in the house are approximately equal for each age group. This provides us one way in which to check the mean = variance assumption.

If there is a problem with this assumption, most often we see variances much larger than means. Here, as expected, we see more variability as age increases. However, it appears that the variance is smaller than the mean for lower ages, while the variance is greater than the mean for higher ages. Thus, there is some evidence of a violation of the mean=variance assumption (Assumption 3), although any violations are modest.

The Poisson regression model also implies that log(\(\lambda_i\)), not the mean household size \(\lambda_i\), is a linear function of age; i.e., \(log(\lambda_i)=\beta_0+\beta_1\textrm{age}_i\). Therefore, to check the linearity assumption (Assumption 4) for Poisson regression, we would like to plot log(\(\lambda_i\)) by age. Unfortunately, \(\lambda_i\) is unknown. Our best guess of \(\lambda_i\) is the observed mean number in the household for each age (level of \(X\)). Because these means are computed for observed data, they are referred to as empirical means. Taking the logs of the empirical means and plotting by age provides a way to assess the linearity assumption. The smoothed curve added to Figure @ref(fig:ageXnhouse) suggests that there is a curvilinear relationship between age and the log of the mean household size, implying that adding a quadratic term should be considered. This finding is consistent with the researchers’ hypothesis that there is an age at which a maximum household size occurs. It is worth noting that we are not modeling the log of the empirical means, rather it is the log of the true rate that is modeled. Looking at empirical means, however, does provide an idea of the form of the relationship between log(\(\lambda)\) and \(x_i\).

Fitting a GLM

modela = glm(total ~ age, family = poisson, data = fHH1)

Fitting a GLM

modela = glm(total ~ age, family = poisson, data = fHH1)
coef(summary(modela))

                Estimate   Std. Error   z value      Pr(>|z|)
(Intercept)  1.549942225 0.0502754106 30.829032 1.070156e-208
age         -0.004705881 0.0009363388 -5.025832  5.012548e-07

Question

How can we interpret these coefficients?

Deviance

modela_dev <- summary(modela)$deviance
modela_disp <- summary(modela)$dispersion

c(Deviance=modela_dev, Dispersion=modela_disp)

  Deviance Dispersion 
  2337.089      1.000

Confidence Intervals

confint(modela)

Confidence Intervals

confint(modela)

                   2.5 %       97.5 %
(Intercept)  1.451170100  1.648249185
age         -0.006543163 -0.002872717

exp(confint(modela))

                2.5 %    97.5 %
(Intercept) 4.2681057 5.1978713
age         0.9934782 0.9971314

Question

How should these be interpreted? What does significance look like?

These results suggest that by exponentiating the coefficient on age we obtain the multiplicative factor by which the mean count changes. In this case, the mean number in the house changes by a factor of \(e^{-0.0047}=0.995\) or decreases by 0.5% (since \(1-.995 = .005\)) with each additional year older the household head is; or, we predict a 0.47% increase in mean household size for a 1-year decrease in age of the household head (since \(1/.995=1.0047\)). The quantity on the left-hand side of Equation @ref(eq:rateRatio) is referred to as a rate ratio or relative risk, and it represents a percent change in the response for a unit change in X. In fact, for regression models in general, whenever a variable (response or explanatory) is logged, we make interpretations about multiplicative effects on that variable, while with unlogged variables we can reach our usual interpretations about additive effects.

Exponentiating the endpoints yields a confidence interval for the relative risk; i.e., the percent change in household size for each additional year older. Thus \((e^{-0.0065},e^{-0.0029})=(0.993,0.997)\) suggests that we are 95% confident that the mean number in the house decreases between 0.7% and 0.3% for each additional year older the head of household is. It is best to construct a confidence interval for the coefficient and then exponentiate the endpoints because the estimated coefficients more closely follow a normal distribution than the exponentiated coefficients. There are other approaches to constructing intervals in these circumstances, including profile likelihood, the delta method, and bootstrapping, and we will discuss some of those approaches later. In this case, for instance, the profile likelihood interval is nearly identical to the Wald-type (normal theory) confidence interval above.

If there is no association between age and household size, there is no change in household size for each additional year, so \(\lambda_X\) is equal to \(\lambda_{X+1}\) and the ratio \(\lambda_{X+1}/\lambda_X\) is 1. In other words, if there is no association between age and household size, then \(\beta_1=0\) and \(e^{\beta_1}=1\). Note that our interval for \(e^{\beta_1}\), (0.993,0.997), does not include 1, so the model with age is preferred to a model without age; i.e., age is significantly associated with household size. Note that we could have similarly confirmed that our interval for \(\beta_1\) does not include 0 to show the significance of age as a predictor of household size.

Your Turn

Exploration: Singapore Auto

Data: https://prof.mkjanssen.org/glm/data/SingaporeAuto.csv

Description: https://prof.mkjanssen.org/glm/data/SingaporeAutoDescriptions.pdf

Hands-on with GLMs I

Explorations

Exploring blocks

Case Study: Household Size in the Philippines

Setting the Stage

Importing Data

Questions

Exploratory Data Analysis

Fitting a GLM

Fitting a GLM

Deviance

Confidence Intervals

Confidence Intervals

Your Turn

Exploration: Singapore Auto

Exploring `blocks`