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
|
densityplot.lambdaReg <- function(x, by = "column", col,
xlim, ylim,
main = "", sub = NULL, xlab,
ylab, lty = par("lty"), lwd = par("lwd"), ...) {
readpars <- par(no.readonly = TRUE)
idx <- dimnames(x$lambda)
lidx <- sapply(idx, length)
names(lidx) <- names(idx) <- c("rows", "columns")
dens <- array(NA, dim = lidx, dimnames = idx)
for (i in idx[[2]]) {
for (j in idx[[1]])
dens[j,i] <- dnorm(x$lambda[j,i], mean = x$lambda[j,i],
sd = x$se[j,i])
}
if (missing(xlim)) {
tmpL <- x$lambda - 3*x$se
tmpH <- x$lambda + 3*x$se
xlim <- c(min(tmpL), max(tmpH))
}
if (by == "row") {
par(mfrow = c(lidx[1], 1), ...)
if (missing(ylim)) {
ylim <- apply(dens, 1, max)
}
if (missing(xlab)) {
xlab <- paste("Percentage", unlist(idx[[1]]))
}
else if (length(xlab) != length(idx[[1]])) {
warning(paste("xlab needs to be", lidx[1], "long"))
xlab <- rep(xlab, lidx[1])[1:lidx[1]]
}
if (missing(col)) {
col <- rainbow(lidx[2])
}
else {
if (length(col) < lidx[2])
col <- rep(col, idx[[2]])[1:lidx[2]]
if (length(col) > lidx[2])
col <- col[1:lidx[2]]
}
names(col) <- idx[[2]]
names(xlab) <- idx[[1]]
for (ii in idx[[1]]) {
x1 <- seq(tmpL[ii,1], tmpH[ii,1], by = 0.001)
d1 <- dnorm(x1, mean = x$lambda[ii,1], sd = x$se[ii,1])
plot(x1, d1, type = "l", xlim = xlim, col = col[1],
ylim = c(0, max(ylim[ii])), main = main, sub = sub,
xlab = xlab[ii], ylab = "Density", lty = lty,
lwd = lwd)
for (jj in idx[[2]][2:lidx[2]]) {
x1 <- seq(tmpL[ii,jj], tmpH[ii,jj], by = 0.001)
d1 <- dnorm(x1, mean = x$lambda[ii,jj], sd = x$se[ii,jj])
lines(x1, d1, lty = lty, lwd = lwd, col = col[jj])
}
abline(v = 0, col = "grey50")
abline(v = 1, col = "grey50")
}
}
if (by == "column") {
par(mfrow = c(lidx[2], 1), ...)
if (missing(ylim)) {
ylim <- apply(dens, 2, max)
}
if (missing(xlab)) {
xlab <- paste("Percentage", unlist(idx[[2]]))
}
else {
if (length(xlab) != lidx[1]) {
warning(paste("xlab needs to be", lidx[1], "long"))
xlab <- rep(xlab, lidx[2])[1:lidx[2]]
}
}
if (missing(col)) {
col <- rainbow(lidx[1])
}
else {
if (length(col) < lidx[1])
col <- rep(col, idx[[1]])[1:lidx[1]]
if (length(col) > lidx[1])
col <- col[1:lidx[1]]
}
names(col) <- idx[[1]]
names(xlab) <- idx[[2]]
for (ii in idx[[2]]) {
x1 <- seq(tmpL[1,ii], tmpH[1,ii], by = 0.001)
d1 <- dnorm(x1, mean = x$lambda[1,ii], sd = x$se[1,ii])
plot(x1, d1, type = "l", xlim = xlim, col = col[1],
ylim = c(0, max(ylim[ii])), main = main, sub = sub,
xlab = xlab[ii], ylab = "Density", lty = lty,
lwd = lwd)
for (jj in idx[[1]][2:lidx[1]]) {
x1 <- seq(tmpL[jj,ii], tmpH[jj,ii], by = 0.001)
d1 <- dnorm(x1, mean = x$lambda[jj,ii], sd = x$se[jj,ii])
lines(x1, d1, lty = lty, lwd = lwd, col = col[jj])
}
abline(v = 0, col = "grey50")
abline(v = 1, col = "grey50")
}
}
return(invisible(x))
par(readpars)
}
|
