How Harmful is Smoking?
MATH 4939
January 2026
Last update: January 11 2026 20:38
1 Setup
This is an example of an This R script is designed to be run either interactively one line at a time, or as a Rmarkdown script that creates an html file with interactive 3D widgets.
In a clean version of a report, this portion of the Rmarkdown script would be hidden by using a ‘include=FALSE’ option for ‘chunk’. That is the header line for the chunk would read:
#+ setup, include=FALSE
knitr::opts_chunk$set(comment = ' ') # For interactive use, the following sets the working directory
# to this file's location.
#
setwd(this.path::here())
#
# It's an alternative to using menus:
# Session|Set Working Directory|To Source File Location
#
cex <- if(interactive()) 1.5 else 1 # setting a character expansion factor
widget <- function(){
#
# Does nothing when run interactively,
# displays interactive 3d widgets when the script
# is compiled into an HTML document.
#
if(interactive()) invisible(NULL)
else {
spin(fov=10)
Id3d(c("Greece", "Japan", "United States", "Canada"))
rglwidget()
}
}
summ <- function(fit) {
#
# Display a short summary of regression coefficient
# estimates for a linear model.
#
waldf(fit)[,c('coef','se','p-value','t-value','DF')]
}
library(p3d) # remotes::install_github('gmonette/p3d') Loading required package: rgllibrary(spida2) # remotes::install_github('gmonette/spida2')
Attaching package: 'spida2' The following objects are masked from 'package:p3d':
cell, center, ConjComp, dell, disp, ell, ell.conj, ellbox, ellplus,
ellpt, ellptc, ellpts, ellptsc, elltan, elltanc, na.include, uv The following object is masked from 'package:base':
grepv#
# For a clean run you can suppress messages when loading a package:
suppressPackageStartupMessages(library(car)) # for qqPlot
data(Smoking4) # in the spida2 package
?Smoking4Snapshot of univariate distributions in a dataset using uniform quantile plots: i.e. data values on the y-axis and uniform quantiles on the x-axis.
This shows you the data standing shoulder-to-shoulder from the shortest to the tallest. It’s easy to spot outliers, skewness, kurtosis, gaps, discrete versus continuous, some invalid data values, consistent frequencies for categorical data, etc.
xqplot(Smoking4) # uniform quantile plots
xqplot(Smoking4, ptype = 'n') # normal quantile plots
1.1 Question
Describe some aspects of these variables that are readily seen with uniform quantile plots.
Select rows of the data frame that contain values with for both sexes combined:
dd <- subset(Smoking4, sex == 'BTSX')
dim(dd) [1] 189 34rownames(dd) <- dd$Country # for labelling in 3d plots
names(dd) [1] "Country" "Continent"
[3] "LE" "CigCon"
[5] "LE.q" "Cont"
[7] "Cont2" "HealthExpPC"
[9] "Year" "HE"
[11] "country" "iso3"
[13] "region" "HealthExpPC.Govt.exch"
[15] "HealthExpPC.Tot.ppp" "HealthExpPC.Govt.ppp"
[17] "HealthExpPC.Tot.exch" "total"
[19] "govt" "private"
[21] "sex" "lifeexp.Birth"
[23] "lifeexp.At60" "smoking.tobacco.current"
[25] "smoking.tobacco.daily" "smoking.cig.current"
[27] "smoking.cig.daily" "Pop.Total"
[29] "Pop.MedAge" "Pop.pCntUnder15"
[31] "Pop.pCntOver60" "Pop.pCntAnnGrowth"
[33] "consumption.cigPC" "hiv_prev15_49"dd <- within(
dd,
{
# more information variable names
# although more awkward to type
`Cigarettes/day` <- CigCon/365.25
`Life Expectancy` <- LE
region <- factor(region)
}
)
# Initialize 3d graphics window
Init3d(cex = if(interactive()) 1.5 else 1)
Plot3d(`Life Expectancy` ~
`Cigarettes/day` + HealthExpPC | region, dd,
theta = 0, phi = 0, fov = 0) region col
1 AFR blue
2 AMR #008800
3 EMR orange
4 EUR magenta
5 SEAR darkcyan
6 WPR red Use left mouse to rotate, middle mouse (or scroll) to zoom, right mouse to change perspectiveAxes3d()
Id3d(pad = 1, cex = 1.5) NULLId3d('China') [1] "China"widget()
Least-squares line to estimate relationship between Life Expectancy and Cigarette consumption
fitlin <- lm(`Life Expectancy` ~ `Cigarettes/day`, dd)
summary(fitlin)
Call:
lm(formula = `Life Expectancy` ~ `Cigarettes/day`, data = dd)
Residuals:
Min 1Q Median 3Q Max
-19.4763 -5.7922 0.7627 5.3289 17.8041
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 47.9836 1.3668 35.107 < 2e-16 ***
`Cigarettes/day` 3.1297 0.3269 9.573 7.18e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8.131 on 102 degrees of freedom
(85 observations deleted due to missingness)
Multiple R-squared: 0.4732, Adjusted R-squared: 0.4681
F-statistic: 91.64 on 1 and 102 DF, p-value: 7.184e-16qqPlot(fitlin)
Angola Norway
5 127Short summary
summ(fitlin) coef se p-value t-value DF
(Intercept) 47.983620 1.3667689 9.129311e-59 35.107339 102
`Cigarettes/day` 3.129681 0.3269322 7.184013e-16 9.572874 102Fit3d(fitlin, lwd = 3, resid = TRUE)
widget()
1.2 Questions
- What determines the position of the least-squares regression line?
- What seems wrong with this model?
Let’s try to include curvature with a quadratic model:
fitq <- lm(`Life Expectancy` ~
`Cigarettes/day` + I(`Cigarettes/day`^2), dd)
summ(fitq) coef se p-value t-value DF
(Intercept) 43.5421197 1.62945686 1.349121e-47 26.721861 101
`Cigarettes/day` 6.6164243 0.86341445 1.140729e-11 7.663092 101
I(`Cigarettes/day`^2) -0.4232845 0.09819861 3.792242e-05 -4.310494 101qqPlot(fitq)
Angola Panama
5 131Fit3d(fitq, col = 'red', lwd = 3)
widget()
Comparing the two models
anova(fitlin, fitq) Analysis of Variance Table
Model 1: `Life Expectancy` ~ `Cigarettes/day`
Model 2: `Life Expectancy` ~ `Cigarettes/day` + I(`Cigarettes/day`^2)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 102 6743.6
2 101 5695.8 1 1047.8 18.58 3.792e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summ(fitq) coef se p-value t-value DF
(Intercept) 43.5421197 1.62945686 1.349121e-47 26.721861 101
`Cigarettes/day` 6.6164243 0.86341445 1.140729e-11 7.663092 101
I(`Cigarettes/day`^2) -0.4232845 0.09819861 3.792242e-05 -4.310494 101summ(fitlin) coef se p-value t-value DF
(Intercept) 47.983620 1.3667689 9.129311e-59 35.107339 102
`Cigarettes/day` 3.129681 0.3269322 7.184013e-16 9.572874 1021.3 Questions
- Is ‘fitq’ a linear model?
- Now that we have a quadratic model, it looks like the the ‘effect of smoking’ is much larger? Can you explain why?
- What is the purpose of the “I” in “I(
Cigarettes/day^2)”? - What is the estimated ‘effect of smoking’ (in Cigarettes/day) according to the regression output for ‘fitq’?
summ(fitq) coef se p-value t-value DF
(Intercept) 43.5421197 1.62945686 1.349121e-47 26.721861 101
`Cigarettes/day` 6.6164243 0.86341445 1.140729e-11 7.663092 101
I(`Cigarettes/day`^2) -0.4232845 0.09819861 3.792242e-05 -4.310494 1012 “Controlling” for another variable
Since we are going to replot this data over and over again, I’m going to define a function to work like a macro to make that easier.
This is very bad programming style for serious work, such as when writing code for a package.
f <- function() {
Plot3d(`Life Expectancy` ~
`Cigarettes/day` + HealthExpPC | region, dd,
theta = 0, phi = 0, fov = 0, ylim = c(30,75))
Axes3d()
}
f() region col
1 AFR blue
2 AMR #008800
3 EMR orange
4 EUR magenta
5 SEAR darkcyan
6 WPR red Use left mouse to rotate, middle mouse (or scroll) to zoom, right mouse to change perspectiveExplore how Health Expenditures per capita comes into the picture.
What do you see?
Id3d(pad = 1, cex = 2) NULLwidget()
2.1 Question
The distribution of HealthExpPC is highly ________?
fitlin2 <- lm(`Life Expectancy` ~
`Cigarettes/day` + HealthExpPC, dd)
Fit3d(fitlin2, col = 'green')
Fit3d(fitlin)
widget()
summ(fitlin2) coef se p-value t-value DF
(Intercept) 48.063444647 1.2152254642 2.786550e-63 39.551051 101
`Cigarettes/day` 2.384473425 0.3229312088 4.484208e-11 7.383843 101
HealthExpPC 0.003139478 0.0005928222 6.946486e-07 5.295817 101summ(fitlin) coef se p-value t-value DF
(Intercept) 47.983620 1.3667689 9.129311e-59 35.107339 102
`Cigarettes/day` 3.129681 0.3269322 7.184013e-16 9.572874 1022.2 Questions
- Using the model ‘fitlin2’, which variable (cigarette consumption or health expenditures per capita) provides stronger evidence of its positive effect on Life Expectancy?
- Does controlling for health expenditures with the simple additive model ‘fitlin2’ lead to an estimate of the effect of smoking that seems closer to its causal effect?
- How could we improve this model?
3 Accounting for non-linearity in the relationship with HealthExpPC
fitlog2 <- lm(`Life Expectancy` ~
(
`Cigarettes/day` +
I(`Cigarettes/day`^2)
) + # additive
log(HealthExpPC), dd)
fitlog2int <- lm(`Life Expectancy` ~
(
`Cigarettes/day` +
I(`Cigarettes/day`^2)
) * # interactive
log(HealthExpPC), dd)
summ(fitlin) coef se p-value t-value DF
(Intercept) 47.983620 1.3667689 9.129311e-59 35.107339 102
`Cigarettes/day` 3.129681 0.3269322 7.184013e-16 9.572874 102summ(fitq) coef se p-value t-value DF
(Intercept) 43.5421197 1.62945686 1.349121e-47 26.721861 101
`Cigarettes/day` 6.6164243 0.86341445 1.140729e-11 7.663092 101
I(`Cigarettes/day`^2) -0.4232845 0.09819861 3.792242e-05 -4.310494 101summ(fitlog2) coef se p-value t-value DF
(Intercept) 32.192205 1.4601490 3.489711e-40 22.047205 100
`Cigarettes/day` 2.826029 0.6583501 4.090443e-05 4.292593 100
I(`Cigarettes/day`^2) -0.225268 0.0670480 1.105183e-03 -3.359801 100
log(HealthExpPC) 4.081924 0.3553103 5.590550e-20 11.488335 100summ(fitlog2int) coef se p-value
(Intercept) 21.9386136 2.64524971 5.972133e-13
`Cigarettes/day` 9.8782485 2.18819780 1.770068e-05
I(`Cigarettes/day`^2) -1.0281430 0.36330499 5.647787e-03
log(HealthExpPC) 6.5554857 0.66220778 2.000524e-16
`Cigarettes/day`:log(HealthExpPC) -1.3921132 0.33955197 8.536686e-05
I(`Cigarettes/day`^2):log(HealthExpPC) 0.1443672 0.05077192 5.432081e-03
t-value DF
(Intercept) 8.293589 98
`Cigarettes/day` 4.514331 98
I(`Cigarettes/day`^2) -2.829972 98
log(HealthExpPC) 9.899439 98
`Cigarettes/day`:log(HealthExpPC) -4.099853 98
I(`Cigarettes/day`^2):log(HealthExpPC) 2.843445 98Since these models are nested we can compare them using sequential likelihood ratio tests (LRTs or Wilk’s tests)
anova(fitlin, fitq, fitlog2, fitlog2int) Analysis of Variance Table
Model 1: `Life Expectancy` ~ `Cigarettes/day`
Model 2: `Life Expectancy` ~ `Cigarettes/day` + I(`Cigarettes/day`^2)
Model 3: `Life Expectancy` ~ (`Cigarettes/day` + I(`Cigarettes/day`^2)) +
log(HealthExpPC)
Model 4: `Life Expectancy` ~ (`Cigarettes/day` + I(`Cigarettes/day`^2)) *
log(HealthExpPC)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 102 6743.6
2 101 5695.8 1 1047.8 50.461 1.961e-10 ***
3 100 2455.3 1 3240.5 156.058 < 2.2e-16 ***
4 98 2035.0 2 420.3 10.121 0.000101 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Each model seems to improve on the previous one as evaluated by statistical significance.
f() region col
1 AFR blue
2 AMR #008800
3 EMR orange
4 EUR magenta
5 SEAR darkcyan
6 WPR red Use left mouse to rotate, middle mouse (or scroll) to zoom, right mouse to change perspectiveFit3d(fitlog2) Warning in log(HealthExpPC): NaNs producedFit3d(fitlog2int, col = 'red') Warning in log(HealthExpPC): NaNs producedwidget()
f() region col
1 AFR blue
2 AMR #008800
3 EMR orange
4 EUR magenta
5 SEAR darkcyan
6 WPR red Use left mouse to rotate, middle mouse (or scroll) to zoom, right mouse to change perspectiveplotmath3d(6,25,0,
text = expression(
E(LE) == beta[0] + beta[1] *C +
beta[2]*C^2 + beta[3]*log(H) + beta[4]*C*' '*log(H) +
beta[5]*C^2*log(H)
),
cex = if(interactive()) 2 else 0.8)
Fit3d(fitlog2int, col = 'red') Warning in log(HealthExpPC): NaNs producedwidget()
3.1 Questions
- Fit a model that controls for Health Expenditure per capita using a different transformation than a log. Does it seem to work better? Why or why not?
- What assumption about the relationship between Life Expectancy and Health Expenditure per capita is implied by the use of log(HEPC) versus HEPC itself?
4 Estimating the ‘effect’ of C
The statistical ‘effect’ of Cigarettes/day in the interaction model is the partial derivative of \(E(LE)\) with respect to ‘Cigarettes/day’
\[E(LE) = \beta_0 + \beta_1 C + \beta_2 C^2 + \beta_3 \log(H) + \beta_4 C \log(H) + \beta_5 C^2 \log(H)\]
We are fitting a linear model which means, to a statistician, that it is a **linear function* of the \(\beta\) parameters.
What is the effect of C in this model?
It is the partial derivative of \(E(LE)\) with respect to \(C\).
\[\frac{\partial}{\partial C}E(LE) = \beta_1 + \beta_2 2 \times C + \beta_4 \log(H) + \beta_5 2 \times C \times \log(H)\]
With a simple additive model (a model that is linear in the parameters and in the predictors), the effect of C is simply the coefficient of \(C\).
In any other model it is a function that depends on some of the values of the predictors.
How do you estimate the effect of C in a linear model?
If the model is linear in the \(\beta\)s then any derivative with respect to predictors must also be linear in the \(\beta\)s.
This means that a partial derivative \(\eta = \frac{\partial}{\partial C}E(LE)\) at particular values of the predictors can be expressed in the form
\[\eta = L \beta\]
and we get the estimated partial derivative with:
\[\hat{\eta} = L \hat{\beta}\]
where, in our example, \(L\) is a row matrix so that multiplying it by the parameter vector gives us the desired linear combination of the parameters:
\[\eta = L \beta = \begin{pmatrix} L_0 & L_1 & L_2 & L_3 & L_4 & L_5 \end{pmatrix} \begin{pmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ \beta_3 \\ \beta_4 \\ \beta_5 \end{pmatrix}\]
For the example above, we substitute the desired values for \(C\) and \(H\) in the row vector:
\[L = \begin{pmatrix}0 & 1 & 2C & 0 & \log(H) & 2C \log(H) \end{pmatrix}\]
With \(C = 6\) and \(H = 1000\), we create the \(L\) matrix:
L <- matrix(c(0, 1, 2*6, 0, log(1000), 2 * 6 * log(1000)), nrow = 1)
L [,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0 1 12 0 6.907755 82.89306then use the wald function in the spida2 package:
wald(fitlog2int, L) numDF denDF F-value p-value
1 1 98 0.1183191 0.7316
Estimate Std.Error DF t-value p-value Lower 0.95 Upper 0.95
[1,] -0.108806 0.316319 98 -0.343975 0.7316 -0.73653 0.518919or the lht function in the car package:
lht(fitlog2int, L)
Linear hypothesis test:
Model 1: restricted model
Model 2: `Life Expectancy` ~ (`Cigarettes/day` + I(`Cigarettes/day`^2)) *
log(HealthExpPC)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 99 2037.4
2 98 2035.0 1 2.4569 0.1183 0.7316to carry out the estimation of \(\eta\) as well as calculating the estimated SE of \(\hat{\eta}\), etc.
To simultaneously estimate effects at two points:
L <- rbind(
"C = 6, H = 1,000" = c(0, 1, 2*6, 0, log(1000), 2 * 6 * log(1000)),
"C = 8, H = 2,000" = c(0, 1, 2*8, 0, log(2000), 2 * 8 * log(2000))
)
L [,1] [,2] [,3] [,4] [,5] [,6]
C = 6, H = 1,000 0 1 12 0 6.907755 82.89306
C = 8, H = 2,000 0 1 16 0 7.600902 121.61444wald(fitlog2int, L) numDF denDF F-value p-value
1 2 98 0.8358975 0.43655
Estimate Std.Error DF t-value p-value Lower 0.95 Upper 0.95
C = 6, H = 1,000 -0.108806 0.316319 98 -0.343975 0.73160 -0.736530 0.518919
C = 8, H = 2,000 0.403779 0.586405 98 0.688567 0.49272 -0.759922 1.567481Note that first row of the output provided information for the test of the joint hypothesis: \[H_0: \eta_1 = \eta_2 = 0 \mathrm{\ versus\ } H_A: \mathrm{\ at\ least\ one\ of\ } \eta_1 \mathrm{\ or\ } \eta_2 \neq 0.\]
To do the same for lots of points:
ee <- expression(list(0, 1, 2*C, 0, log(H), 2 * C * log(H)))
ee expression(list(0, 1, 2 * C, 0, log(H), 2 * C * log(H)))Create a data frame with desired combinations of values of C and H
df <- expand.grid(C = seq(0, 10, by = 2),
H = c(2000, 3000))
df C H
1 0 2000
2 2 2000
3 4 2000
4 6 2000
5 8 2000
6 10 2000
7 0 3000
8 2 3000
9 4 3000
10 6 3000
11 8 3000
12 10 3000Llist <- eval(ee, df)
Llist [[1]]
[1] 0
[[2]]
[1] 1
[[3]]
[1] 0 4 8 12 16 20 0 4 8 12 16 20
[[4]]
[1] 0
[[5]]
[1] 7.600902 7.600902 7.600902 7.600902 7.600902 7.600902 8.006368 8.006368
[9] 8.006368 8.006368 8.006368 8.006368
[[6]]
[1] 0.00000 30.40361 60.80722 91.21083 121.61444 152.01805 0.00000
[8] 32.02547 64.05094 96.07641 128.10188 160.12735L <- do.call(cbind, Llist)
L [,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0 1 0 0 7.600902 0.00000
[2,] 0 1 4 0 7.600902 30.40361
[3,] 0 1 8 0 7.600902 60.80722
[4,] 0 1 12 0 7.600902 91.21083
[5,] 0 1 16 0 7.600902 121.61444
[6,] 0 1 20 0 7.600902 152.01805
[7,] 0 1 0 0 8.006368 0.00000
[8,] 0 1 4 0 8.006368 32.02547
[9,] 0 1 8 0 8.006368 64.05094
[10,] 0 1 12 0 8.006368 96.07641
[11,] 0 1 16 0 8.006368 128.10188
[12,] 0 1 20 0 8.006368 160.12735wald(fitlog2int, L) numDF denDF F-value p-value
1 4 98 11.05301 <.00001
Estimate Std.Error DF t-value p-value Lower 0.95 Upper 0.95
[1,] -0.703068 1.047190 98 -0.671385 0.50355 -2.781183 1.375047
[2,] -0.426356 0.717070 98 -0.594581 0.55349 -1.849357 0.996645
[3,] -0.149644 0.441127 98 -0.339232 0.73516 -1.025047 0.725758
[4,] 0.127067 0.371242 98 0.342276 0.73288 -0.609650 0.863785
[5,] 0.403779 0.586405 98 0.688567 0.49272 -0.759922 1.567481
[6,] 0.680491 0.901523 98 0.754824 0.45216 -1.108552 2.469534
[7,] -1.267521 1.128701 98 -1.122992 0.26418 -3.507391 0.972348
[8,] -0.756666 0.774947 98 -0.976410 0.33127 -2.294524 0.781192
[9,] -0.245811 0.496426 98 -0.495161 0.62160 -1.230953 0.739331
[10,] 0.265044 0.460057 98 0.576112 0.56586 -0.647924 1.178013
[11,] 0.775900 0.704403 98 1.101499 0.27338 -0.621965 2.173764
[12,] 1.286755 1.048824 98 1.226854 0.22282 -0.794603 3.368112as a data frame:
waldf(fitlog2int, L) coef se U2 L2 p-value t-value DF L.1
1 -0.7030680 1.0471900 1.3913120 -2.7974480 0.5035546 -0.6713853 98 0.000000
2 -0.4263562 0.7170695 1.0077828 -1.8604953 0.5534936 -0.5945814 98 0.000000
3 -0.1496445 0.4411273 0.7326101 -1.0318990 0.7351609 -0.3392319 98 0.000000
4 0.1270673 0.3712420 0.8695512 -0.6154166 0.7328760 0.3422763 98 0.000000
5 0.4037791 0.5864049 1.5765889 -0.7690307 0.4927226 0.6885671 98 0.000000
6 0.6804909 0.9015229 2.4835367 -1.1225550 0.4521649 0.7548237 98 0.000000
7 -1.2675213 1.1287006 0.9898799 -3.5249225 0.2641847 -1.1229916 98 0.000000
8 -0.7566661 0.7749473 0.7932285 -2.3065607 0.3312656 -0.9764098 98 0.000000
9 -0.2458109 0.4964264 0.7470419 -1.2386638 0.6215950 -0.4951608 98 0.000000
10 0.2650443 0.4600572 1.1851587 -0.6550701 0.5658605 0.5761116 98 0.000000
11 0.7758995 0.7044031 2.1847056 -0.6329066 0.2733774 1.1014994 98 0.000000
12 1.2867547 1.0488242 3.3844032 -0.8108937 0.2228166 1.2268545 98 0.000000
L.2 L.3 L.4 L.5 L.6
1 1.000000 0.000000 0.000000 7.600902 0.000000
2 1.000000 4.000000 0.000000 7.600902 30.403610
3 1.000000 8.000000 0.000000 7.600902 60.807220
4 1.000000 12.000000 0.000000 7.600902 91.210830
5 1.000000 16.000000 0.000000 7.600902 121.614439
6 1.000000 20.000000 0.000000 7.600902 152.018049
7 1.000000 0.000000 0.000000 8.006368 0.000000
8 1.000000 4.000000 0.000000 8.006368 32.025470
9 1.000000 8.000000 0.000000 8.006368 64.050941
10 1.000000 12.000000 0.000000 8.006368 96.076411
11 1.000000 16.000000 0.000000 8.006368 128.101881
12 1.000000 20.000000 0.000000 8.006368 160.127351combine with the original:
df2 <- cbind(waldf(fitlog2int, L), df)
df2[,c('coef','se','C','H')] coef se C H
1 -0.7030680 1.0471900 0 2000
2 -0.4263562 0.7170695 2 2000
3 -0.1496445 0.4411273 4 2000
4 0.1270673 0.3712420 6 2000
5 0.4037791 0.5864049 8 2000
6 0.6804909 0.9015229 10 2000
7 -1.2675213 1.1287006 0 3000
8 -0.7566661 0.7749473 2 3000
9 -0.2458109 0.4964264 4 3000
10 0.2650443 0.4600572 6 3000
11 0.7758995 0.7044031 8 3000
12 1.2867547 1.0488242 10 3000If you only want a few partial derivatives, e.g. “what is the estimated effect of Cig/day when Cig/day = 4 and HEPC = 2000?” the quick way to create the \(L\) matrix is to differentiate each term in the output of the regression coefficients.
For example, with
summ(fitlog2int) coef se p-value
(Intercept) 21.9386136 2.64524971 5.972133e-13
`Cigarettes/day` 9.8782485 2.18819780 1.770068e-05
I(`Cigarettes/day`^2) -1.0281430 0.36330499 5.647787e-03
log(HealthExpPC) 6.5554857 0.66220778 2.000524e-16
`Cigarettes/day`:log(HealthExpPC) -1.3921132 0.33955197 8.536686e-05
I(`Cigarettes/day`^2):log(HealthExpPC) 0.1443672 0.05077192 5.432081e-03
t-value DF
(Intercept) 8.293589 98
`Cigarettes/day` 4.514331 98
I(`Cigarettes/day`^2) -2.829972 98
log(HealthExpPC) 9.899439 98
`Cigarettes/day`:log(HealthExpPC) -4.099853 98
I(`Cigarettes/day`^2):log(HealthExpPC) 2.843445 98to get the effect of \(H\) when
C = 2 and H = 1000
write the coefficients of the partial derivative with respect to H for each term into a hypothesis matrix:
L <- rbind(
c(
0, # d/d(H) Intercept
0, # d/d(H) `Cigarettes/day`
0, # d/d(H) I(`Cigarettes/day`^2)
1/1000, # d/d(H) log(HealthExpPC)
2/1000, # d/d(H) `Cigarettes/day`:log(HealthExpPC)
2^2/1000 # d/d(H) I(`Cigarettes/day`^2):log(HealthExpPC)
)
)
L [,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0 0 0 0.001 0.002 0.004wald(fitlog2int, L) numDF denDF F-value p-value
1 1 98 96.13471 <.00001
Estimate Std.Error DF t-value p-value Lower 0.95 Upper 0.95
[1,] 0.004349 0.000444 98 9.804831 <.00001 0.003469 0.005229Equivalently, transform HealthExpPC
dd $ logHEPC <- log(dd$HealthExpPC)
Plot3d(`Life Expectancy` ~
`Cigarettes/day` + logHEPC | region, dd,
theta = 0, phi = 0, fov = 0, ylim = c(30,75)) region col
1 AFR blue
2 AMR #008800
3 EMR orange
4 EUR magenta
5 SEAR darkcyan
6 WPR red Use left mouse to rotate, middle mouse (or scroll) to zoom, right mouse to change perspectiveAxes3d()
Id3d(pad = 1, cex = 1.5) NULLdd$Region <- as.factor(dd$region)
fit2 <- lm(`Life Expectancy` ~
`Cigarettes/day` + logHEPC,
subset(dd, region != 'AFR'))
Fit3d(fit2)
widget()
Models adjusting for Region
Plot3d(`Life Expectancy` ~
`Cigarettes/day` + logHEPC | region, dd,
theta = 0, phi = 0, fov = 0, ylim = c(30,75)) region col
1 AFR blue
2 AMR #008800
3 EMR orange
4 EUR magenta
5 SEAR darkcyan
6 WPR red Use left mouse to rotate, middle mouse (or scroll) to zoom, right mouse to change perspectivefitlogRegion <- lm(`Life Expectancy` ~
(`Cigarettes/day` + I(`Cigarettes/day`^2)) +
logHEPC + region, dd)
Fit3d(fitlogRegion)
widget()
4.1 Questions
- What lessons can we learn from this example?
- Explore some questions you find interesting using this data set. For example, the variables ‘govt’ and ‘private’ show the amount spent on Heath Expenses per capita through the government and the private sector. Do these two routes have a different impact? Or, compare Life Expectancy for women and for men. Etc.
5 Hints
1(b): It depends! Linear in what? The ‘betas’ (parameters) or the ‘xs’ (predictors)? When statisticians talk about ‘linear’ models, they almost always mean “linear in the parameters”.
1(d): Always think of ‘effects’ as PARTIAL DERIVATIVES
- Take the fitted model and evaluate the partial
derivative with respect to the variable of interest.
- In a simple additive model this will be a constant
for all values of the predictors.
- If the model is more interesting, then the
partial derivative with respect to a variable
may depend on the value of that
variable and on the values of other variables.
- WORK OUT THE FORMULA FOR THE PARTIAL DERIVATIVE.
- If a model is linear in the betas (parameters)
then partial derivatives with respect to xs (predictors)
will also be linear in the betas. Why?
- The fact that they can be expressed as linear functions
of betas means that they can be estimated with
Wald tests by specifying a hypothesis matrix to
be applied to the vector of estimated betas.
6 Summary
The fact that two variables X and Y are associated in a data set does not mean that actively changing the level of X (intervening with X) would result in a change in the level of Y similar to the difference in Y observed in the data set.
To estimate the causal effect of X on Y, we can try to control for “confounding factors”, but our ability to estimate a causal effect will depend on a number of factors:
- Have we identified and measured accurately a sufficient set of confounding factors? (more on this soon).
- Have we identified the correct form of the relationship of between the outcome and confounding factors?
- Have we avoided controlling for ‘mediating’ factors that are responsible for the causal effect of X on Y.
- By correlation we mean the STATISTICAL ASSOCIATION between variables in a data set. Although the word ‘effect’ connotes causation in everyday speech, statisticians use it in contexts that do not imply causation. This is an unfortunate source of miscommunication, a problem to be aware of when interpreting statistical results.
- By causation (‘X causes Y’) we mean how much a variable Y is expected to differ when the variable X is changed through intervention (‘do X’), in contrast with passive observation of a difference in X (‘see X’).
- An important approach to causality, developed over the past four decades, uses the notion of ‘counterfactuals’. The causal effect of X on Y, supposing that X has two levels: 0 and 1, say, is to consider the difference in the expected value of Y if the same subject were to receive either level 0 or level 1 of X. For example, level 1 could be receiving a COVID vaccine and level 0 the injection of a placebo instead. If Y is 1 if the subject gets COVID and 0 if not, then the regression coefficient of Y on X measures the difference in the probability of COVID ‘caused’ by the vaccine. (Hernán and Robins (2020))
- The problem is that a subject can’t get both the vaccine and the placebo, so we can’t directly measure the difference in Y for a particular subject. Thus one of the two conditions is a counterfactual for each subject.
- One approach to the study of causality uses directed acyclic graphs or DAGs, often referred to as causal graphs when used for causal inference (Pearl and Mackenzie (2018)).
Correlation is not necessarily causation but sometimes it is. The challenge is being able to tell when.
6.1 Questions
- Can you think of other examples in which the statistical association between two variables might be expected to have a different sign than the causal relationship. How can incorrect conclusions about causal effects lead to harmful public policies?
- Explore relationships in this data among more recent health and demographic variables. For example, do the data suggest possible differences in the effect of private versus government expenditures on heath care as reflected by life expectancy at birth or life expectancy at 60? Are there differences by sex? Is the difference consistent in different countries? What might account for differences?