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

[Hugh Morgan]

Thanks for this, but predict(model, newdata=nd, level=1) seems to be exactly the same as predict(model, newdata=nd, level=0).  predict(model, newdata=nd, level=0:1) (from ?predict.lme) just adds more columns (1 replicating the batch column, 2 with the same values you get with either level=0 or level=1.

Hugh

________________________________________
From: r-sig-mixed-models-bounces at r-project.org [r-sig-mixed-models-bounces at r-project.org] on behalf of Alan Haynes [aghaynes at gmail.com]
Sent: Tuesday, February 19, 2013 9:39 AM
To: r-sig-mixed-models at r-project.org
Subject: [R-sig-ME] Fwd:  Predicted density function from mixed model?

And to the list...

Hi Hugh,

If you want the individual level predictions then you could include the
individual (Batch) in the nd dataframe and use

predict(model, newdata=nd, level=1)

that might be more meaningful than using rnorm. If i remember correctly,
you can get SE estimates for the values too...see ?predict.lme

HTH

Alan

--------------------------------------------------
Email: aghaynes at gmail.com
Mobile: +41763389128
Skype: aghaynes


On 19 February 2013 09:38, Hugh Morgan <h.morgan at har.mrc.ac.uk> wrote:

>  Hi Alan,
>
> Thank you very much for this.  The pred function is what I was looking
> for.  I THINK what I actually want I can get with:
>
> > summary(model)
> Linear mixed-effects model fit by maximum likelihood
> ...
> Random effects:
>  Formula: ~1 | Batch
>         (Intercept)  Residual
> StdDev:  0.09476296 0.1319891
>
> Fixed effects: Weight ~ Genotype + Gender + Genotype * Gender
>                               Value  Std.Error  DF   t-value p-value
> (Intercept)               0.5628385 0.02804277 159 20.070716  0.0000
> ...
>
> Taking:
> Random effects StdDev:  0.09476296
> Fixed effects Intercept Std.Error: 0.02804277
>
> I can then randomly generate these effects:
>
> nd$batchEffect <- rnorm(nrow(nd), mean = 0, sd = 0.09476296)
> nd$randomEffect <- rnorm(nrow(nd), mean = 0, sd = 0.02804277)
>
> And add them to the single values generated by pred:
>
> nd$predicedDistribution <- (nd$pred + nd$batchEffect + nd$randomEffect)
>
> Plotting these gives something like the real distribution, so it is
> possible I am right:
>
> plot(density(subset(nd, Genotype=="mutant")$predicedDistribution))
> plot(density(subset(nd, Genotype=="control")$predicedDistribution))
>
>
>  Incidentally, for your model, Weight ~ Genotype + Gender + Genotype*Gender is identical to Weight ~ Genotype*Gender because the * includes both main effects and interaction. If you want to specify just an interaction use : (e.g. Genotype:Gender)
>
>
> Thank you for this.  As is obvious I am just learning this.
>
> Hugh
>
>  ------------------------------
> *From:* Alan Haynes [aghaynes at gmail.com]
> *Sent:* Tuesday, February 19, 2013 7:27 AM
> *To:* Hugh Morgan
> *Cc:* r-sig-mixed-models at r-project.org
> *Subject:* Re: [R-sig-ME] Predicted density function from mixed model?
>
>   Hi Hugh,
>
>  I *think* this is probably what youre referring to:
>
>  # set up a new data frame with
> nd <- data.frame(Genotype=testDataset$Genotype, Gender=testDataset$Gender)
>  # predict new values omitting the random effect (pop. level)
>  nd$pred <- predict(model, newdata=nd, level=0)
>
> plot(density(subset(nd, Genotype=="mutant")$pred))
> plot(density(subset(nd, Genotype=="control")$pred))
>
>  However it doesn't mean much because the predicted values only take one
> of 4 numbers depending on male/female/ & control/mutant:
>
> unique(nd$pred)
> #0.6441061 0.5628385 0.6422222 0.7777778
>
>
> plot(density(subset(nd, Genotype=="mutant")$pred))
>
>
> stripchart(subset(nd, Genotype=="mutant")$pred, method="jitter", add=TRUE)
>
>
> plot(density(subset(nd, Genotype=="control")$pred))
>
>
> stripchart(subset(nd, Genotype=="control")$pred, method="jitter", add=TRUE)
>
>
> Did you perhaps mean the residuals? These are accessed via the residuals(model) command which you can then run density on.
>
>
>
>
>
>
>
>
>
> Incidentally, for your model, Weight ~ Genotype + Gender + Genotype*Gender is identical to Weight ~ Genotype*Gender because the * includes both main effects and interaction. If you want to specify just an interaction use : (e.g. Genotype:Gender)
>
>
>
>
>
>
>
>
>
>
> HTH
>
> Alan
>
>
>  --------------------------------------------------
> Email: aghaynes at gmail.com
> Mobile: +41763389128
> Skype: aghaynes
>
>
> On 18 February 2013 20:35, Hugh Morgan <h.morgan at har.mrc.ac.uk> wrote:
>
>> 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
>>
>>
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>>
>
> 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
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

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