Exercises B

Statistics III - CdL SSE

Author
Affiliation

Tommaso Rigon

Università degli Studi di Milano-Bicocca

Homepage

The theoretical exercises below are quite difficult. In the exam, you can expect simplified versions; in their present form, most of them would be too challenging.

The data-analysis exercises, on the other hand, are broadly representative of what you may encounter in the final examination.

Most of these exercises are adapted, sometimes with minor modifications, from the textbooks Salvan et al. (2020) and Agresti (2015). You may consult these textbooks for additional exercises.

Data analysis

The Seed dataset in the MLGdata package comes from an experiment designed to assess whether, and to what extent, the amount of fertilizer affects seed germination. Twenty seeds were used, and the following variables were recorded for each seed:

  • fert indicates the amount of fertilizer,
  • x is a binary variable equal to 1 if the seed germinated, and 0 otherwise.

Import the data and then:

  1. Identify the statistical units, the response variable, and the covariate, specifying the type of each variable (continuous quantitative, discrete quantitative, nominal qualitative, or ordinal qualitative).

  2. Conduct an exploratory analysis to assess the relationship between the response and the covariate.

  3. Using pen and paper, specify an appropriate GLM with the canonical link function.

  4. Fit the GLM specified in (c) using R.

  5. Report the estimates and confidence intervals for the coefficients. Provide an interpretation of the obtained values.

  6. For each element in the output of summary() for the fitted glm object, indicate which quantity is being computed and match it to the corresponding formula in the slides.

The Wool dataset in the MLGdata package comes from a 3^3 factorial experiment on worsted yarn. The response y is the number of loading cycles until failure, and the three experimental factors are specimen length (x1), loading-cycle amplitude (x2), and load (x3). The data were analyzed in the classic paper by Box and Cox (1964):

Box, G. E. P. and Cox, D. R. (1964). An analysis of transformations (with discussion). Journal of the Royal Statistical Society, Series B, 26(2), 211–252.

For each of the three variables x1, x2, and x3, three levels were fixed:

  • Length: 250, 300, 350 mm (coded as -1, 0, 1)
  • Loading-cycle amplitude: 8, 9, 10 mm (coded as -1, 0, 1)
  • Load: 40, 45, 50 g (coded as -1, 0, 1)

Import the data and then:

  1. Identify the statistical units, the response variable, and the covariates, specifying the type of each variable (continuous quantitative, discrete quantitative, nominal qualitative, or ordinal qualitative).

  2. Using pen and paper, specify a normal linear model for the logarithmic transformation of the response.

  3. Fit in R the linear model from the previous point.

  4. Assess the goodness of fit of the model and compare the usual standard errors with heteroskedasticity-robust standard errors. Discuss whether a transformation of the response other than the logarithm might be more appropriate.

  5. Write down the expression of the estimated curve.

  6. For a test with length = 300 mm, loading-cycle amplitude = 10 mm, and load = 40 g, obtain a 95% confidence interval for the conditional mean of the log-response. Transform the interval back to the original scale and explain which feature of the conditional response distribution it estimates. For the same experimental conditions, obtain a 95% prediction interval for the response.

  7. For the same data, considering the untransformed response and using pen and paper, specify a generalized linear model with Gamma response and logarithmic link function.

  8. Fit in R the generalized linear model from the previous point.

  9. Write down the expression of the estimated curve.

  10. Report the estimates and confidence intervals for the coefficients. Provide an interpretation of the obtained values.

  11. For each element in the output of summary() for the fitted glm object, indicate which quantity is being computed and match it to the corresponding formula in the slides.

  12. Evaluate the goodness of fit of the Gamma model.

  13. Using the fitted Gamma model, obtain a 95% confidence interval for the mean response in an experiment with length = 300 mm, loading-cycle amplitude = 10 mm, and load = 40 g.

  14. Compare the results of the analysis based on the normal linear model with those of the analysis based on the Gamma model.

The dde dataset contains observations on a sample of n = 2312 pregnant women and can be downloaded here. The dataset was also analyzed in Statistica I. The original data are described in this study:

Longnecker, M. P., Klebanoff, M. A., Zhou, H., and Brock, J. W. (2001). Association between maternal serum concentration of the DDT metabolite DDE and preterm and small-for-gestational-age babies at birth. The Lancet, 358(9276), 110–114.

The variable DDE measures the concentration of DDE, a metabolite of DDT, in maternal serum during the third trimester of pregnancy (in \mu\mathrm{g}/\mathrm{L}). The variable GAD (gestational age at delivery) records the day of pregnancy on which delivery occurred.

DDT is extremely effective against malaria mosquitoes and is therefore widely used in areas where malaria is endemic. At the same time, DDT may pose a health risk, especially for pregnant women.

Research question. We are interested in the association between DDE and GAD, particularly in estimating the probability of preterm delivery (before 37 completed weeks, i.e., before day 259).

Import the data and then:

  1. Specify several reasonable models for the main research question. Select the most appropriate one using predictive accuracy and diagnostic tools.

  2. Is DDE concentration associated with gestational age at delivery?

  3. Predict GAD for a woman with a DDE concentration of 100\ \mu\mathrm{g}/\mathrm{L} and compare it with the prediction at 0\ \mu\mathrm{g}/\mathrm{L}. Note that the latter value lies just outside the observed range and therefore entails slight extrapolation.

  4. Compute the probabilities of preterm birth for all the observed values, namely \psi_i = \mathbb{P}(Y_i < 259 ; \bm{x}_i), \qquad i=1,\dots,n. Do they differ across models? Hint: use pnorm() for a normal model and pgamma() for a Gamma GLM. Pay attention to the parametrization.

  5. Estimate the probability of preterm delivery, \psi, at a DDE concentration of 100\ \mu\mathrm{g}/\mathrm{L} and compare it with the estimated probability at 0\ \mu\mathrm{g}/\mathrm{L}.

  6. ☠️ Optional. Obtain a prediction interval for GAD at a DDE concentration of 100\ \mu\mathrm{g}/\mathrm{L} and compare it with the interval at 0\ \mu\mathrm{g}/\mathrm{L}.

What affects the selling price of a house? The Houses dataset is available here and comprises n = 100 observations on recent home sales in Gainesville, Florida.

Variables listed are selling price (in thousands of dollars, price), size of house (in square feet, size), annual property tax bill (in dollars, taxes), number of bedrooms (beds), number of bathrooms (baths), and whether the house is new (new, with 1 = new and 0 = not new). The original source of the dataset is the textbook:

Agresti, A. (2015). Foundations of Linear and Generalized Linear Models. Wiley.

Because all 100 observations come from a single city, they do not support general conclusions about housing markets. For illustrative purposes, however, we treat them as a random sample from a conceptual population of home sales in this market and analyze how selling price is related to the recorded characteristics. The goal is to predict the sale price from the other variables.

Import the data and then:

  1. Summarize the data with descriptive statistics and plots.

  2. Using forward, backward, or stepwise selection with all five predictors, together with judgments about practical significance, select and interpret a (generalized) linear model for selling price.

  3. Check whether results depend on any influential observations.

This exercise is intentionally left without strict instructions, so you can freely explore the data and develop your own modeling strategy. As is often the case, multiple approaches can yield reasonable results.

Theoretical

Let Y be a random variable with an inverse Gaussian distribution, with support (0, \infty) and probability density function

p(y ; \xi, \lambda) = \left( \frac{\lambda}{2 \pi} \right)^{1/2}y^{-3/2}e^{\sqrt{\lambda\xi}}\exp\left\{-\frac{1}{2}\left(\frac{\lambda}{y} + \xi y\right)\right\}, \qquad y > 0, \; \xi > 0, \; \lambda > 0. Show that this distribution belongs to the exponential-dispersion family and identify its characteristic elements: the canonical parameter, the functions a_i(\cdot), b(\cdot), and c(\cdot), the dispersion parameter, and the variance function.

Let Y have a negative binomial distribution, representing the number of independent Bernoulli trials with constant success probability \pi \in (0,1) required to obtain k successes. Its support is S = \{k, k+1, \dots\}, and its probability mass function is p(y; \pi) = \binom{y-1}{k-1} \pi^k (1-\pi)^{\,y-k}, \qquad y \in S.

  1. Assuming that k is known, verify that this distribution belongs to the exponential-dispersion family and identify its characteristic elements: the canonical parameter, the functions a_i(\cdot), b(\cdot), and c(\cdot), the dispersion parameter, and the variance function.

  2. Using the properties of exponential families, recover the well-known relations
    \mathbb{E}(Y) = \frac{k}{\pi}, \qquad \operatorname{Var}(Y) = \frac{k(1-\pi)}{\pi^2}.

  3. Show that the variance function is quadratic in \mu.

  4. Does the distribution still constitute an exponential-dispersion family if k is treated as unknown?

Consider the general likelihood equations for a GLM:

  1. Specialize the general formulas and re-obtain the likelihood equations for the binomial regression model (Beetles data with p = 2) presented in the slides, using the canonical link function.

  2. Specialize the general formulas and re-obtain the likelihood equations for the Poisson regression model (Aids data with p = 2) presented in the slides, using the canonical link function.

  3. Obtain the likelihood equations for a binomial regression model (Beetles data with p = 2) using the Cauchy link, defined as g(\mu) = \tan\{\pi(\mu - 1/2)\}. This link maps (0,1) onto \mathbb{R}.

Hint: Parts i. and ii. require almost no calculations, only substitutions. They are intended to familiarize you with the notation and the general formulas. Part iii. is more involved; use the derivative of the Cauchy link, g'(\mu) = \pi/\sin^2(\pi\mu).

Explicitly derive the contribution d_i of a single observation to the deviance of a Gamma generalized linear model. Then write down the deviance D(\hat{\bm{\mu}}; \boldsymbol{y}) for a sample of size n.

Let Y_1,\dots,Y_n be independent random variables, with corresponding values x_1,\dots,x_n of a univariate quantitative covariate. Consider a GLM with no intercept and p=1, i.e., \eta_i = \beta x_i, \qquad i=1,\dots,n, under each of the following exponential-dispersion models, specified through their mean-variance pairs:

  1. \mathbb{E}(Y_i)=\mu_i and \operatorname{Var}(Y_i)=\phi, where \mu_i \in \mathbb{R} and the constant \phi>0 is known.
  2. \mathbb{E}(Y_i)=\operatorname{Var}(Y_i)=\mu_i, where \mu_i \in (0,\infty).
  3. \mathbb{E}(Y_i)=\mu_i and \operatorname{Var}(Y_i)=\mu_i(1-\mu_i), where \mu_i \in (0,1).

For each model, use the canonical link and:

  1. Specify the statistical model in terms of the “standard” probability distribution and write the log-likelihood as a function of \beta.

  2. Derive the score function \ell_*(\beta) and both the observed and expected Fisher information.

  3. Derive the equation that characterizes the maximum likelihood estimator \hat\beta, obtain a closed-form expression when possible, and give an approximation to its sampling distribution.

  4. For a fixed covariate value x_i, construct an approximate 1-\alpha confidence interval for the linear predictor \eta_i = \beta x_i. Then derive the corresponding confidence interval for \mu_i. Indicate in which of the three cases these intervals are also exact 1-\alpha intervals.

Suppose that Y_1,\dots,Y_{n_A+n_B} are independent Poisson random variables whose means satisfy g(\mu_i) = \beta_1 + \beta_2 x_i, where x_i = 1 for i = 1, \dots, n_A (group A) and x_i = 0 for i = n_A + 1, \dots, n_A + n_B (group B); thus, x_i is a dummy variable.

  1. Show that, for the log link g(\mu)=\log(\mu), the GLM likelihood equations imply that the fitted group means \hat\mu_A and \hat\mu_B equal the corresponding sample means.

  2. Using the likelihood equations, show that the same result holds for any differentiable link function for this Poisson model, provided that an interior maximum likelihood estimate exists.

  3. Using the likelihood equations, show that the same result holds for any GLM of the form g(\mu_i) = \beta_1 + \beta_2 x_i with a binary indicator predictor, under the same regularity condition. Note: If the weights \omega_i are not constant, then \hat\mu_A and \hat\mu_B equal the corresponding weighted sample means.

In a generalized linear model that uses a non-canonical link function, explain why it need not be true that
\sum_{i=1}^n \hat{\mu}_i = \sum_{i=1}^n y_i. Consequently, the response residuals need not have mean zero.

Then explain why, in a GLM with a canonical link and constant weights \omega_i, an intercept must be included for this equality to be guaranteed.

When selecting explanatory variables for a model, why is it inadequate to choose the model with the smallest deviance (or, for a linear model, the largest R^2)?

References

Agresti, A. (2015), Foundations of Linear and Generalized Linear Models, Wiley.
Salvan, A., Sartori, N., and Pace, L. (2020), Modelli lineari generalizzati, Springer.