Hi,

I see.  I used subset(nd, Genotype=='mutant') to get a readable group and that is the same [1].  I think this is because they all have the same batch.  You do still appear to have all values that share a Genotype, Gender and batch.  I think this is exactly what I need, and I shall just add the individual variance:

nd$randomEffect <- rnorm(nrow(nd), mean = 0, sd = 0.1319891)
# perhaps sd should be 0.09476296
nd$predictedDistribution <- (nd$pred + nd$randomEffect)

Thank you for all your help,

Hugh

________________________________
From: Alan Haynes [aghaynes@gmail.com]
Sent: Tuesday, February 19, 2013 4:25 PM
To: Hugh Morgan
Cc: r-sig-mixed-models@r-project.org
Subject: Re: [R-sig-ME] Fwd: Predicted density function from mixed model?

Hi again Hugh,

I just ran the code with level=1 and level=0 and got two different results. level=1 arent massively different, but theyre not the same at all. level=0:1 does both population and batch level predictions, thus 2 columns.

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")

nd <- data.frame(Genotype=testDataset$Genotype, Gender=testDataset$Gender, Batch=testDataset$Batch)
# predict new values omitting the random effect (pop. level)
nd$pred <- predict(model, newdata=nd, level=0)
head(nd$pred)
[1] 0.6441061 0.6441061 0.5628385 0.6441061 0.6441061 0.5628385
nd$pred2 <- predict(model, newdata=nd, level=1)
head(nd$pred2)
[1] 0.5845188 0.5845188 0.5032512 0.5845188 0.5845188 0.5032512


HTH

Alan

--------------------------------------------------
Email: aghaynes@gmail.com<mailto:aghaynes@gmail.com>
Mobile: +41763389128
Skype: aghaynes


On 19 February 2013 15:12, Hugh Morgan <h.morgan@har.mrc.ac.uk<mailto:h.morgan@har.mrc.ac.uk>> wrote:
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@r-project.org<mailto:r-sig-mixed-models-bounces@r-project.org> [r-sig-mixed-models-bounces@r-project.org<mailto:r-sig-mixed-models-bounces@r-project.org>] on behalf of Alan Haynes [aghaynes@gmail.com<mailto:aghaynes@gmail.com>]
Sent: Tuesday, February 19, 2013 9:39 AM
To: r-sig-mixed-models@r-project.org<mailto:r-sig-mixed-models@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@gmail.com<mailto:aghaynes@gmail.com>
Mobile: +41763389128<tel:%2B41763389128>
Skype: aghaynes


On 19 February 2013 09:38, Hugh Morgan <h.morgan@har.mrc.ac.uk<mailto:h.morgan@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@gmail.com<mailto:aghaynes@gmail.com>]
> *Sent:* Tuesday, February 19, 2013 7:27 AM
> *To:* Hugh Morgan
> *Cc:* r-sig-mixed-models@r-project.org<mailto:r-sig-mixed-models@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@gmail.com<mailto:aghaynes@gmail.com>
> Mobile: +41763389128<tel:%2B41763389128>
> Skype: aghaynes
>
>
> On 18 February 2013 20:35, Hugh Morgan <h.morgan@har.mrc.ac.uk<mailto:h.morgan@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@r-project.org<mailto:R-sig-mixed-models@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@r-project.org<mailto:R-sig-mixed-models@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



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]]

