R help - Aug 2023 - simulating future observations from heteroscedastic fits

Hello, All:


	  I want to simulate future observations from fits to heteroscedastic 
data. A simple example is as follows:


(DF3_2 <- data.frame(y=c(1:3, 10*(1:3)),
                      gp=factor(rep(1:2, e=3))))


# I want to fit 4 models
# and simulate future observations from all 4:


fit11 <- lm(y~1, DF3_2)
fit21 <- lm(y~gp, DF3_2)

library(nlme)
(fit12 <- lme(y~1, data=DF3_2, random=~1|gp))
(fit22 <- lme(y~gp, data=DF3_2, random=~1|gp))
library(lme4)
(fit12r <- lmer(y~1+(1|gp), data=DF3_2, REML=FALSE))
(fit22r <- lmer(y~gp+(1|gp), data=DF3_2, REML=FALSE))

# I can simulate what I want for fit11 and fit21
# as follows:

simPred <- function(object, nSims=2){
   pred <- predict(object, DF3_2[6,], se.fit=TRUE,
                   interval='prediction')
   with(pred, fit[1, 'fit'] +
          se.fit*rt(nSims, df))
}
simPred(fit11)
simPred(fit21)


# How can I do the same with either fit12 and fit22
# or fit12r and fit22r?


	  Thanks,
	  Spencer Graves

On 8/31/23 1:27 AM, Spencer Graves wrote:
> Hello, All:
> 
> 
>  ????? I want to simulate future observations from fits to 
> heteroscedastic data. A simple example is as follows:
> 
> 
> (DF3_2 <- data.frame(y=c(1:3, 10*(1:3)),
>  ???????????????????? gp=factor(rep(1:2, e=3))))
> 
> 
> # I want to fit 4 models
> # and simulate future observations from all 4:
> 
> 
> fit11 <- lm(y~1, DF3_2)
> fit21 <- lm(y~gp, DF3_2)
> 
> library(nlme)
> (fit12 <- lme(y~1, data=DF3_2, random=~1|gp))
> (fit22 <- lme(y~gp, data=DF3_2, random=~1|gp))
> library(lme4)
> (fit12r <- lmer(y~1+(1|gp), data=DF3_2, REML=FALSE))
> (fit22r <- lmer(y~gp+(1|gp), data=DF3_2, REML=FALSE))
> 
> # I can simulate what I want for fit11 and fit21
> # as follows:
> 
> simPred <- function(object, nSims=2){
>  ? pred <- predict(object, DF3_2[6,], se.fit=TRUE,
>  ????????????????? interval='prediction')
>  ? with(pred, fit[1, 'fit'] +
>  ???????? se.fit*rt(nSims, df))
> }
> simPred(fit11)
> simPred(fit21)
> 
> 
> # How can I do the same with either fit12 and fit22
> # or fit12r and fit22r?
> 
I think I found it:


simPred4 <- function(object, nSims=2){
# class(object) = lmeeMod
   sim <- simulate(object, nsim=nSims,
           newdata=DF3_2[6,])
   sim
}

simPred4(fit12r)
simPred4(fit22r)
> 
>  ????? Thanks,
>  ????? Spencer Graves
> 
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.