Upload the R packages needed to illustate
library("MVA")
## Loading required package: HSAUR2
## Loading required package: tools
library("ellipse")
##
## Attaching package: 'ellipse'
## The following object is masked from 'package:graphics':
##
## pairs
set.seed(280875)
library("lattice")
lattice.options(default.theme =
function()
standard.theme("pdf", color = FALSE))
The data show life expectancy in years by country, age, and sex. The data come from Keyfitz and Flieger (1971) and relate to life expectancies in the 1960s.
d <- c(0.447,
0.422, 0.619,
0.435, 0.604, 0.583,
0.114, 0.068, 0.053, 0.115,
0.203, 0.146, 0.139, 0.258, 0.349,
0.091, 0.103, 0.110, 0.122, 0.209, 0.221,
0.082, 0.063, 0.066, 0.097, 0.321, 0.355, 0.201,
0.513, 0.445, 0.365, 0.482, 0.186, 0.315, 0.150, 0.154,
0.304, 0.318, 0.240, 0.368, 0.303, 0.377, 0.163, 0.219, 0.534,
0.245, 0.203, 0.183, 0.255, 0.272, 0.323, 0.310, 0.288, 0.301, 0.302,
0.101, 0.088, 0.074, 0.139, 0.279, 0.367, 0.232, 0.320, 0.204, 0.368, 0.340,
0.245, 0.199, 0.184, 0.293, 0.278, 0.545, 0.232, 0.314, 0.394, 0.467, 0.392, 0.511)
druguse <- diag(13) / 2
druguse[upper.tri(druguse)] <- d
druguse <- druguse + t(druguse)
rownames(druguse) <- colnames(druguse) <- c("cigarettes", "beer", "wine", "liquor", "cocaine",
"tranquillizers", "drug store medication", "heroin",
"marijuana", "hashish", "inhalants", "hallucinogenics", "amphetamine")
"life" <- structure(.Data = list(c(63., 34., 38., 59., 56., 62., 50., 65., 56., 69., 65., 64., 56., 60., 61., 49., 59., 63.,59., 65., 65., 64., 64., 67., 61., 68., 67., 65., 59., 58., 57.)
, c(51., 29., 30., 42., 38., 44., 39., 44., 46., 47., 48., 50., 44., 44., 45., 40., 42., 44., 44., 48., 48., 63.,
43., 45., 40., 46., 45., 46., 43., 44., 46.)
, c(30., 13., 17., 20., 18., 24., 20., 22., 24., 24., 26., 28., 25., 22., 22., 22., 22., 23., 24., 28., 26., 21.,
21., 23., 21., 23., 23., 24., 23., 24., 28.)
, c(13., 5., 7., 6., 7., 7., 7., 7., 11., 8., 9., 11., 10., 6., 8., 9., 6., 8., 8., 14., 9., 7., 6., 8., 10., 8.,
8., 9., 10., 9., 9.)
, c(67., 38., 38., 64., 62., 69., 55., 72., 63., 75., 68., 66., 61., 65., 65., 51., 61., 67., 63., 68., 67., 68.,
68., 74., 67., 75., 74., 71., 66., 62., 60.)
, c(54., 32., 34., 46., 46., 50., 43., 50., 54., 53., 50., 51., 48., 45., 49., 41., 43., 48., 46., 51., 49., 47.,
47., 51., 46., 52., 51., 51., 49., 47., 49.)
, c(34., 17., 20., 25., 25., 28., 23., 27., 33., 29., 27., 29., 27., 25., 27., 23., 22., 26., 25., 29., 27., 25.,
24., 28., 25., 29., 28., 28., 27., 25., 28.)
, c(15., 6., 7., 8., 10., 14., 8., 9., 19., 10., 10., 11., 12., 9., 10., 8., 7., 9., 8., 13., 10., 9., 8., 10., 11.,
10., 10., 10., 12., 10., 11.)
)
, class = "data.frame"
, names = c("m0", "m25", "m50", "m75", "w0", "w25", "w50", "w75")
, row.names = c("Algeria", "Cameroon", "Madagascar", "Mauritius", "Reunion", "Seychelles", "South Africa (C)", "South Africa (W)",
"Tunisia", "Canada", "Costa Rica", "Dominican Rep.", "El Salvador", "Greenland", "Grenada", "Guatemala",
"Honduras", "Jamaica", "Mexico", "Nicaragua", "Panama", "Trinidad (62)", "Trinidad (67)",
"United States (66)", "United States (NW66)", "United States (W66)", "United States (67)", "Argentina",
"Chile", "Colombia", "Ecuador")
)
toLatex(HSAURtable(life), pcol = 1, rownames = TRUE,
caption = "Life expectancies for different countries by age and gender.",
label = "ch:EFA:life:tab")
## \index{life data@\Robject{life} data}
## \begin{center}
## \begin{longtable}{l rrrrrrrr }
## \caption{\Robject{life} data. Life expectancies for different countries by age and gender. \label{ch:EFA:life:tab}}
## \\
## \hline
## & \Robject{m0} & \Robject{m25} & \Robject{m50} & \Robject{m75} & \Robject{w0} & \Robject{w25} & \Robject{w50} & \Robject{w75} \\ \hline
## \endfirsthead
## \caption[]{\Robject{life} data (continued).} \\
## \hline
## & \Robject{m0} & \Robject{m25} & \Robject{m50} & \Robject{m75} & \Robject{w0} & \Robject{w25} & \Robject{w50} & \Robject{w75} \\ \hline
## \endhead
## Algeria & 63 & 51 & 30 & 13 & 67 & 54 & 34 & 15 \\
## Cameroon & 34 & 29 & 13 & 5 & 38 & 32 & 17 & 6 \\
## Madagascar & 38 & 30 & 17 & 7 & 38 & 34 & 20 & 7 \\
## Mauritius & 59 & 42 & 20 & 6 & 64 & 46 & 25 & 8 \\
## Reunion & 56 & 38 & 18 & 7 & 62 & 46 & 25 & 10 \\
## Seychelles & 62 & 44 & 24 & 7 & 69 & 50 & 28 & 14 \\
## South Africa (C) & 50 & 39 & 20 & 7 & 55 & 43 & 23 & 8 \\
## South Africa (W) & 65 & 44 & 22 & 7 & 72 & 50 & 27 & 9 \\
## Tunisia & 56 & 46 & 24 & 11 & 63 & 54 & 33 & 19 \\
## Canada & 69 & 47 & 24 & 8 & 75 & 53 & 29 & 10 \\
## Costa Rica & 65 & 48 & 26 & 9 & 68 & 50 & 27 & 10 \\
## Dominican Rep. & 64 & 50 & 28 & 11 & 66 & 51 & 29 & 11 \\
## El Salvador & 56 & 44 & 25 & 10 & 61 & 48 & 27 & 12 \\
## Greenland & 60 & 44 & 22 & 6 & 65 & 45 & 25 & 9 \\
## Grenada & 61 & 45 & 22 & 8 & 65 & 49 & 27 & 10 \\
## Guatemala & 49 & 40 & 22 & 9 & 51 & 41 & 23 & 8 \\
## Honduras & 59 & 42 & 22 & 6 & 61 & 43 & 22 & 7 \\
## Jamaica & 63 & 44 & 23 & 8 & 67 & 48 & 26 & 9 \\
## Mexico & 59 & 44 & 24 & 8 & 63 & 46 & 25 & 8 \\
## Nicaragua & 65 & 48 & 28 & 14 & 68 & 51 & 29 & 13 \\
## Panama & 65 & 48 & 26 & 9 & 67 & 49 & 27 & 10 \\
## Trinidad (62) & 64 & 63 & 21 & 7 & 68 & 47 & 25 & 9 \\
## Trinidad (67) & 64 & 43 & 21 & 6 & 68 & 47 & 24 & 8 \\
## United States (66) & 67 & 45 & 23 & 8 & 74 & 51 & 28 & 10 \\
## United States (NW66) & 61 & 40 & 21 & 10 & 67 & 46 & 25 & 11 \\
## United States (W66) & 68 & 46 & 23 & 8 & 75 & 52 & 29 & 10 \\
## United States (67) & 67 & 45 & 23 & 8 & 74 & 51 & 28 & 10 \\
## Argentina & 65 & 46 & 24 & 9 & 71 & 51 & 28 & 10 \\
## Chile & 59 & 43 & 23 & 10 & 66 & 49 & 27 & 12 \\
## Colombia & 58 & 44 & 24 & 9 & 62 & 47 & 25 & 10 \\
## Ecuador & 57 & 46 & 28 & 9 & 60 & 49 & 28 & 11 \\
## \hline
## \end{longtable}
## \end{center}
sapply(1:3, function(f)
factanal(life, factors = f, method ="mle")$PVAL)
## objective objective objective
## 1.879555e-24 1.911514e-05 4.578204e-01
factanal(life, factors = 3, method ="mle")
##
## Call:
## factanal(x = life, factors = 3, method = "mle")
##
## Uniquenesses:
## m0 m25 m50 m75 w0 w25 w50 w75
## 0.005 0.362 0.066 0.288 0.005 0.011 0.020 0.146
##
## Loadings:
## Factor1 Factor2 Factor3
## m0 0.964 0.122 0.226
## m25 0.646 0.169 0.438
## m50 0.430 0.354 0.790
## m75 0.525 0.656
## w0 0.970 0.217
## w25 0.764 0.556 0.310
## w50 0.536 0.729 0.401
## w75 0.156 0.867 0.280
##
## Factor1 Factor2 Factor3
## SS loadings 3.375 2.082 1.640
## Proportion Var 0.422 0.260 0.205
## Cumulative Var 0.422 0.682 0.887
##
## Test of the hypothesis that 3 factors are sufficient.
## The chi square statistic is 6.73 on 7 degrees of freedom.
## The p-value is 0.458
(scores <- factanal(life, factors = 3, method = "mle",
scores = "regression")$scores)
## Factor1 Factor2 Factor3
## Algeria -0.258062561 1.90095771 1.91581631
## Cameroon -2.782495791 -0.72340014 -1.84772224
## Madagascar -2.806428187 -0.81158820 -0.01210318
## Mauritius 0.141004934 -0.29028454 -0.85862443
## Reunion -0.196352142 0.47429917 -1.55046466
## Seychelles 0.367371307 0.82902375 -0.55214085
## South Africa (C) -1.028567629 -0.08065792 -0.65421971
## South Africa (W) 0.946193522 0.06400408 -0.91995289
## Tunisia -0.862493550 3.59177195 -0.36442148
## Canada 1.245304248 0.29564122 -0.27342781
## Costa Rica 0.508736247 -0.50500435 1.01328707
## Dominican Rep. 0.106044085 0.01111171 1.83871599
## El Salvador -0.608155779 0.65100820 0.48836431
## Greenland 0.235114220 -0.69123901 -0.38558654
## Grenada 0.132008172 0.25241049 -0.15220645
## Guatemala -1.450336359 -0.67765804 0.65911906
## Honduras 0.043253249 -1.85175707 0.30633182
## Jamaica 0.462124701 -0.51918493 0.08032855
## Mexico -0.052332675 -0.72020002 0.44417800
## Nicaragua 0.268974443 0.08407227 1.70568388
## Panama 0.442333434 -0.73778272 1.25218728
## Trinidad (62) 0.711367053 -0.95989475 -0.21545329
## Trinidad (67) 0.787286051 -1.10729029 -0.51958264
## United States (66) 1.128331259 0.16389896 -0.68177046
## United States (NW66) 0.400058903 -0.36230253 -0.74299137
## United States (W66) 1.214345385 0.40877239 -0.69225320
## United States (67) 1.128331259 0.16389896 -0.68177046
## Argentina 0.731344988 0.24811968 -0.12817725
## Chile 0.009751528 0.75222637 -0.49198911
## Colombia -0.240602517 -0.29543613 0.42919600
## Ecuador -0.723451797 0.44246371 1.59164974
cex <- 0.8
plot(scores[,1], scores[,2], type = "n", xlab = "Factor 1", ylab = "Factor 2")
text(scores[,1], scores[,2], abbreviate(rownames(life), 5), cex = cex)
plot(scores[,1], scores[,3], type = "n", xlab = "Factor 1", ylab = "Factor 3")
text(scores[,1], scores[,3], abbreviate(rownames(life), 5), cex = cex)
plot(scores[,2], scores[,3], type = "n", xlab = "Factor 2", ylab = "Factor 3")
text(scores[,2], scores[,3], abbreviate(rownames(life), 5), cex = cex)
The majority of adult and adolescent Americans regularly use psychoactive substances during an increasing proportion of their lifetimes. Various forms of licit and illicit psychoactive substance use are prevalent, suggesting that patterns of psychoactive substance taking are a major part of the individual’s behavioural repertory and have pervasive implications for the performance of other behaviours. In an investigation of these phenomena, Huba, Wingard, and Bentler (1981) collected data on drug usage rates for 1634 students in the seventh to ninth grades in 11 schools in the greater metropolitan area of Los Angeles. Each participant completed a questionnaire about the number of times a particular substance had ever been used
ord <- order.dendrogram(as.dendrogram(hclust(dist(druguse))))
panel.corrgram <-
function(x, y, z, subscripts, at,
level = 0.9, label = FALSE, ...)
{
require("ellipse", quietly = TRUE)
x <- as.numeric(x)[subscripts]
y <- as.numeric(y)[subscripts]
z <- as.numeric(z)[subscripts]
zcol <- level.colors(z, at = at, col.regions = grey.colors, ...)
for (i in seq(along = z)) {
ell <- ellipse(z[i], level = level, npoints = 50,
scale = c(.2, .2), centre = c(x[i], y[i]))
panel.polygon(ell, col = zcol[i], border = zcol[i], ...)
}
if (label)
panel.text(x = x, y = y, lab = 100 * round(z, 2), cex = 0.8,
col = ifelse(z < 0, "white", "black"))
}
print(levelplot(druguse[ord, ord], at = do.breaks(c(-1.01, 1.01), 20),
xlab = NULL, ylab = NULL, colorkey = list(space = "top"),
scales = list(x = list(rot = 90)),
panel = panel.corrgram, label = TRUE))
sapply(1:6, function(nf)
factanal(covmat = druguse, factors = nf,
method = "mle", n.obs = 1634)$PVAL)
## objective objective objective objective objective objective
## 0.000000e+00 9.786000e-70 7.363910e-28 1.794578e-11 3.891743e-06 9.752967e-02
(factanal(covmat = druguse, factors = 6,
method = "mle", n.obs = 1634))
##
## Call:
## factanal(factors = 6, covmat = druguse, n.obs = 1634, method = "mle")
##
## Uniquenesses:
## cigarettes beer wine
## 0.563 0.368 0.374
## liquor cocaine tranquillizers
## 0.412 0.681 0.522
## drug store medication heroin marijuana
## 0.785 0.669 0.318
## hashish inhalants hallucinogenics
## 0.005 0.541 0.620
## amphetamine
## 0.005
##
## Loadings:
## Factor1 Factor2 Factor3 Factor4 Factor5 Factor6
## cigarettes 0.494 0.407 0.110
## beer 0.776 0.112
## wine 0.786
## liquor 0.720 0.121 0.103 0.115 0.160
## cocaine 0.519 0.132 0.158
## tranquillizers 0.130 0.564 0.321 0.105 0.143
## drug store medication 0.255 0.372
## heroin 0.532 0.101 0.190
## marijuana 0.429 0.158 0.152 0.259 0.609 0.110
## hashish 0.244 0.276 0.186 0.881 0.194 0.100
## inhalants 0.166 0.308 0.150 0.140 0.537
## hallucinogenics 0.387 0.335 0.186 0.288
## amphetamine 0.151 0.336 0.886 0.145 0.137 0.187
##
## Factor1 Factor2 Factor3 Factor4 Factor5 Factor6
## SS loadings 2.301 1.415 1.116 0.964 0.676 0.666
## Proportion Var 0.177 0.109 0.086 0.074 0.052 0.051
## Cumulative Var 0.177 0.286 0.372 0.446 0.498 0.549
##
## Test of the hypothesis that 6 factors are sufficient.
## The chi square statistic is 22.41 on 15 degrees of freedom.
## The p-value is 0.0975
pfun <- function(nf) {
fa <- factanal(covmat = druguse, factors = nf,
method = "mle", n.obs = 1634)
est <- tcrossprod(fa$loadings) + diag(fa$uniquenesses)
ret <- round(druguse - est, 3)
colnames(ret) <- rownames(ret) <-
abbreviate(rownames(ret), 3)
ret
}
pfun(6)
## cgr ber win lqr ccn trn dsm hrn mrj hsh inh
## cgr 0.000 -0.001 0.014 -0.018 0.010 0.001 -0.020 -0.004 0.001 0 0.010
## ber -0.001 0.000 -0.002 0.004 0.004 -0.011 -0.001 0.007 0.002 0 -0.004
## win 0.014 -0.002 0.000 -0.001 -0.001 -0.005 0.008 0.008 -0.004 0 -0.007
## lqr -0.018 0.004 -0.001 0.000 -0.008 0.021 -0.006 -0.018 0.003 0 0.012
## ccn 0.010 0.004 -0.001 -0.008 0.000 0.000 0.008 0.004 -0.004 0 -0.003
## trn 0.001 -0.011 -0.005 0.021 0.000 0.000 0.006 -0.004 -0.004 0 0.002
## dsm -0.020 -0.001 0.008 -0.006 0.008 0.006 0.000 -0.015 0.008 0 0.004
## hrn -0.004 0.007 0.008 -0.018 0.004 -0.004 -0.015 0.000 0.006 0 -0.002
## mrj 0.001 0.002 -0.004 0.003 -0.004 -0.004 0.008 0.006 0.000 0 -0.006
## hsh 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0 0.000
## inh 0.010 -0.004 -0.007 0.012 -0.003 0.002 0.004 -0.002 -0.006 0 0.000
## hll -0.005 0.005 -0.001 -0.005 -0.008 -0.008 -0.002 0.020 0.003 0 -0.002
## amp 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0 0.000
## hll amp
## cgr -0.005 0
## ber 0.005 0
## win -0.001 0
## lqr -0.005 0
## ccn -0.008 0
## trn -0.008 0
## dsm -0.002 0
## hrn 0.020 0
## mrj 0.003 0
## hsh 0.000 0
## inh -0.002 0
## hll 0.000 0
## amp 0.000 0
pfun(3)
## cgr ber win lqr ccn trn dsm hrn mrj hsh
## cgr 0.000 -0.001 0.009 -0.013 0.011 0.009 -0.011 -0.004 0.003 -0.027
## ber -0.001 0.000 -0.002 0.002 0.002 -0.014 0.000 0.005 -0.001 0.019
## win 0.009 -0.002 0.000 0.000 -0.002 -0.004 0.012 0.013 0.001 -0.017
## lqr -0.013 0.002 0.000 0.000 -0.008 0.024 -0.017 -0.020 -0.001 0.014
## ccn 0.011 0.002 -0.002 -0.008 0.000 0.031 0.038 0.082 -0.002 0.041
## trn 0.009 -0.014 -0.004 0.024 0.031 0.000 -0.021 0.026 -0.002 -0.016
## dsm -0.011 0.000 0.012 -0.017 0.038 -0.021 0.000 0.021 0.007 -0.040
## hrn -0.004 0.005 0.013 -0.020 0.082 0.026 0.021 0.000 0.006 -0.035
## mrj 0.003 -0.001 0.001 -0.001 -0.002 -0.002 0.007 0.006 0.000 0.001
## hsh -0.027 0.019 -0.017 0.014 0.041 -0.016 -0.040 -0.035 0.001 0.000
## inh 0.039 -0.002 -0.007 -0.002 0.023 -0.038 0.113 0.031 0.003 -0.035
## hll -0.017 0.009 0.004 -0.015 -0.030 -0.058 0.000 -0.005 -0.002 0.034
## amp 0.002 -0.007 0.002 0.006 -0.075 0.044 -0.038 -0.049 -0.002 0.010
## inh hll amp
## cgr 0.039 -0.017 0.002
## ber -0.002 0.009 -0.007
## win -0.007 0.004 0.002
## lqr -0.002 -0.015 0.006
## ccn 0.023 -0.030 -0.075
## trn -0.038 -0.058 0.044
## dsm 0.113 0.000 -0.038
## hrn 0.031 -0.005 -0.049
## mrj 0.003 -0.002 -0.002
## hsh -0.035 0.034 0.010
## inh 0.000 0.007 -0.015
## hll 0.007 0.000 0.041
## amp -0.015 0.041 0.000
pfun(4)
## cgr ber win lqr ccn trn dsm hrn mrj hsh
## cgr 0.000 -0.001 0.008 -0.012 0.009 0.008 -0.015 -0.007 0.001 -0.023
## ber -0.001 0.000 -0.001 0.001 0.000 -0.016 -0.002 0.003 -0.001 0.018
## win 0.008 -0.001 0.000 0.000 -0.001 -0.005 0.012 0.014 0.001 -0.020
## lqr -0.012 0.001 0.000 0.000 -0.004 0.029 -0.015 -0.015 -0.001 0.018
## ccn 0.009 0.000 -0.001 -0.004 0.000 0.024 -0.014 0.007 -0.003 0.035
## trn 0.008 -0.016 -0.005 0.029 0.024 0.000 -0.020 0.027 -0.001 0.001
## dsm -0.015 -0.002 0.012 -0.015 -0.014 -0.020 0.000 -0.018 0.003 -0.042
## hrn -0.007 0.003 0.014 -0.015 0.007 0.027 -0.018 0.000 0.003 -0.037
## mrj 0.001 -0.001 0.001 -0.001 -0.003 -0.001 0.003 0.003 0.000 0.000
## hsh -0.023 0.018 -0.020 0.018 0.035 0.001 -0.042 -0.037 0.000 0.000
## inh 0.037 -0.005 -0.008 0.001 -0.022 -0.032 0.090 -0.001 0.001 -0.031
## hll -0.020 0.006 0.001 -0.010 -0.028 -0.028 0.008 0.005 -0.002 0.055
## amp 0.000 0.000 0.000 -0.001 0.000 0.001 0.000 0.000 0.000 -0.001
## inh hll amp
## cgr 0.037 -0.020 0.000
## ber -0.005 0.006 0.000
## win -0.008 0.001 0.000
## lqr 0.001 -0.010 -0.001
## ccn -0.022 -0.028 0.000
## trn -0.032 -0.028 0.001
## dsm 0.090 0.008 0.000
## hrn -0.001 0.005 0.000
## mrj 0.001 -0.002 0.000
## hsh -0.031 0.055 -0.001
## inh 0.000 0.021 0.000
## hll 0.021 0.000 0.000
## amp 0.000 0.000 0.000