Skip to content

MCMCglmm prior distributions

9 messages · Jarrod Hadfield, Boby Mathew

Original
Boby Mathew
# 
Dear MCMCglmm users,

Is it possible to use double exponential priors(Laplace) in MCMCglmm?

Thanks for the helps.

regards,
Boby
Jarrod Hadfield
# 
Hi Boby,

In short - no. I haven't tried this (or thought about it much), but  
you could treat each fixed effect as a single random effect with its  
own associated variance component. Presumably, you could then specify  
the prior for the variance component in a way that induces a prior  
t-distribution on the effect. Like the Laplace it has fatter tails  
than the Normal, but it lacks the peakiness and won't give some of the  
nice features of the LASSO.

Cheers,

Jarrod



Quoting Boby Mathew <bobyboby at gmail.com> on Thu, 16 Oct 2014 16:06:13 +0200:

  
    

      
      
    

  


  


      

    

    
3 days later

    

      
      

        


  

        

      Boby Mathew
    

    

      
      #
    

  


  

      
Dear Jarrod Hadfield,

Here I have attached a small code with the simulation code. I want to
estimate the effect of 'b' here. As you suggested I treated fixed effect as
random and gave own variance. But I am not sure this is the right way.

Could you please check whether the implementation is right?

regards,
Boby

mark=100; line=150

x=round(matrix(runif(mark*line),nrow=mark))
b=rep(0,mark)
b[8]=3; b[80]=5;b[90]=5;

noise=rnorm(line,0,sqrt(1))



Line=1:line

y = b%*%x + noise

Z=t(x)

library(MCMCglmm)

data=data.frame(Phe=t(y),animal=Line)

data$animal=as.factor(data$animal)


prior2.2 <- list(G = list(G1 = list(V = 1, n = 0.002)), R = list(V = 1, n =
0.002))

mod_mcmc=MCMCglmm(Phe~1,random=~idv(Z),pr=T,data=data,nitt=50000,thin=500,burnin=10000,prior=prior2.2)

val=colMeans (mod_mcmc$Sol)




On Thu, Oct 16, 2014 at 5:51 PM, Jarrod Hadfield <j.hadfield at ed.ac.uk>
wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Jarrod Hadfield
    

    

      
      #
    

  


  

      
Hi,

This gives the b's a common variance. There is no point giving them  
individual variances unless you want to treat them as `fixed' but  
place a t-prior rather than a normal prior on each effect.

Jarrod




Quoting Boby Mathew <bobyboby at gmail.com> on Mon, 20 Oct 2014 14:42:02 +0200:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Boby Mathew
    

    

      
      #
    

  


  

      
Dear jarrod hadfield,

How can I pace  individual variances with a t-prior for the random effects?
Could you please provide me an example?

thanks for the help.

regards,
Boby

On Mon, Oct 20, 2014 at 2:58 PM, Jarrod Hadfield <j.hadfield at ed.ac.uk>
wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Jarrod Hadfield
    

    

      
      #
    

  


  

      
Hi Boby,

If you mean a common t-prior (with estimated scale) for the random  
effects then you cannot. All that you can do is place a fixed t-prior  
on the `fixed' effects.

Cheers,

Jarrod




Quoting Boby Mathew <bobyboby at gmail.com> on Mon, 20 Oct 2014 15:52:24 +0200:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Boby Mathew
    

    

      
      #
    

  


  

      
Hello Jarrod Hadfield,

Thank you so much for your help.  I have the winbugs code for estimation of
my marker effect but the problem is winbugs is too slow and cannot handle
large datasets.

my models is "y(observation) =Z%*%beta + noise"

here 'Z is the markers matrix coded 0 and 1 and beta is marker effect  I
used MCMCglmm for this model and included the marker matrix(Z) as random
(random=~idv(Z)).

I simulated some data and MCMCglmm was giving good results when the number
of markers were less than the number of observation. Does MCMCglmm can
handle this type of models?



I have attached the winbugs code here  for the reference.


model{
for(i in 1:n){
    y[i]~dnorm(mu[i],prec)
    mu[i]<- inprod(x[i,], beta[])
        }
for (j in 1: p){
    beta[j]~dnorm(0,tau[j])
    tau[j]<-1/var[j]
    var[j]~dgamma(0.1,0.1)
    }
sd~dunif(0,10)
sigma2<-sd*sd
prec<-1/sigma2

    }

Once again thanks for the help

regards,
Boby

On Mon, Oct 20, 2014 at 4:17 PM, Jarrod Hadfield <j.hadfield at ed.ac.uk>
wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Jarrod Hadfield
    

    

      
      #
    

  


  

      
Hi,

If the beta's are drawn from a normal distribution with the same  
variance then the idv(Z) model is appropriate. I thought you would get  
an error message with your code because Z was not in the data.frame.  
If this was not the case could you send me your sessionInfo()?

Cheers,

JArrod







Quoting Boby Mathew <bobyboby at gmail.com> on Mon, 20 Oct 2014 17:09:52 +0200:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Boby Mathew
    

    

      
      #
    

  


  

      
Hello Jarrod Hadfield,

Below I have attached the sessionInfo().  So if  beta's~ N(0,sigma^2)  then
the G structure assign prior for sigma^2~Inv-Gamma(a,b). Is it right?.

Also is it possible to assign  beta's some other distributions?

Thanks for the help.

regards,
Boby



R version 3.1.0 (2014-04-10)
Platform: x86_64-pc-linux-gnu (64-bit)

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C
 [9] LC_ADDRESS=C               LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] splines   stats     graphics  grDevices utils     datasets  methods
[8] base

other attached packages:
[1] INLA_0.0-1392038736 sp_1.0-15           MCMCglmm_2.21
[4] ape_3.1-4           coda_0.16-1         lattice_0.20-29
[7] Matrix_1.1-4

loaded via a namespace (and not attached):
[1] corpcor_1.6.7 grid_3.1.0    nlme_3.1-118  tcltk_3.1.0   tensorA_0.36
[6] tools_3.1.0


On Mon, Oct 20, 2014 at 7:03 PM, Jarrod Hadfield <j.hadfield at ed.ac.uk>
wrote: