关于R2WinBUGS数据的读入问题

dongshengjay 发表于 2013-10-25 16:33
It's OK!

我试试
谢谢你

好像算起来时间会很长

dongshengjay 发表于 2013-10-25 16:30
library(R2WinBUGS)
library(lattice)
library(coda)

现在mass是因变量，length是自变量，ngroups是区组，那pop是什么？是mass和length的总个数吗？

忘记告诉你了，winbugs很快就计算结束了，结束后需要关掉winbugs，r软件才能进行后续的运算

xingzhaoh 发表于 2013-10-25 17:23
现在mass是因变量，length是自变量，ngroups是区组，那pop是什么？是mass和length的总个数吗？

楼主应该搞清楚多水平模型的概念，pop再次应该是水平1的变量

Thank you to noting me. But to be honest, it is a long time since I use WinBUGS, thus, I can't help you. Sorry for this.

dongshengjay 发表于 2013-10-25 20:42
楼主应该搞清楚多水平模型的概念，pop再次应该是水平1的变量

大侠，你能帮我是试试吗？怎么把pop改成水平1的变量（mass是因变量，length是自变量，ngroups是区组）
数据改成什么样子才能算出来？这几天都没做出来，如果可以的话你用附件的数据试试，告诉我一下怎么做，现在头都大了，怎么都做不出来了。

dongshengjay 发表于 2013-10-25 20:42
楼主应该搞清楚多水平模型的概念，pop再次应该是水平1的变量

出自Introduction to WinBUGS for Ecologists这本书的12章（麻烦您了，你费心看看怎么可以做出来，如果这个麻烦的话，你那有现成的例子让我用也可以，就是利用winbug算混合模型，截距、斜率、截距和斜率分别作为随机变量时怎么拟合出来，感觉他给的例子太麻烦，你那有现成的就更好了）
Chapter 12: Linear mixed-effects model
全部代码

1. n.groups <- 56                                # Number of populations
   
2. n.sample <- 10                                # Number of vipers in each pop
   
3. n <- n.groups * n.sample                 # Total number of data points
   
4. pop <- gl(n = n.groups, k = n.sample)         # Indicator for population
   
5. 
   
6. # Body length (cm)
   
7. original.length <- runif(n, 45, 70)
   
8. mn <- mean(original.length)
   
9. sd <- sd(original.length)
   
10. cat("Mean and sd used to normalise.original length:", mn, sd, "\n\n")
    
11. length <- (original.length - mn) / sd
    
12. hist(length, col = "grey")
    
13. 
    
14. Xmat <- model.matrix(~pop*length-1-length)
    
15. print(Xmat[1:21,], dig = 2)                 # Print top 21 rows
    
16. 
    
17. intercept.mean <- 230                        # mu_alpha
    
18. intercept.sd <- 20                        # sigma_alpha
    
19. slope.mean <- 60                        # mu_beta
    
20. slope.sd <- 30                                # sigma_beta
    
21. 
    
22. intercept.effects<-rnorm(n = n.groups, mean = intercept.mean, sd = intercept.sd)
    
23. slope.effects <- rnorm(n = n.groups, mean = slope.mean, sd = slope.sd)
    
24. all.effects <- c(intercept.effects, slope.effects) # Put them all together
    
25. 
    
26. lin.pred <- Xmat[,] %*% all.effects        # Value of lin.predictor
    
27. eps <- rnorm(n = n, mean = 0, sd = 30)        # residuals
    
28. mass <- lin.pred + eps                        # response = lin.pred + residual
    
29. 
    
30. hist(mass, col = "grey")                # Inspect what we’ve created
    
31. 
    
32. library("lattice")
    
33. xyplot(mass ~ length | pop)
    
34. 
    
35. 
    
36. 
    
37. ### 12.3. Analysis under a random-intercepts model
    
38. 
    
39. 
    
40. ### 12.3.1. REML analysis using R
    
41. library('lme4')
    
42. lme.fit1 <- lmer(mass ~ length + (1 | pop), REML = TRUE)
    
43. lme.fit1
    
44. 
    
45. 
    
46. ### 12.3.2. Bayesian analysis using WinBUGS
    
47. # Write model
    
48. sink("lme.model1.txt")
    
49. cat("
    
50. model {
    
51. 
    
52. # Priors
    
53. for (i in 1:ngroups){               
    
54.     alpha[i] ~ dnorm(mu.int, tau.int)        # Random intercepts
    
55. }
    
56. 
    
57. mu.int ~ dnorm(0, 0.001)                # Mean hyperparameter for random intercepts
    
58. tau.int <- 1 / (sigma.int * sigma.int)
    
59. sigma.int ~ dunif(0, 100)                # SD hyperparameter for random intercepts
    
60. 
    
61. beta ~ dnorm(0, 0.001)                        # Common slope
    
62. tau <- 1 / ( sigma * sigma)        # Residual precision
    
63. sigma ~ dunif(0, 100)                        # Residual standard deviation
    
64. 
    
65. # Likelihood
    
66. for (i in 1:n) {
    
67.     mass[i] ~ dnorm(mu[i], tau)                # The random variable
    
68.     mu[i] <- alpha[pop[i]] + beta* length[i] # Expectation
    
69. }
    
70. }
    
71. ",fill=TRUE)
    
72. sink()
    
73. 
    
74. # Bundle data
    
75. win.data <- list(mass = as.numeric(mass), pop = as.numeric(pop),
    
76. length = length, ngroups = max(as.numeric(pop)), n = n)
    
77. 
    
78. # Inits function
    
79. inits <- function(){list(alpha = rnorm(n.groups, 0, 2), beta = rnorm(1, 1, 1),
    
80. mu.int = rnorm(1, 0, 1), sigma.int = rlnorm(1), sigma = rlnorm(1))}
    
81. 
    
82. # Parameters to estimate
    
83. parameters <- c("alpha", "beta", "mu.int", "sigma.int", "sigma")
    
84. 
    
85. # MCMC settings
    
86. ni <- 2000
    
87. nb <- 500
    
88. nt <- 2
    
89. nc <- 3
    
90. 
    
91. # Start Gibbs sampling
    
92. out <- bugs(win.data, inits, parameters, "lme.model1.txt", n.thin=nt,
    
93. n.chains=nc, n.burnin=nb, n.iter=ni, debug = TRUE)
    
94. 
    
95. # Inspect results
    
96. print(out, dig = 3)

复制代码

截距及系数的随机效应模型，与单个截距效应模型相似，在指定先验后，将pop的脚标加上即可mu[i] <- alpha[pop[i]] + beta[pop[i]]* length[i] ，这本书的的p159-12.4.2的模型即可，并不复杂，你可以模仿
如果pop对应你的数据每个不同的分组。

dongshengjay 发表于 2013-10-28 16:10
截距及系数的随机效应模型，与单个截距效应模型相似，在指定先验后，将pop的脚标加上即可mu

还是没这么明白，pop和group怎么设置，您能帮用我的data1 数据试试这个例子吗？我没看懂你的回复，基础太差，希望您可以用我的data，做出个成功的例子，教教我，谢谢您

library(R2WinBUGS)
library(lattice)
library(coda)
dat=read.csv("c:/data.csv")
attach(dat)
setwd("c:/")

sink("lme.model2.txt")
cat("
model {

# Priors
for (i in 1:6){               
    alpha[i] ~ dnorm(mu.int, tau.int)        # Random intercepts
    beta[i] ~ dnorm(mu.slope, tau.slope)# Random slopes
}

mu.int ~ dnorm(0, 0.001)                # Mean hyperparameter for random intercepts
tau.int <- 1 / (sigma.int * sigma.int)
sigma.int ~ dunif(0, 100)                # SD hyperparameter for random intercepts

mu.slope ~ dnorm(0, 0.001)                # Mean hyperparameter for random slopes
tau.slope <- 1 / (sigma.slope * sigma.slope)
sigma.slope ~ dunif(0, 100)                # SD hyperparameter for slopes

tau <- 1 / ( sigma * sigma)                # Residual precision
sigma ~ dunif(0, 100)                        # Residual standard deviation

# Likelihood
for (i in 1:n) {
    mass[i] ~ dnorm(mu[i], tau)
    mu[i] <- alpha[pop[i]] + beta[pop[i]]* length[i]
}
}
",fill=TRUE)
sink()

# Bundle data
win.data <- list(mass = as.numeric(mass), pop = pop,
length = length, n =length(mass))

# Inits function
inits <- function(){ list(alpha = rnorm(6, 0, 2), beta = rnorm(6, 10, 2),
mu.int = rnorm(1, 0, 1), sigma.int = rlnorm(1), mu.slope = rnorm(1, 0, 1),
sigma.slope = rlnorm(1), sigma = rlnorm(1))}

# Parameters to estimate
parameters <- c("alpha", "beta", "mu.int", "sigma.int", "mu.slope", "sigma.slope", "sigma")

# MCMC settings
ni <- 2000
nb <- 500
nt <- 2
nc <- 3

# Start Gibbs sampling
out <- bugs(win.data, inits, parameters, "lme.model2.txt", n.thin=nt,
n.chains=nc, n.burnin=nb, n.iter=ni, debug = TRUE)

print(out, dig = 2)