r - plot mixed effects model in ggplot -

i new mixed effect models , need please. have plotted below graph in ggplot:

ggplot(tempef,aes(trtyear,co2effect,group=myc,col=myc)) +    facet_grid(~n) +   geom_smooth(method="lm",se=t,size=1) +   geom_point(alpha = 0.3) +    geom_hline(yintercept=0, linetype="dashed") +   theme_bw() 

however, represent mixed effects model instead of lmin geom_smooth, can include siteas random effect.

the model following:

library(lme4) tempef$trtyear <- as.numeric(tempef$trtyear) mod <- lmer(r ~ myc * n * trtyear + (1|site), data=tempef) 

i have included trtyear(year of treatment) because interested in patterns of effect, may increase or decrease on time groups.

next best attempt far extract plotting variables out of model, extracted values trtyear= 5, 10 , 15.

library(effects) ef <- effect("myc:n:trtyear", mod) x <- as.data.frame(ef) > x    myc     n trtyear        fit         se       lower     upper 1    nlow       5 0.04100963 0.04049789 -0.03854476 0.1205640 2  ecm  nlow       5 0.41727928 0.07342289  0.27304676 0.5615118 3   nhigh       5 0.20562700 0.04060572  0.12586080 0.2853932 4  ecm nhigh       5 0.24754017 0.27647151 -0.29556267 0.7906430 5    nlow      10 0.08913042 0.03751783  0.01543008 0.1628307 6  ecm  nlow      10 0.42211957 0.15631679  0.11504963 0.7291895 7   nhigh      10 0.30411129 0.03615213  0.23309376 0.3751288 8  ecm nhigh      10 0.29540744 0.76966410 -1.21652689 1.8073418 9    nlow      15 0.13725120 0.06325159  0.01299927 0.2615031 10 ecm  nlow      15 0.42695986 0.27301163 -0.10934636 0.9632661 11  nhigh      15 0.40259559 0.05990085  0.28492587 0.5202653 12 ecm nhigh      15 0.34327471 1.29676632 -2.20410343 2.8906529 

suggestions different approach represent analysis welcome. thought question better suited stackoverflow because it’s technicalities in r rather statistics behind. thanks

you can represent model variety of different ways. easiest plot data various parameters using different plotting tools (color, shape, line type, facet), did example except random effect site. model residuals can plotted communicate results. @mrflick commented, depends on want communicate. if want add confidence/prediction bands around estimates, you'll have dig deeper , consider bigger statistical issues (example1, example2).

here's example taking yours bit further.
also, in comment said didn't provide reproducible example because data not belong you. doesn't mean can't provide example out of made data. please consider future posts can faster answers.

#make data: tempef <- data.frame(   n = rep(c("nlow", "nhigh"), each=300),   myc = rep(c("am", "ecm"), each=150, times=2),   trtyear = runif(600, 1, 15),   site = rep(c("a","b","c","d","e"), each=10, times=12)   #5 sites   )  # make response data tempef$r <- 2*tempef$trtyear +                                -8*as.numeric(tempef$myc=="ecm") +             4*as.numeric(tempef$n=="nlow") +             0.1*tempef$trtyear * as.numeric(tempef$n=="nlow") +             0.2*tempef$trtyear*as.numeric(tempef$myc=="ecm") +            -11*as.numeric(tempef$myc=="ecm")*as.numeric(tempef$n=="nlow")+              0.5*tempef$trtyear*as.numeric(tempef$myc=="ecm")*as.numeric(tempef$n=="nlow")+             as.numeric(tempef$site) +  #random intercepts; intercepts increase 1            tempef$trtyear/10*rnorm(600, mean=0, sd=2)    #add noise  library(lme4) model <- lmer(r ~ myc * n * trtyear + (1|site), data=tempef) tempef$fit <- predict(model)   #add model fits dataframe 

incidentally, model fit data compared coefficients above:

model  #linear mixed model fit reml ['lmermod'] #formula: r ~ myc * n * trtyear + (1 | site) #   data: tempef #reml criterion @ convergence: 2461.705 #random effects: # groups   name        std.dev. # site     (intercept) 1.684    # residual             1.825    #number of obs: 600, groups:  site, 5 #fixed effects: #         (intercept)                mycecm                 nnlow                #             3.03411              -7.92755               4.30380                #             trtyear          mycecm:nnlow        mycecm:trtyear   #             1.98889             -11.64218               0.18589   #       nnlow:trtyear  mycecm:nnlow:trtyear   #             0.07781               0.60224       

adapting example show model outputs overlaid on data

   library(ggplot2)     ggplot(tempef,aes(trtyear, r, group=interaction(site, myc), col=site, shape=myc )) +        facet_grid(~n) +       geom_line(aes(y=fit, lty=myc), size=0.8) +       geom_point(alpha = 0.3) +        geom_hline(yintercept=0, linetype="dashed") +       theme_bw() 

notice did change colour myc site, , linetype myc.

i hope example gives ideas how visualize mixed effects model.