simulating a linear random intercept model with exogeneity assumption
dear ir. Thierry Onkelinx thank you very much ir. Thierry Onkelinx?i appreciate your help mervat
On Wednesday, July 12, 2017, 4:02:03 PM GMT+2, Thierry Onkelinx <thierry.onkelinx at inbo.be> wrote:
Dear Mervat,
Here is how I would write it. No need for the multivariate rnorm since you set all covariances are zero, making them independent. All rnorm() are by default independent so need for a special command to make them independent.?
library(dplyr)grp <- 10nindiv <- 20b0 <- 2b <- 0.2u <- rnorm(grp, mean = 0, sd = 0.05)e <- rnorm(grp * nindiv, mean = 0, sd = 1)expand.grid(? group = seq_len(grp),? individual = seq_len(nindiv)) %>%? mutate(? ? X = rnorm(n(), mean = 0, sd = 1),? ? id = interaction(group, individual),? ? mu = b0 + b * X + u[group],? ? y = mu + e? )
Best regards,
ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and Forest
team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
Kliniekstraat 25
1070 Anderlecht
Belgium
To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data. ~ John Tukey
2017-07-12 10:51 GMT+02:00 mervat mohamed <mervat_moh2006 at yahoo.com>:
dear ir. Thierry Onkelinx
thank you for your advise , my program code is as follows:
## The code of the program
# yj[i] = b0j+ b1*xj[i]+ej[i]
# b0j = b0 + u0j, u0j ~ N(0,1)??????????????random intercept
# b1 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ????? ? ? ? ? ?????? fixed slope
# -- load the used package
library(lme4)
set.seed(1)
????n=200000????? ?? ? ? ? ? ? ? ? ? ? ?? ? ##??the total size of the simulated population
????grp=100??????? ? ? ? ? ? ? ? ? ? ? ? ? ?? ##??the number of groups of the simulated population
????nindiv=2000??????????????? ???????????? ##??the number of observations in each group
????b0=2?????????? ? ? ? ? ? ? ? ? ? ? ? ??? ? ##??intercept
????b1=0.2?????? ? ? ? ? ? ? ? ? ?? ? ? ?? ?? ##??slope of X
??? u = rnorm(grp,0,0.05)
????e = rnorm(nindiv*grp,1)
????M.x=rnorm(grp,3,0.05)
????sigma<-diag(grp)????? ? ? ? ? ???? ## identity matrix
????require(MASS)
????x=mvrnorm(nindiv,M.x, sigma)
????X=as.vector(x)
????data.set = data.frame(group.id = rep(1:grp, each=nindiv),
??????????????????u = rep(u, each=nindiv),
??????????????????M.x = rep(M.x, each=nindiv),
??????????????????indiv.id = 1:(nindiv*grp),
??????????????????X = X,
??????????????????e = e)
??? data.set$Y = b0 + b1*data.set$X + b2*data.set$mean.X+ data.set$u + data.set$e
????indiv<- data.frame(data.set$group.id, data.set$X,data.set$Y)????? ??? names(indiv)<-c("group.id","X" ,"Y")
my problem is i do not know how can i write the command of the independence between X and the error terms of the two levels of the model (e and u ) in my program code.
and i appreciate your help
thanks mervat?
On Wednesday, July 12, 2017, 9:46:10 AM GMT+2, Thierry Onkelinx <thierry.onkelinx at inbo.be> wrote:
Dear Mervat, This the r-sig-mixed-models mailing list. Not the write-code-for-me mailing list. Please show us what you have tried. This can easily be solved using expand.grid() and rnorm(). Best regards, ir. Thierry Onkelinx Instituut voor natuur- en bosonderzoek / Research Institute for Nature and Forest team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance Kliniekstraat 25 1070 Anderlecht Belgium To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of. ~ Sir Ronald Aylmer Fisher The plural of anecdote is not data. ~ Roger Brinner The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data. ~ John Tukey 2017-07-11 13:28 GMT+02:00 mervat mohamed via R-sig-mixed-models <r-sig-mixed-models at r-project. org>: ----- Forwarded Message -----From: mervat mohamed <mervat_moh2006 at yahoo.com>To: mailman-owner at r-project.org <mailman-owner at r-project.org> Sent: Sunday, July 2, 2017, 5:39:42 AM GMT+2Subject: Fw: [R-sig-ME] Fw: simulating a linear random intercept model with exogeneity assumption ? ? ?On Wednesday, June 21, 2017 8:27 PM, mervat mohamed via R-sig-mixed-models <r-sig-mixed-models at r-project. org> wrote: ? Show original message? ? On Monday, June 19, 2017 1:36 PM, mervat mohamed via R-sig-mixed-models <r-sig-mixed-models at r-project. org> wrote: ?Hi r-sig group I want to know how to write the code of simulating a linear random intercept model with following assumptions using R programing: the model: yij=b.+b1 xij + u.j + eij, ?where j refer to the group number and i refer to the observation number in the j group model assumptions: ? 1- xij ~ N (3,1.5)? 2- u.j ~ N (0,1)? 3 - eij ~ N (0,1)? 4 - cov (xij , u.j)=0? 5 - cov (xij , eij)=0 ? 6- cov (u.j , eij)= 0 thnx for help mervat ??? [[alternative HTML version deleted]] ______________________________ _________________ R-sig-mixed-models at r-project. org mailing list https://stat.ethz.ch/mailman/ listinfo/r-sig-mixed-models ? ______________________________ _________________ R-sig-mixed-models at r-project. org mailing list https://stat.ethz.ch/mailman/ listinfo/r-sig-mixed-models ? ? ______________________________ _________________ R-sig-mixed-models at r-project. org mailing list https://stat.ethz.ch/mailman/ listinfo/r-sig-mixed-models