如何用R语言写随机系数模型程序

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

单变量随机系数模型怎么用R写程序？
截距和斜率均受随机效应影响的单变量随机系数回归模型. 假定有 n 个个体被观测, 第 i 个个体被观测 mi次, 第 j 次观测结果为
yij= μi+ xijβi+ εij， i = 1, 2, · · · , n; j = 1, 2, · · · , mi. 其中 μi～ (μ, σ2μ), βi～ (β, σ2β), εij～ (0, σ2)

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：随机系数 如何用 R语言 随机效应 回归模型

1. set.seed(123)
   
2. 
   
3. i <- 1:100 # number of group
   
4. j <- rpois(n=length(i), lambda=100)+5  # use your m_i sequence
   
5. 
   
6. # mean and std of stochastic parameter
   
7. miu_mean <- rnorm(n=length(i), mean = 1, sd=2)
   
8. miu_sigma <- rpois(n=length(i), lambda=2)+1
   
9. beta_mean <- rnorm(n=length(i), mean =5, sd = 10)
   
10. beta_sigma <- rpois(n=length(i), lambda=4)+1
    
11. 
    
12. # each group have different parameter
    
13. miu_i <- vector()
    
14. beta_i <- vector()
    
15. for (index in 1:length(i)){
    
16.     temp_miu <- rnorm(n=1, mean = miu_mean[index], sd = miu_sigma[index])
    
17.     miu_i <- c(miu_i, temp_miu)
    
18.     temp_beta <- rnorm(n=1, mean = beta_mean[index], sd = beta_sigma[index])
    
19.     beta_i <- c(beta_i, temp_beta)
    
20. }
    
21. 
    
22. 
    
23. x_ij <- rnorm(n=sum(j), mean = 3)
    
24. epsilon_ij <- rnorm(n=sum(j))
    
25. miu_i <- rep(miu_i, j)
    
26. beta_i <- rep(beta_i, j)
    
27. 
    
28. y_ij <- miu_i + x_ij*beta_i + epsilon_ij

复制代码

deem 发表于 2017-8-10 18:30

非常非常感谢！我尝试一下

deem 发表于 2017-8-10 18:30

再请问一下，模型做完之后，怎样画出随机效应的 BLUE ?  以及怎样得出μi～ (μ, σ2μ), βi～ (β, σ2β), εij～ (0, σ2)这三个正态分布的具体的参数来？拜托了

wliusdau 发表于 2017-8-10 19:46
再请问一下，模型做完之后，怎样画出随机效应的 BLUE ?  以及怎样得出μi～ (μ, σ2μ), βi～ (β, σ2 ...

什么意思，我只是按照模型simulate了一个数据。如果还要估计参数的话：

1. require(lme4)
   
2. sim_dat <- data.frame(id=rep(i,j), y = y_ij, x = x_ij)
   
3. lm <- lmer(y ~ x + (1+x|id), data = sim_dat)
   
4. summary(lm)

复制代码

deem 发表于 2017-8-10 23:21
什么意思，我只是按照模型simulate了一个数据。如果还要估计参数的话：

非常感谢，解决了困扰我很久的问题。假定有 n 个个体被观测, 第 i 个个体被观测 mi次, 第 j 次观测结果为
yij= μi+ xijβi+ εij， i = 1, 2, · · · , n; j = 1, 2, · · · , mi. 其中 μi～ (μ, σ2μ), βi～ (β, σ2β), εij～ (0, σ2)
依然以这个为例，比如我按照这个程序做出一个随机系数模型来,固定效应：y=458.67*x-814.51，随机效应： μi～ (μ, 54456), βi～ (β, 15336), εij～ (0, 188)，假定他们都是正态分布，而且我的结果里面 i = 1, 2, · · · , 18; j = 1, 2, · · · , 21. 那么我想请问的是：您上面simulate了一个数据，我想按照我做出来的这个模型随便拿出一个x值然后来simulate y的值，simulate 1000次以上，这样plot（y）应该就可以看出预测值y的分布范围，但是μ和β的均值也没给出，我也不知道这该怎么设定，这样的话程序该怎么写，所以还请您多多指导，不胜感激！

wliusdau 发表于 2017-8-11 14:31
非常感谢，解决了困扰我很久的问题。假定有 n 个个体被观测, 第 i 个个体被观测 mi次, 第 j 次观测结果为 ...

如果告诉你y=458.67*x-814.51，那就是说β=458.67，μ=-814.51

deem 发表于 2017-8-11 14:38
如果告诉你y=458.67*x-814.51，那就是说β=458.67，μ=-814.51

惭愧，这次懂了，非常感谢！

deem 发表于 2017-8-11 14:38
如果告诉你y=458.67*x-814.51，那就是说β=458.67，μ=-814.51

大哥，再麻烦问一下就像上面我做出来的这个模型，我能不能用什么语句快速随机生成y1，y2，……y1000个预测值呢？