[R-sig-ME] Simulate type I error

[Phillip Alday]

Since I haven't seen an answer go by, I'll try to give you a quick
answer and hope that it moves you along. (Thanks for including an MWE,
but I don't have time at the moment to check it. :/ )

When looking at type I error -- false positives -- you're looking for
significant results when the null hypothesis holds. So the thing you
need to simulate is the null hypothesis because you then you know it
holds! So the usual way to do this (using the likelihood-ratio test)
would be:


1. simulate data for the null hypothesis
2. fit the null model m0
3. fit the alternative model m1
4. compute the LRT for m0 vs m1 and store the p-value (not just the
significance!)
5. goto (1) and repeat many times
6. take the p-values you collected and examine them. The simplest thing
in in R would be mean(ps < 0.05) to get your type I error rate. But you
can also look at the standard error or confidence interval on that
estimate (it's statistics all the way down), or plot the distribution of
p-values or .... just a lot of things.
7. If you want to do this for lots of different simulation settings,
then you can do a plot at the end with the simulation parameters on one
axis and the error-rate + associated CI to see how much your error rate
would change under different ground truths.

By the way, the simr package does something for different sample sizes
to create a "power curve" of how the power changes. That could also be a
nice place to look for inspiration on visualization.

On 27/1/22 12:00 pm, Chao Liu wrote:
> Dear R-sig-mixed-models community,
> 
> I would like to simulate type I error for a random-effects model I
> generated.
> 
> The statistic of interest is standard deviations of the random intercept
> and random slope. Specifically, for random intercept, H_{0}: lambda_{0} =2
> and H_{1}: lambda_{0} not equal to 2; for random slope, H_{0}: lambda_{1}
> =1 and H_{1}: lambda_{1} not equal to 1. I assume the test would be a
> likelihood ratio test but please correct me if I am wrong. How do I assess
> type I error for the random-effects model I specified below:
> 
> set.seed(323)
> #The following code is to specify the structure and parameters of the
> random-effects model
> dtfunc = function(nsub){
>   time = 0:9
>   rt = c()
>   time.all = rep(time, nsub)
>   subid.all = as.factor(rep(1:nsub, each = length(time)))
> 
>   # Step 1:  Specify the lambdas.
>   G = matrix(c(2^2, 0, 0, 1^2), nrow = 2)
>   int.mean = 251
>   slope.mean = 10
>   sub.ints.slopes = mvrnorm(nsub, c(int.mean, slope.mean), G)
>   sub.ints = sub.ints.slopes[,1]
>   time.slopes = sub.ints.slopes[,2]
> 
>   # Step 2:  Use the intercepts and slopes to generate RT data
>   sigma = 30
>   for (i in 1:nsub){
>     rt.vec = sub.ints[i] + time.slopes[i]*time + rnorm(length(time), sd =
> sigma)
>     rt = c(rt, rt.vec)
>   }
> 
>   dat = data.frame(rt, time.all, subid.all)
>   return(dat)
> }
> 
> #Here I run one random-effects model
> set.seed(10)
> dat = dtfunc(16)
> lmer(rt~time.all + (1+time.all |subid.all), dat)
> 
> Assuming the test for significance is likelihood ratio test and so in the
> end, I want to see if I run the test 1000 times, what is the probability of
> rejecting the null hypothesis when it is TRUE. Also, how do I plot the
> behavior of type I error if I change the values of standard deviations?
> 
> Any help is appreciated!
> 
> Best,
> 
> Chao
> 
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>