[R-sig-ME] Fw: Fw: simulating a linear random intercept model with exogeneity assumption
Thierry Onkelinx
thierry.onkelinx at inbo.be
Wed Jul 12 16:01:59 CEST 2017
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 <- 10
nindiv <- 20
b0 <- 2
b <- 0.2
u <- 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 op 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 op 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 op r-project.org>:
>
>
>
>
> ----- Forwarded Message -----From: mervat mohamed <
> mervat_moh2006 op yahoo.com>To: mailman-owner op r-project.org <
> mailman-owner op 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 op r-project. org
> <r-sig-mixed-models op 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 op r-project. org
> <r-sig-mixed-models op 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 op r-project. org <R-sig-mixed-models op r-project.org>
> mailing list
> https://stat.ethz.ch/mailman/ listinfo/r-sig-mixed-models
> <https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models>
>
> ______________________________ _________________
> R-sig-mixed-models op r-project. org <R-sig-mixed-models op r-project.org>
> mailing list
> https://stat.ethz.ch/mailman/ listinfo/r-sig-mixed-models
> <https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models>
>
>
> ______________________________ _________________
> R-sig-mixed-models op r-project. org <R-sig-mixed-models op r-project.org>
> mailing list
> https://stat.ethz.ch/mailman/ listinfo/r-sig-mixed-models
> <https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models>
>
>
>
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list