marPredict: Unbiased predictions for fixed variables
In MaciejDanko/marPredict: Marginal predictions for generalized mixed effect models
Description
Usage
Arguments
Value
Author(s)
Examples
Description
Unbiased predictions for fixed variables
Usage
1
2
3
4
5
Arguments
object
a object of class inherited from: glmer{lme4},
gamm{mgcv}, gam{mgcv}, gamm4{gamm4}, glmer.nb{lme4}, lmer{lme4}, lme{nlme},
glm{stats}, lm{stats}, or glm.nb{MASS}.
newdata
an optional data frame in which to look for variables with which to predict. If omitted, the fitted values are used.
type
the type of prediction required, either "link" or "response".
The default is on the scale of the linear predictors; the alternative "response" is on the scale of the response variable.
see also predict.glm{stats}.
ci.fit
a logical indicting if to bootstrap percentile confidence intervals.
alpha, n.sims
a significance level and a number of repetitions for the bootstrap method.
average.missing
a logical indicating if to perform population averaging over the variables missing from the newdata.
sim.missing
a logical indicating if to simulate uncertenity of the missing variables via the bootstrap method. Miningfull only when average.missing is TRUE.
sim.RE
a logical indicating if to simulate uncertenity of the random effects via the bootstrap method.
approx.RE
logical indicating if to use infinite-subjects approximation for random efect contribution to the marginalized model predictions.
method
a method of averaging over random effects.
The default "at.each.cat" estimate avarage random effects at each cobination of levels for categorical variable(s) present in the newdata;
the alternative "overall" estimate average random effects for total population.
use.offset
a logical indicating if to include averaged offset into model predictions.
as.rate
a logical indicating if to return predicitions as rates. Setting this TRUE makes sense only for count models (e.g., Poisson, neg.bin).
extended.boot.info
a logical indicating if to return median for bootstrap percentile method.
Value
a list with model predicitions with the following components:
fit
a list with biased (random effects excluded) and unbiased(averaged random effects included) predicitions.
CI.fit
a optional list with bootstraped confidence intervals for biased and unbiased predicitions.
boot
a optional list with bootstraped matrices for biased and unbiased predicitions.
offset
the offset used for model predictions.
Author(s)
Maciej J. Danko
Examples
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
#############################################################################
# EXAMPLE 1
#############################################################################
fit1 <- lme4::glmer(Y~X1*X2+(1|ID),family=poisson(), data=marPredict_e1)
DATA <- with(marPredict_e1, tapply(Y, data.frame(X1=X1,X2=X2), mean))
posX <- as.vector(barplot(DATA,beside=TRUE,ylim=c(0,1050),
ylab='Mean counts (response scale)'))
text(posX,-40,rep(rownames(DATA),3),xpd=TRUE)
axis(1,at=c(-5,100),labels = FALSE)
newdata1 <- expand.grid(X1=levels(marPredict_e1$X1),
X2=levels(marPredict_e1$X2))
# calculate the biased and the unbiased predictions
# on the response scale
cYA <- marPredict(object=fit1,
newdata=as.data.frame(newdata1),
type='response',
ci.fit = TRUE,
method = 'at.each.cat')
# using overall method
cYB <- marPredict(object=fit1,
newdata=as.data.frame(newdata1),
type='response',
ci.fit = TRUE,
method = 'overall')
# plot the predicited responses
lines(posX-0.3,cYA$fit$biased,pch=20,type='p',col='red',cex=1.5)
lines(posX+0.3,cYA$fit$unbiased,pch=20,type='p',
col=rgb(0,180,0,maxColorValue = 255),cex=1.5)
lines(posX,cYB$fit$unbiased,pch=20,type='p',col='darkorange',cex=1.5)
# plot the confidence intervals (the bootstrap percentile method)
for (i in seq_along(posX)) lines(c(posX[i],posX[i])+0.3,
cYA$CI.fit$unbiased[i,],
col=rgb(0,180,0,maxColorValue = 255))
for (i in seq_along(posX)) lines(c(posX[i],posX[i])-0.3,
cYA$CI.fit$biased[i,],
col='red')
for (i in seq_along(posX)) lines(c(posX[i],posX[i]),
cYB$CI.fit$unbiased[i,],
col='darkorange')
legend('topleft',c('Simulated data','RE excluded',
'RE averaged',
'RE category-averaged'),
pch=c(NA,19,19,19),pt.cex=c(1,1.5,1.5,1.5),
fill=c(gray.colors(3)[2],NA,NA,NA),
col=c(NA,2,'darkorange',rgb(0,180,0,maxColorValue = 255)),
bty='n', border=c('black',NA,NA,NA))
#############################################################################
# EXAMPLE 2
#############################################################################
fit2 <- gamm4::gamm4(Y~ X1 + s(C1), random = ~(1|ID),
family = poisson(), data = marPredict_e2)
newdata2 <- expand.grid(C1=seq(-0.5,0.5,0.05), X1=levels(marPredict_e2$X1))
cY2 <- marPredict(fit2, newdata2, type = 'response', ci.fit = TRUE)
ina <- which(newdata2$X1=='a'); inb <- which(newdata2$X1=='b')
inc <- which(newdata2$X1=='c')
xa <- newdata2$C1[ina]; xb <- newdata2$C1[inb]; xc <- newdata2$C1[inc]
CIplot(xa, cY2$fit$unbiased[ina], cY2$CI.fit$unbiased[ina,],
density=40, angle=0,first = TRUE,ylim=c(0,60),col='blue',lwd=4,
xaxt='s', yaxt='s', ylab='Mean counts (response scale)', xlab='C1')
CIplot(xa, cY2$fit$unbiased[inb], cY2$CI.fit$unbiased[inb,],
density=40, angle=0,col='red',lwd=4)
CIplot(xa, cY2$fit$unbiased[inc], cY2$CI.fit$unbiased[inc,],
density=40, angle=0,col=rgb(0,180,0,maxColorValue = 255),lwd=4)
CIplot(xa, cY2$fit$biased[ina], cY2$CI.fit$biased[ina,],
density=20, angle=90,col='blue',lty=2)
CIplot(xa, cY2$fit$biased[inb], cY2$CI.fit$biased[inb,],
density=20, angle=90,col='red',lty=2)
CIplot(xa, cY2$fit$biased[inc], cY2$CI.fit$biased[inc,],
density=20, angle=90,col=rgb(0,180,0,maxColorValue = 255),lty=2)
legend('top',c('biased','unbiased','a','b','c'), lty=c(2,1,1,1,1),
lwd=c(1,4,2,2,2),col=c(1,1,'blue','red',
rgb(0,180,0,maxColorValue = 255)),bty='n')
#############################################################################
# EXAMPLE 3
#############################################################################
ina <- which(marPredict_e3$X1=='a'); inb <- which(marPredict_e3$X1=='b')
inc <- which(marPredict_e3$X1=='c')
plot(marPredict_e3$C1[ina],binomial()$linkinv(marPredict_e3$linkY)[ina],
ylim=c(0,1), cex=0.2, col = 'blue')
lines(marPredict_e3$C1[inb],binomial()$linkinv(marPredict_e3$linkY)[inb],
type='p', cex=0.2, col = 'red')
lines(marPredict_e3$C1[inc],binomial()$linkinv(marPredict_e3$linkY)[inc],
type='p', cex=0.2, col = rgb(0,180,0,maxColorValue = 255))
fit3 <- glmer(Y ~ X1 + C1 + (1|ID), family=binomial(),
data=marPredict_e3, nAGQ = 20)
cY3 <- marPredict(fit3, newdata3, type = 'response', ci.fit = TRUE)
newdata3 <- expand.grid(C1=seq(-0.5,0.5,0.05), X1=levels(marPredict_e3$X1))
ina <- which(newdata3$X1=='a'); inb <- which(newdata3$X1=='b')
inc <- which(newdata3$X1=='c')
xa <- newdata3$C1[ina]; xb <- newdata3$C1[inb]; xc <- newdata3$C1[inc]
CIplot(xa, cY3$fit$unbiased[ina], cY3$CI.fit$unbiased[ina,],
density=40, angle=0,first = TRUE,ylim=c(0,1),col='blue',lwd=4,
xaxt='s', yaxt='s', ylab='Probabilyty (response scale)',
xlab='C1')
CIplot(xa, cY3$fit$unbiased[inb], cY3$CI.fit$unbiased[inb,],
density=40, angle=0,col='red',lwd=4)
CIplot(xa, cY3$fit$unbiased[inc], cY3$CI.fit$unbiased[inc,],
density=40, angle=0,col=rgb(0,180,0,maxColorValue = 255),
lwd=4)
CIplot(xa, cY3$fit$biased[ina], cY3$CI.fit$biased[ina,],
density=20, angle=90,col='blue',lty=2)
CIplot(xa, cY3$fit$biased[inb], cY3$CI.fit$biased[inb,],
density=20, angle=90,col='red',lty=2)
CIplot(xa, cY3$fit$biased[inc], cY3$CI.fit$biased[inc,],
density=20, angle=90,col=rgb(0,180,0,maxColorValue = 255),
lty=2)
legend('top',c('Biased pred. response', 'Unbiased pred. response',
'a','b','c'), lty=c(2,1,1,1,1), lwd=c(1,4,2,2,2),
col=c(1,1,'blue','red', rgb(0,180,0,maxColorValue = 255)), bty='n')
#############################################################################
# EXAMPLE 4
#############################################################################
#dontrun
#cat(1)
#donttest
cat(1)
#
MaciejDanko/marPredict documentation built on June 22, 2019, 12:42 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Usage Arguments Value Author(s) Examples
Description
Unbiased predictions for fixed variables
Usage
1 2 3 4 5 |
Arguments
object |
a object of class inherited from: |
newdata |
an optional data frame in which to look for variables with which to predict. If omitted, the fitted values are used. |
type |
the type of prediction required, either " |
ci.fit |
a logical indicting if to bootstrap percentile confidence intervals. |
alpha, n.sims |
a significance level and a number of repetitions for the bootstrap method. |
average.missing |
a logical indicating if to perform population averaging over the variables missing from the |
sim.missing |
a logical indicating if to simulate uncertenity of the missing variables via the bootstrap method. Miningfull only when |
sim.RE |
a logical indicating if to simulate uncertenity of the random effects via the bootstrap method. |
approx.RE |
logical indicating if to use infinite-subjects approximation for random efect contribution to the marginalized model predictions. |
method |
a method of averaging over random effects.
The default " |
use.offset |
a logical indicating if to include averaged offset into model predictions. |
as.rate |
a logical indicating if to return predicitions as rates. Setting this |
extended.boot.info |
a logical indicating if to return median for bootstrap percentile method. |
Value
a list with model predicitions with the following components:
fit |
a list with |
CI.fit |
a optional list with bootstraped confidence intervals for |
boot |
a optional list with bootstraped matrices for |
offset |
the offset used for model predictions. |
Author(s)
Maciej J. Danko
Examples
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 | #############################################################################
# EXAMPLE 1
#############################################################################
fit1 <- lme4::glmer(Y~X1*X2+(1|ID),family=poisson(), data=marPredict_e1)
DATA <- with(marPredict_e1, tapply(Y, data.frame(X1=X1,X2=X2), mean))
posX <- as.vector(barplot(DATA,beside=TRUE,ylim=c(0,1050),
ylab='Mean counts (response scale)'))
text(posX,-40,rep(rownames(DATA),3),xpd=TRUE)
axis(1,at=c(-5,100),labels = FALSE)
newdata1 <- expand.grid(X1=levels(marPredict_e1$X1),
X2=levels(marPredict_e1$X2))
# calculate the biased and the unbiased predictions
# on the response scale
cYA <- marPredict(object=fit1,
newdata=as.data.frame(newdata1),
type='response',
ci.fit = TRUE,
method = 'at.each.cat')
# using overall method
cYB <- marPredict(object=fit1,
newdata=as.data.frame(newdata1),
type='response',
ci.fit = TRUE,
method = 'overall')
# plot the predicited responses
lines(posX-0.3,cYA$fit$biased,pch=20,type='p',col='red',cex=1.5)
lines(posX+0.3,cYA$fit$unbiased,pch=20,type='p',
col=rgb(0,180,0,maxColorValue = 255),cex=1.5)
lines(posX,cYB$fit$unbiased,pch=20,type='p',col='darkorange',cex=1.5)
# plot the confidence intervals (the bootstrap percentile method)
for (i in seq_along(posX)) lines(c(posX[i],posX[i])+0.3,
cYA$CI.fit$unbiased[i,],
col=rgb(0,180,0,maxColorValue = 255))
for (i in seq_along(posX)) lines(c(posX[i],posX[i])-0.3,
cYA$CI.fit$biased[i,],
col='red')
for (i in seq_along(posX)) lines(c(posX[i],posX[i]),
cYB$CI.fit$unbiased[i,],
col='darkorange')
legend('topleft',c('Simulated data','RE excluded',
'RE averaged',
'RE category-averaged'),
pch=c(NA,19,19,19),pt.cex=c(1,1.5,1.5,1.5),
fill=c(gray.colors(3)[2],NA,NA,NA),
col=c(NA,2,'darkorange',rgb(0,180,0,maxColorValue = 255)),
bty='n', border=c('black',NA,NA,NA))
#############################################################################
# EXAMPLE 2
#############################################################################
fit2 <- gamm4::gamm4(Y~ X1 + s(C1), random = ~(1|ID),
family = poisson(), data = marPredict_e2)
newdata2 <- expand.grid(C1=seq(-0.5,0.5,0.05), X1=levels(marPredict_e2$X1))
cY2 <- marPredict(fit2, newdata2, type = 'response', ci.fit = TRUE)
ina <- which(newdata2$X1=='a'); inb <- which(newdata2$X1=='b')
inc <- which(newdata2$X1=='c')
xa <- newdata2$C1[ina]; xb <- newdata2$C1[inb]; xc <- newdata2$C1[inc]
CIplot(xa, cY2$fit$unbiased[ina], cY2$CI.fit$unbiased[ina,],
density=40, angle=0,first = TRUE,ylim=c(0,60),col='blue',lwd=4,
xaxt='s', yaxt='s', ylab='Mean counts (response scale)', xlab='C1')
CIplot(xa, cY2$fit$unbiased[inb], cY2$CI.fit$unbiased[inb,],
density=40, angle=0,col='red',lwd=4)
CIplot(xa, cY2$fit$unbiased[inc], cY2$CI.fit$unbiased[inc,],
density=40, angle=0,col=rgb(0,180,0,maxColorValue = 255),lwd=4)
CIplot(xa, cY2$fit$biased[ina], cY2$CI.fit$biased[ina,],
density=20, angle=90,col='blue',lty=2)
CIplot(xa, cY2$fit$biased[inb], cY2$CI.fit$biased[inb,],
density=20, angle=90,col='red',lty=2)
CIplot(xa, cY2$fit$biased[inc], cY2$CI.fit$biased[inc,],
density=20, angle=90,col=rgb(0,180,0,maxColorValue = 255),lty=2)
legend('top',c('biased','unbiased','a','b','c'), lty=c(2,1,1,1,1),
lwd=c(1,4,2,2,2),col=c(1,1,'blue','red',
rgb(0,180,0,maxColorValue = 255)),bty='n')
#############################################################################
# EXAMPLE 3
#############################################################################
ina <- which(marPredict_e3$X1=='a'); inb <- which(marPredict_e3$X1=='b')
inc <- which(marPredict_e3$X1=='c')
plot(marPredict_e3$C1[ina],binomial()$linkinv(marPredict_e3$linkY)[ina],
ylim=c(0,1), cex=0.2, col = 'blue')
lines(marPredict_e3$C1[inb],binomial()$linkinv(marPredict_e3$linkY)[inb],
type='p', cex=0.2, col = 'red')
lines(marPredict_e3$C1[inc],binomial()$linkinv(marPredict_e3$linkY)[inc],
type='p', cex=0.2, col = rgb(0,180,0,maxColorValue = 255))
fit3 <- glmer(Y ~ X1 + C1 + (1|ID), family=binomial(),
data=marPredict_e3, nAGQ = 20)
cY3 <- marPredict(fit3, newdata3, type = 'response', ci.fit = TRUE)
newdata3 <- expand.grid(C1=seq(-0.5,0.5,0.05), X1=levels(marPredict_e3$X1))
ina <- which(newdata3$X1=='a'); inb <- which(newdata3$X1=='b')
inc <- which(newdata3$X1=='c')
xa <- newdata3$C1[ina]; xb <- newdata3$C1[inb]; xc <- newdata3$C1[inc]
CIplot(xa, cY3$fit$unbiased[ina], cY3$CI.fit$unbiased[ina,],
density=40, angle=0,first = TRUE,ylim=c(0,1),col='blue',lwd=4,
xaxt='s', yaxt='s', ylab='Probabilyty (response scale)',
xlab='C1')
CIplot(xa, cY3$fit$unbiased[inb], cY3$CI.fit$unbiased[inb,],
density=40, angle=0,col='red',lwd=4)
CIplot(xa, cY3$fit$unbiased[inc], cY3$CI.fit$unbiased[inc,],
density=40, angle=0,col=rgb(0,180,0,maxColorValue = 255),
lwd=4)
CIplot(xa, cY3$fit$biased[ina], cY3$CI.fit$biased[ina,],
density=20, angle=90,col='blue',lty=2)
CIplot(xa, cY3$fit$biased[inb], cY3$CI.fit$biased[inb,],
density=20, angle=90,col='red',lty=2)
CIplot(xa, cY3$fit$biased[inc], cY3$CI.fit$biased[inc,],
density=20, angle=90,col=rgb(0,180,0,maxColorValue = 255),
lty=2)
legend('top',c('Biased pred. response', 'Unbiased pred. response',
'a','b','c'), lty=c(2,1,1,1,1), lwd=c(1,4,2,2,2),
col=c(1,1,'blue','red', rgb(0,180,0,maxColorValue = 255)), bty='n')
#############################################################################
# EXAMPLE 4
#############################################################################
#dontrun
#cat(1)
#donttest
cat(1)
#
|
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.