modèle d'effets mixtes de tracé dans ggplot

Je suis nouveau à effets mixtes modèles et j'ai besoin de votre aide s'il vous plaît.
J'ai tracé le graphique ci-dessous dans 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()

modèle d'effets mixtes de tracé dans ggplot

Cependant, je tiens à représenter un modèle à effets mixtes au lieu de lmdans geom_smoothafin que je puisse inclure SITEcomme un effet aléatoire.

Le modèle serait le suivant:

library(lme4)
tempEf$TRTYEAR <- as.numeric(tempEf$TRTYEAR)
mod <- lmer(r ~ Myc * N * TRTYEAR + (1|SITE), data=tempEf)

J'ai inclus TRTYEAR(année de traitement) parce que je suis également intéressé par les modèles de l'effet, qui peut augmenter ou diminuer avec le temps pour certains groupes.

Suivant est ma meilleure tentative à ce jour pour extraire le traçage des variables de sortie du modèle, mais seulement extrait les valeurs de TRTYEAR= 5, 10 et 15.

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

Suggestions pour une approche complètement différente pour représenter cette analyse sont les bienvenus. Je pensais que cette question est mieux adapté pour stackoverflow parce que c'est sur les aspects techniques dans la R plutôt que les statistiques derrière. Grâce

source d'informationauteur fede_luppi

  1. 18

    Vous pouvez représenter votre modèle d'une variété de façons différentes. Le plus simple est de tracer des données par les différents paramètres à l'aide de différents outils de dessins (couleur, forme, type de ligne, facette), qui est ce que vous avez fait avec votre exemple à l'exception de l'effet aléatoire site. Les valeurs résiduelles du modèle peut également être tracés pour communiquer les résultats. Comme @MrFlick commenté, ça dépend de ce que vous voulez communiquer. Si vous souhaitez ajouter de la confiance/la prédiction des bandes autour de vos estimations, vous aurez à creuser plus profondément et d'envisager de grandes questions statistiques (exemple1exemple2).

    Voici un exemple de prendre la vôtre juste un peu plus loin.

    Aussi, dans votre commentaire vous avez dit que vous n'avez pas fourni un exemple reproductible, car les données ne vous appartiennent pas. Cela ne signifie pas que vous ne pouvez pas fournir un exemple de composé de données. Veuillez prendre en compte que pour les futurs postes de sorte que vous pouvez obtenir plus rapidement des réponses.

    #Make up 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 up some 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 will increase by 1
           tempEf$TRTYEAR/10*rnorm(600, mean=0, sd=2)    #Add some noise

library(lme4)
model <- lmer(r ~ Myc * N * TRTYEAR + (1|site), data=tempEf)
tempEf$fit <- predict(model)   #Add model fits to dataframe

    D'ailleurs, l'ajustement du modèle aux données et par rapport aux coefficients ci-dessus:

    model

#Linear mixed model fit by REML ['lmerMod']
#Formula: r ~ Myc * N * TRTYEAR + (1 | site)
#   Data: tempEf
#REML criterion at 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

    L'adaptation de votre exemple pour montrer les sorties du modèle superposé sur les données 

       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()

    Avis que je n'ai fait changer votre couleur de Myc à siteet de type de ligne à Myc.
    modèle d'effets mixtes de tracé dans ggplot

    J'espère que cet exemple donne quelques idées sur la façon de visualiser votre modèle à effets mixtes.