library(asreml) set.seed(1234) k=10 #no. of cluster n=100 #no. of obs. per cluster N=n*k #total no. of obs. Z <- diag(k)%x% rep(1,n) indx <-c(sort(rep(1:k,n))) clust <- factor(indx) #######SIMULATED PARAMETERS mu=1; b=1 #Fixed effects var.a_mean=2 #Variance of random effects in the mean model mu.v=1; b.v=1 #Fixed effects in the variance model var.a_var=1 #Variance of random effects in the variance model ############################# ###SIMULATE VALUES x <- rbinom(N,1,0.5) a <- rnorm(k,0,sqrt(var.a_mean)) a.v <- rnorm(k,0,sqrt(var.a_var)) eta.v <- mu.v + b.v*x + Z%*%a.v e <- rnorm(N,0,sqrt(exp(eta.v))) y <- mu + b*x + Z%*%a + e #### ############################# ###Analysis #DHGLM calculations mean.w <- c(rep(1,N)) res.var = 0 conv.crit = 0.00001 max.iter = 20 i.iter = 0 while (i.iterconv.crit) { i.iter=i.iter+1 mean.model <- asreml(y~x,random = ~clust, weights=mean.w, control = asreml.control(maxiter=100)) asreml.hv <- mean.model$hat res.var <- mean.model$sigma2 #Residual variance res<-resid(mean.model, type = "response") hv.true <- asreml.hv*mean.w/res.var #Hat values # These are checks for extreme cases tol.val=0.0001 for (i in 1:N) { if (abs(res[i])(1-tol.val)) hv.true[i]=(1-tol.val) } y_d<-(res^2)/(1-hv.true) #Response for variance model var.w <- (1-hv.true)/2 #Weights for variance model var.model <- asreml(y_d~x,random = ~clust, weights=var.w, family=asreml.Gamma(link=log), control = asreml.control(maxiter=100)) mean.w <- 1/fitted(var.model) #Update weights for mean model } ############################ ####Variance component estimates ####The residual variance for the mean model should be 1.0 at convergence print("Results from the mean model") print(summary(mean.model)$varcomp) print("Results from the variance model") print(summary(var.model)$varcomp) ########################################### # 2013 February 1, some WinBUGS code added for those interested #If you have WinBUGS installed you can compare the results to those from Gibbs Sampler #WARNING: It can take a long time to run it and it is difficult to get mixing between chains #If you want to do it anyway set compare2winbugs = TRUE and specify suitable paths compare2winbugs = FALSE model.file.path = "C:/BUGS/" working.directory.path = "C:/BUGS/" winbugs.program.path = "C:/Program Files/WinBUGS14/" #No need to change the code from here. if (compare2winbugs) { library("R2WinBUGS") model.file.name = paste(model.file.path,"dhglm_test.bugs", sep="") ##### #Specifying the BUGS model model.bugs="model{ for(i in 1:n){ mu[i] <- b0 + b1*x[i] + u[id[i]] logs2[i] <- c0 + c1*x[i] + v[id[i]] tau.y[i] <- 1/exp(logs2[i]) y[i]~dnorm(mu[i], tau.y[i]) } for (j in 1:k) { u[j]~dnorm(0,tau.u) v[j]~dnorm(0,tau.v) } b0~dnorm(0,0.001) b1~dnorm(0,0.001) c0~dnorm(0,0.1) c1~dnorm(0,0.1) tau.u <- 1/sig2.u tau.v <- 1/sig2.v sig2.u~dunif(0.01,10) sig2.v~dunif(0.01,10) }" write.table(model.bugs, file=model.file.name, quote = FALSE, row.names = FALSE, col.names = FALSE) ### Y <- list(y=as.numeric(y), n=N, k=k, id=indx, x=x) ################### model1 <- bugs(data=Y, inits=NULL, n.burnin=18000, n.iter=20000, parameters.to.save=c("b0", "c0", "b1", "c1", "sig2.u", "sig2.v"), model.file=model.file.name, working.directory=working.directory.path, bugs.directory=winbugs.program.path, debug=T) par(mfrow=c(3,2)) plot(density(model1$sims.list$b0)) plot(density(model1$sims.list$c0)) plot(density(model1$sims.list$b1)) plot(density(model1$sims.list$c1)) plot(density(model1$sims.list$sig2.u)) plot(density(model1$sims.list$sig2.v)) }