[R-sig-ME] Predicted density function from mixed model?

[Hugh Morgan]

Dear All,

I want to plot the predicted probability density function derived from a model fitted with lme or gls next to the actual density function. I hope to be able to use this to present to the user, having the double function of helping to explain what statistics I am performing and some estimation of the model fit. I can get the actual density distribution (using the function density(x, ...)) but is there an easy way of getting the one predicted by the model? I expect I could derive it in code, by feeding the actual numbers for the various parameters used in the model, but I am hoping someone has done it for me.


For an example, using the attached csv file:


library(nlme)

# Note, 8.5k test file, I have also tried attaching it
testDataset <- read.csv("http://sanger_impc.drivehq.com/exampleDataset.csv")

model <- lme(Weight~Genotype+Gender+Genotype*Gender, random=~1|Batch, testDataset, na.action="na.omit", method="ML")


I get the density distribution with:


testDensity <- density(testDataset$Weight)


And then plot it


plot(testDensity)


and add the predicted density:


# Does not work

plot(density(model))


In the end of the day I would like to split the plots into control and test group, like :


plot(density(subset(testDataset, Genotype=="mutant")$Weight))

plot(density(subset(testDataset, Genotype=="control")$Weight))


I really appreciate any help, and I will understand if the answer is I have to write it myself.


Thanks in advance,


Hugh



This email may have a PROTECTIVE MARKING, for an explanation please see: http://www.mrc.ac.uk/About/Informationandstandards/Documentmarking/index.htm


This email may have a PROTECTIVE MARKING, for an explanation please see:
http://www.mrc.ac.uk/About/Informationandstandards/Documentmarking/index.htm