Hi folks, Warning: I don't know if the result I am getting makes sense, so this may be a statistics question.
The fitted values from my binomial lmer mixed model seem to consistently overestimate the cell means, and I don't know why. I assume I am doing something stupid. Below I include code, and a binary image of the data is available at this link: http://www.cas.muohio.edu/~stevenmh/tf.RdataBin This was done with `Matrix' version 0.995-10 and `lme4' version 0.995-2. and R v. 2.3.1 on a Mac, OS 10.4.6. The binomial model below ("mod") was reduced from a more complex one by first using AIC, BIC and LRT for "random" effects, and then relying on Helmert contrasts and AIC, BIC, and LRT to simplify fixed effects. Maybe this was wrong? > load("tf.RdataBin") > library(lme4) > options(contrasts=c("contr.helmert", "contr.poly")) > mod <- lmer(tfb ~ reg+nutrient+amd +reg:nutrient+ + (1|rack) + (1|popu) + (1|gen), data=dat.tb2, family=binomial, method="Laplace") > summary(mod) Generalized linear mixed model fit using Laplace Formula: tfb ~ reg + nutrient + amd + reg:nutrient + (1 | rack) + (1 | popu) + (1 | gen) Data: dat.tb2 Family: binomial(logit link) AIC BIC logLik deviance 402.53 446.64 -191.26 382.53 Random effects: Groups Name Variance Std.Dev. gen (Intercept) 0.385 0.621 popu (Intercept) 0.548 0.741 rack (Intercept) 0.401 0.633 number of obs: 609, groups: gen, 24; popu, 9; rack, 2 Estimated scale (compare to 1) 0.80656 Fixed effects: Estimate Std. Error z value Pr(>|z|) (Intercept) 2.391 0.574 4.17 3.1e-05 reg1 0.842 0.452 1.86 0.06252 reg2 0.800 0.241 3.32 0.00091 nutrient1 0.788 0.197 4.00 6.3e-05 amd1 -0.540 0.139 -3.88 0.00010 reg1:nutrient1 0.500 0.227 2.21 0.02734 reg2:nutrient1 -0.176 0.146 -1.21 0.22794 Correlation of Fixed Effects: (Intr) reg1 reg2 ntrnt1 amd1 rg1:n1 reg1 0.169 reg2 -0.066 -0.191 nutrient1 0.178 0.231 -0.034 amd1 -0.074 -0.044 -0.052 -0.078 reg1:ntrnt1 0.157 0.307 -0.180 0.562 -0.002 reg2:ntrnt1 -0.028 -0.154 0.236 0.141 0.033 -0.378 > X <- mod @ X > fitted <- X %*% fixef(mod) > > unlogitH <- function(x) {( 1 + exp(-x) )^-1} > (result <- data.frame(Raw.Data=with(dat.tb2, + tapply(tfb, list(reg:nutrient:amd), + mean ) ), + Fitted.Estimates=with(dat.tb2, + tapply(fitted, list(reg:nutrient:amd), + function(x) unlogitH(mean(x)) ) ) )) Raw.Data Fitted.Estimates SW:1:unclipped 0.50877 0.69520 SW:1:clipped 0.41304 0.43669 SW:8:unclipped 0.67273 0.85231 SW:8:clipped 0.52830 0.66233 NL:1:unclipped 0.88889 0.81887 NL:1:clipped 0.53571 0.60578 NL:8:unclipped 0.96552 0.98830 NL:8:clipped 0.96154 0.96635 SP:1:unclipped 0.98649 0.98361 SP:1:clipped 0.92537 0.95328 SP:8:unclipped 1.00000 0.99308 SP:8:clipped 0.95890 0.97992 > ### Perhaps the cell SP:8:clipped = 1.0 is messing up the fit? > pdf("RawAndFitted.pdf") > par(mar=c(8,3,2,2), las=2) > barplot(t(result), beside=TRUE ) > box(); title("Fractions of Plants Producing Fruits") > legend("topleft", c("Raw Data", "Fitted Values"), + fill=gray.colors(2), bty="n" ) > dev.off() quartz 2 > _ platform powerpc-apple-darwin8.6.0 arch powerpc os darwin8.6.0 system powerpc, darwin8.6.0 status major 2 minor 3.1 year 2006 month 06 day 01 svn rev 38247 language R version.string Version 2.3.1 (2006-06-01) > Dr. M. Hank H. Stevens, Assistant Professor 338 Pearson Hall Botany Department Miami University Oxford, OH 45056 Office: (513) 529-4206 Lab: (513) 529-4262 FAX: (513) 529-4243 http://www.cas.muohio.edu/~stevenmh/ http://www.muohio.edu/ecology/ http://www.muohio.edu/botany/ "E Pluribus Unum" ______________________________________________ [email protected] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
