1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
|
### R code from vignette source 'sparseModels.Rnw'
###################################################
### code chunk number 1: preliminaries
###################################################
options(width=75)
library(grDevices) # for R_DEFAULT_PACKAGES=NULL
library(stats) # ditto
library(utils) # ditto
###################################################
### code chunk number 2: ex1
###################################################
(ff <- factor(strsplit("statistics_is_a_task", "")[[1]], levels=c("_",letters)))
factor(ff) # drops the levels that do not occur
f1 <- ff[, drop=TRUE] # the same, more transparently
###################################################
### code chunk number 3: ex1.2
###################################################
levels(f1)[match(c("c","k"), levels(f1))] <- "ck"
library(Matrix)
Matrix(contrasts(f1)) # "treatment" contrasts by default -- level "_" = baseline
Matrix(contrasts(C(f1, sum)))
Matrix(contrasts(C(f1, helmert)), sparse=TRUE) # S-plus default; much less sparse
###################################################
### code chunk number 4: as_factor_sparse
###################################################
as(f1, "sparseMatrix")
###################################################
### code chunk number 5: contrasts_sub
###################################################
printSpMatrix( t( Matrix(contrasts(f1))[as.character(f1) ,] ),
col.names=TRUE)
###################################################
### code chunk number 6: ex1-model.matrix
###################################################
t( Matrix(model.matrix(~ 0+ f1))) # model with*OUT* intercept
###################################################
### code chunk number 7: chickwts-ex
###################################################
data(chickwts, package = "datasets")
str(chickwts)# a standard R data set, 71 x 2
x.feed <- as(chickwts$feed, "sparseMatrix")
x.feed[ , (1:72)[c(TRUE,FALSE,FALSE)]] ## every 3rd column:
###################################################
### code chunk number 8: warpbreaks-data
###################################################
data(warpbreaks, package = "datasets") # a standard R data set
str(warpbreaks) # 2 x 3 (x 9) balanced two-way with 9 replicates:
xtabs(~ wool + tension, data = warpbreaks)
###################################################
### code chunk number 9: modMat-warpbreaks
###################################################
tmm <- with(warpbreaks,
rbind(as(tension, "sparseMatrix"),
as(wool, "sparseMatrix")[-1,,drop=FALSE]))
print( image(tmm) ) # print(.) the lattice object
###################################################
### code chunk number 10: morley-data
###################################################
data(morley, package = "datasets") # a standard R data set
morley$Expt <- factor(morley$Expt)
morley$Run <- factor(morley$Run)
str(morley)
t.mm <- with(morley,
rbind(as(Expt, "sparseMatrix"),
as(Run, "sparseMatrix")[-1,]))
print( image(t.mm) ) # print(.) the lattice object
###################################################
### code chunk number 11: npk_ex
###################################################
data(npk, package = "MASS")
npk.mf <- model.frame(yield ~ block + N*P*K, data = npk)
## str(npk.mf) # the data frame + "terms" attribute
m.npk <- model.matrix(attr(npk.mf, "terms"), data = npk)
class(M.npk <- Matrix(m.npk))
dim(M.npk)# 24 x 13 sparse Matrix
t(M.npk) # easier to display, column names readably displayed as row.names(t(.))
###################################################
### code chunk number 12: aov-large-ex
###################################################
id <- factor(1:20)
a <- factor(1:2)
b <- factor(1:2)
d <- factor(1:1500)
aDat <- expand.grid(id=id, a=a, b=b, d=d)
aDat$y <- rnorm(length(aDat[, 1])) # generate some random DV data
dim(aDat) # 120'000 x 5 (120'000 = 2*2*1500 * 20 = 6000 * 20)
###################################################
### code chunk number 13: aov-ex-X-sparse
###################################################
d2 <- factor(1:150) # 10 times smaller
tmp2 <- expand.grid(id=id, a=a, b=b, d=d2)
dim(tmp2)
dim(mm <- model.matrix( ~ a*b*d, data=tmp2))
## is 100 times smaller than original example
class(smm <- Matrix(mm)) # automatically coerced to sparse
round(object.size(mm) / object.size(smm), 1)
###################################################
### code chunk number 14: X-sparse-image (eval = FALSE)
###################################################
## image(t(smm), aspect = 1/3, lwd=0, col.regions = "red")
###################################################
### code chunk number 15: X-sparse-image-fake
###################################################
png("sparseModels-X-sparse-image.png", width=6, height=3,
units='in', res=150)
print(
image(t(smm), aspect = 1/3, lwd=0, col.regions = "red")
)
dev.off()
###################################################
### code chunk number 16: X-sparse-mult
###################################################
x <- 1:600
system.time(y <- smm %*% x) ## sparse is much faster
system.time(y. <- mm %*% x) ## than dense
identical(as.matrix(y), y.) ## TRUE
###################################################
### code chunk number 17: sessionInfo
###################################################
toLatex(sessionInfo())
|
