(Updated: September 04 2024 19:54)
Get the lastest versions:
# devtools::install_github('gmonette/p3d')
# devtools::install_github('gmonette/spida2')
library(p3d)
Attaching package: 'p3d' The following object is masked from 'package:knitr':
spinlibrary(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, uvdata(Smoking2)
dd <- Smoking2
head(dd) X Country Cig.pAdult.pYear GDPpCapita Year HE.pCap HE.pcGDP LE.all
Albania 1 Albania 1201 4030 2011 569 6.8 77.41
Algeria 2 Algeria 577 5244 2011 437 5.4 74.02
Angola 3 Angola 397 5148 2011 183 3.3 38.20
Argentina 4 Argentina 1014 10941 2011 1062 7.4 76.95
Armenia 5 Armenia 2083 3305 2011 224 3.8 72.68
Australia 6 Australia 1130 60642 2011 3365 8.5 81.81
LE.male LE.female LE.all2 LE.male2 LE.female2 Continent
Albania 74.82 80.30 76.4 73.4 79.7 Europe
Algeria 72.35 75.77 72.3 70.9 73.7 Africa
Angola 37.24 39.22 NA NA NA Africa
Argentina 73.71 80.36 75.3 71.6 79.1 Americas
Armenia 69.06 76.81 72.0 68.4 75.1 Asia
Australia 79.40 84.35 81.2 78.9 83.6 Oceania
Region LE HE Smoking GDP
Albania Southern Europe 77.41 569 1201 4030
Algeria Northern Africa 74.02 437 577 5244
Angola Middle Africa 38.20 183 397 5148
Argentina South America 76.95 1062 1014 10941
Armenia Western Asia 72.68 224 2083 3305
Australia Australia and New Zealand 81.81 3365 1130 60642dim(dd) [1] 229 19dd$Cigarettes <- dd$Smoking/365.25
dd$Life <- dd$LE
dd$Health_Exp <- dd$HE
Init3d(cex=2)
Plot3d( Life ~ Cigarettes + Health_Exp | Continent, dd) Continent col
1 Africa blue
2 Americas green
3 Asia orange
4 Europe magenta
5 Oceania cyan Use left mouse to rotate, middle mouse (or scroll) to zoom, right mouse to change perspectivespinto()
Axes3d()
################
fg()
Id3d(pad = 1) Move display using left mouse button and middle button or scroll
Press right mouse button and drag a rectangle around points to be identified
Click right mouse button without selecting a point to terminate character(0)fit <- lm( Life ~ Cigarettes, dd)
summary(fit)
Call:
lm(formula = Life ~ Cigarettes, data = dd)
Residuals:
Min 1Q Median 3Q Max
-29.320 -4.521 1.167 5.510 14.560
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 64.8668 1.3066 49.644 < 2e-16 ***
Cigarettes 2.4410 0.4006 6.093 1.41e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8.389 on 119 degrees of freedom
(108 observations deleted due to missingness)
Multiple R-squared: 0.2378, Adjusted R-squared: 0.2314
F-statistic: 37.12 on 1 and 119 DF, p-value: 1.405e-08wald(fit) numDF denDF F-value p-value
2 119 4393.183 <.00001
Estimate Std.Error DF t-value p-value Lower 0.95 Upper 0.95
(Intercept) 64.866758 1.306644 119 49.643793 <.00001 62.279472 67.454043
Cigarettes 2.440957 0.400615 119 6.093026 <.00001 1.647699 3.234214# What does the regression look like?
Fit3d(fit, lwd = 3)
rglwidget()
But the relationship is not really linear
fitsq <- lm( Life ~ Cigarettes+I(Cigarettes^2), dd)
summary(fitsq)
Call:
lm(formula = Life ~ Cigarettes + I(Cigarettes^2), data = dd)
Residuals:
Min 1Q Median 3Q Max
-28.4770 -4.5343 0.6638 4.5436 15.2961
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 58.1288 1.6788 34.625 < 2e-16 ***
Cigarettes 8.9330 1.2166 7.343 2.93e-11 ***
I(Cigarettes^2) -0.9829 0.1760 -5.583 1.53e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7.492 on 118 degrees of freedom
(108 observations deleted due to missingness)
Multiple R-squared: 0.3971, Adjusted R-squared: 0.3868
F-statistic: 38.85 on 2 and 118 DF, p-value: 1.086e-13Fit3d(fitsq, lwd = 3, col = 'red') Error in if (use.groups) {: argument is of length zeroPop3d(4) Error in pop3d(): 'rgl_pop' failed for id 0# Id3d(pad=1)
fg()
# "Controlling" for Health $$
spinto(-90)
fit1 <- lm( Life ~ Cigarettes, dd)
summary(fit1)
Call:
lm(formula = Life ~ Cigarettes, data = dd)
Residuals:
Min 1Q Median 3Q Max
-29.320 -4.521 1.167 5.510 14.560
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 64.8668 1.3066 49.644 < 2e-16 ***
Cigarettes 2.4410 0.4006 6.093 1.41e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8.389 on 119 degrees of freedom
(108 observations deleted due to missingness)
Multiple R-squared: 0.2378, Adjusted R-squared: 0.2314
F-statistic: 37.12 on 1 and 119 DF, p-value: 1.405e-08Fit3d(fit1, alpha = .5) Error in if (use.groups) {: argument is of length zerofitlin <- lm( Life ~ Cigarettes
+ Health_Exp
,dd)
Fit3d(fitlin, col = 'cyan') Error in if (use.groups) {: argument is of length zero#
# Controlled for Health?
#
# - Have we unconfounded the effect of Smoking??
# - Hardly.
#
fith <- lm( Life ~ Cigarettes
+ Health_Exp
+ log( Health_Exp),dd)
summary(fith)
Call:
lm(formula = Life ~ Cigarettes + Health_Exp + log(Health_Exp),
data = dd)
Residuals:
Min 1Q Median 3Q Max
-27.8347 -1.9539 0.7817 3.4422 9.0891
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 33.9683970 4.2075265 8.073 7.80e-13 ***
Cigarettes 0.2115706 0.3608833 0.586 0.5589
Health_Exp -0.0014479 0.0007463 -1.940 0.0549 .
log(Health_Exp) 6.1620855 0.8439484 7.301 4.13e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.987 on 114 degrees of freedom
(111 observations deleted due to missingness)
Multiple R-squared: 0.5977, Adjusted R-squared: 0.5871
F-statistic: 56.46 on 3 and 114 DF, p-value: < 2.2e-16Fit3d(fith, col = 'red') Error in if (use.groups) {: argument is of length zeroApparently:
There is less confounding as we improve the model for the potential confounder, Health_Exp,
But we, presumably, are far from having found a reasonable model to estimate the causal effect of smoking: individually or ecologically.
EXERCISES: