TFM_ANNEX-1_SIMULACIÓ_MODELS_MIXTOS

usePackage <- function(p) {    
    if (!is.element(p, installed.packages()[,1]))
        install.packages(p, dep = TRUE)
    require(p, character.only = TRUE)
}
usePackage("nlme")
usePackage("nlmeU")
usePackage("lme4")
usePackage("ggplot2")
usePackage("gridExtra")
usePackage("ggExtra")
usePackage("kableExtra")
usePackage("RLRsim")
usePackage("mixlm")
usePackage("magrittr")
usePackage("evaluate")

ctrl <- lmeControl(opt='optim') # Mesura de control per l'aplicació dels paquets lme per evitar errors per incorrecte optimització d'iteracions

MODEL MIXT - DATASET HOMOGENI Ho

GENERACIÓ DATASET

Per la generació dels datasets s’utilitza el generador de valors de distribució uniforme mitjançant el paquet Uniform.

Es genera un dataset simulat amb rangs mitjos en quant a individus i temps (posteriorment es valorarà la opció de fer un escombrat en tot el rang d’aquestes 2 variables):

set.seed(7982)
i3 <- 3
i10 <- 10

j4 <- 4
j10 <- 10

(av_i <- round(mean(c(3,10))))

## [1] 6

(av_j <- round(mean(c(4,10))))

## [1] 7

TA <- c(0,3,6,9,12,18,24,30,36,42) # Es genera un vector amb els temps més habituals d'estabilitat del qual es prenen els valors segons la longitud del dataset a utilitzar (es tindrà en compte en el mateix ordre que està posada encara que es podria contemplar agafant valors de manera equilibrada entre el principi i el final)
TA_av <- TA[1:av_j]
Av_DS <- data.frame(Lot = rep(paste("Lot",seq(1:av_i)),each=av_j))
Av_DS$NTemps <- as.factor(rep(paste("Temps",seq(1:av_j)),av_i))
Av_DS$Temps <- rep(TA[1:av_j],av_i)

Per la generació de datasets s’utilitza sempre la mateixa llavor essent aquesta 7982.

GENERACIÓ DATASET Ho1 (NÚVOL DE PUNTS HOMOGENI)

Es generen punts a l’atzar dins un rang reduït (a efectes d’escalat al voltant de 0 entre -1 i 1) mitjançant la funció runif per obtenir un nuvol de punts amb variància homogènia:

set.seed(7982)
Av_DS$Ho1 <- runif(n=av_i*av_j, min=-1, max=1) #Observacions constant de rang petit entre -1 i 1 per simular que no hi Ha tendències ni diferències

kable(head(Av_DS[,c(colnames(Av_DS)[1:3],"Ho1")]), booktabs = T) %>% kable_styling(latex_options =  c("striped","condensed"))

Lot	NTemps	Temps	Ho1
Lot 1	Temps 1	0	-0.1692796
Lot 1	Temps 2	3	0.2031798
Lot 1	Temps 3	6	0.6117839
Lot 1	Temps 4	9	-0.4132929
Lot 1	Temps 5	12	0.4868403
Lot 1	Temps 6	18	-0.1104966

GENERACIÓ DATASET Ho2 (HOMOGENI PERÒ AMB VARIACIÓ DEPENENT DEL LOT)

En aquesta ocasió es genera un dataset amb variància homogènia, però que depengui del lot. Per fer això es considera un rang petit de variància entre medicions de cada lot i una variació de les mitges de cada lot aleatòria entre un ordre superior de magnitud. Per fer es fa servir la funció runif combinada amb alguna de les funcions de la família apply per modificar el rang de cada lot de manera seqüencial:

set.seed(7982)
Av_DS$Ho2 <- as.vector(sapply(1:av_i, function(i){
  m_lot <- runif(1, min=-1, max=1)
  runif(n=av_j, min=m_lot-0.1, max=m_lot+0.1)
  }))

kable(head(Av_DS[,c(colnames(Av_DS)[1:3],"Ho2")]), booktabs = T) %>% kable_styling(latex_options =  c("striped","condensed"))

Lot	NTemps	Temps	Ho2
Lot 1	Temps 1	0	-0.1489616
Lot 1	Temps 2	3	-0.1081012
Lot 1	Temps 3	6	-0.2106089
Lot 1	Temps 4	9	-0.1205956
Lot 1	Temps 5	12	-0.1803292
Lot 1	Temps 6	18	-0.2069329

VISUALITZACIÓ DATASETS GENERATS Ho

Es visualitza gràficament les diferències entre els datasets generats. Per veure amb més claredat les diferències es visualitza primer mitjançant un gràfic tipus scatterplot on es veu l’error estàndard global i la diferenciació per lot, i posteriorment es visualitza mitjançant els diagrames de caixes per apreciar més fàcilment els errors estàndards en el temps:

SCATTERPLOTS Ho

es veu els datasets amb variància homogència.

labsHo_1 <- c("Dataset Ho1 - Núvol de punts","Dataset Ho2 - Homogeneitat en lot")
labsHo_2 <- c("Dataset Ho1 s/ lot","Dataset Ho2 s/ lot")

HoVars <- c("Ho1","Ho2")

Av_DS_plotHo_sc1 <- lapply(1:length(HoVars),function(graph){
  ggplot(Av_DS, aes(x=Temps, y=Av_DS[,HoVars[graph]])) + geom_point() + geom_smooth(linetype="dashed") + 
    labs(title=labsHo_1[graph], y = "Observacions")
  })

Av_DS_plotHo_sc2 <- lapply(1:length(HoVars),function(graph){
  ggplot(Av_DS, aes(x=Temps, y=Av_DS[,HoVars[graph]], color=Lot, group=Lot)) + geom_point() + geom_line() + 
    labs(title=labsHo_2[graph], y = "Observacions")
  })

ydens_DS_plotHo <- lapply(1:length(HoVars),function(graph){
  ggplot(Av_DS, aes(x=Av_DS[,HoVars[graph]])) + geom_density(color="darkblue", fill="lightblue") +
    geom_vline(aes(xintercept=mean(Av_DS[,HoVars[graph]])), color="blue", linetype="dashed", size=1)+ 
    coord_flip() + labs(x = "Observacions")
  })

grid.arrange(Av_DS_plotHo_sc1[[1]], ydens_DS_plotHo[[1]], Av_DS_plotHo_sc1[[2]], ydens_DS_plotHo[[2]],
             ncol=2, widths=c(4, 1))

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

grid.arrange(Av_DS_plotHo_sc2[[1]], ydens_DS_plotHo[[1]], Av_DS_plotHo_sc2[[2]], ydens_DS_plotHo[[2]],
             ncol=2, widths=c(4, 1))

Els gràfics de densitat ens mostren que en el cas del núvol de punts les dades tenen un repartiment bastant homogeni en tot el rang de observacions, mentre que en el cas del repartiment per lot s’observa que es centren on aleatòriament els lots hagin centrat la seva mitja. Tot i que en el gràfic sense distinció de lot es veu poques diferèncie entre els 2 conjunts, en el gràfic posterior es veu que realment sí que es diferencien en el grup que pertany cada dada.

BOXPLOTS Ho

es veu els datasets amb variància homogència.

Av_DS_plotHo_bp1 <- lapply(1:length(HoVars),function(graph){
  ggplot(Av_DS, aes(x=NTemps, y=Av_DS[,HoVars[graph]])) + 
  geom_boxplot()+ geom_dotplot(binaxis='y', stackdir='center', dotsize=1)+
  labs(title=labsHo_1[graph], y = "Observacions")
  })

grid.arrange(Av_DS_plotHo_bp1[[1]], ydens_DS_plotHo[[1]], Av_DS_plotHo_bp1[[2]], ydens_DS_plotHo[[2]],
             ncol=2, widths=c(4, 1))

## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

A part de la diferència d’escala del rang en aquests boxplot tot i haver petites diferències no es pot veure clarament que no hi hagi homogeneitat, tant en el cas del núvol de punts com el de diferenciació per lots.

MODELS LINEALS Ho

Per l’ajust dels models lineals dels datasets en aquest cas es consideren diversos models diferents per comparar entre ells. Es pretén generar un model lineal simple i un model lineal d’efectes mixtos amb la intercepció aleatòria. És important tenir en compte que els datasets generats s’han de considerar independents entre ells i per tant no són directament comparables al ajustar els models, pel que l’anàlisi que es segueix és veure les diferències entre l’ajust dels diferents models a un mateix dataset i posteriorment es compara els diferents datasets, però en termes de la bondat d’ajust dels models aplicats:

ML I/IS Ho

Es genera el model lineal simple i el model amb la variable temps per comparar-los i definir el model de base simple. (s’aprofita per convertir-ho en una funció R per repetir aquest anàlisis en els altres datasets):

Comp_I_IS_LM <- function(var,ds){
  (ML_I_s <- summary(ML_I <- lm(paste(var,"~",1,sep=" "), data = ds)))
  (ML_IS_s <- summary(ML_IS <- lm(paste(var,"~","Temps",sep=" "), data = ds)))
  (ML_I_IS_an <- anova(ML_IS,ML_I))
  
  l <- list(ML_I_s,ML_IS_s,ML_I_IS_an,ML_I,ML_IS)
  names(l) <- c("OLS I","OLS IS","Anova comp. efecte S","OLS I base","OLS IS base")
  return(l)
}

(Comp_I_IS_Ho1 <- Comp_I_IS_LM("Ho1",Av_DS))[1:3]

## $`OLS I`
## 
## Call:
## lm(formula = paste(var, "~", 1, sep = " "), data = ds)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.94876 -0.46492 -0.03899  0.50629  0.99498 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.03198    0.08944  -0.358    0.723
## 
## s: 0.5796 on 41 degrees of freedom
## 
## 
## $`OLS IS`
## 
## Call:
## lm(formula = paste(var, "~", "Temps", sep = " "), data = ds)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.93932 -0.41288 -0.05226  0.47529  0.92559 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.037407   0.148827   0.251    0.803
## Temps       -0.006746   0.011512  -0.586    0.561
## 
## s: 0.5843 on 40 degrees of freedom
## Multiple R-squared: 0.008511,
## Adjusted R-squared: -0.01628 
## F-statistic: 0.3434 on 1 and 40 DF,  p-value: 0.5612 
## 
## 
## $`Anova comp. efecte S`
## Analysis of Variance Table
## 
## Model 1: Ho1 ~ Temps
## Model 2: Ho1 ~ 1
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     40 13.658                           
## 2     41 13.775 -1  -0.11725 0.3434 0.5612

(Comp_I_IS_Ho2 <- Comp_I_IS_LM("Ho2",Av_DS))[1:3]

## $`OLS I`
## 
## Call:
## lm(formula = paste(var, "~", 1, sep = " "), data = ds)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.49918 -0.11234  0.07207  0.17875  0.27545 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -0.3836     0.0349  -10.99 8.54e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## s: 0.2262 on 41 degrees of freedom
## 
## 
## $`OLS IS`
## 
## Call:
## lm(formula = paste(var, "~", "Temps", sep = " "), data = ds)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.49084 -0.12085  0.05421  0.18772  0.26989 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.358327   0.058103  -6.167 2.75e-07 ***
## Temps       -0.002452   0.004494  -0.546    0.588    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## s: 0.2281 on 40 degrees of freedom
## Multiple R-squared: 0.007389,
## Adjusted R-squared: -0.01743 
## F-statistic: 0.2978 on 1 and 40 DF,  p-value: 0.5883 
## 
## 
## $`Anova comp. efecte S`
## Analysis of Variance Table
## 
## Model 1: Ho2 ~ Temps
## Model 2: Ho2 ~ 1
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     40 2.0817                           
## 2     41 2.0972 -1 -0.015497 0.2978 0.5883

Per començar s’observa que, com és d’esperar, afegir la variable temps al model no té sentit i no dona una significació adicional al model per explicar la resposta ja que el resultat de l’Anova corresponent entre models dona un p-valor de 0.5612 (cas Ho1) i 0.5883 (cas Ho2), pel que per regla general sempre és preferible utilitzar el model més simple possible que no afegeix tants factors redundants que poden fer augmentar la variabilitat.

En el model només amb intercepció la prova de significació de la intercepció no té molt de sentit i sol ser més important veure l’error estàndar de l’estimació del paràmetre i l’error estàndar residual. En aquest sentit si es comparen els dos datsets s’observa a primera vista que és més baix en el cas de Ho2:

• Ho1: 0.0894 (Error Estimació); 0.5796 (Error residual)
• Ho2: 0.0349 (Error Estimació); 0.2262 (Error residual)

Al tenir un patró de punts més definit a Ho2 possiblement es produeix per aquest motiu la diferència en l’ajust més simple.

ML RI Ho

Avançant a partir del model simple amb intercepció ML I es fa l’ajust de models amb efecte aleatori (es genera addicionalment la funció):

Comp_I_RI_LM <- function(var,ds) {
  fm1 <- formula(paste(var,"~",1,"+","Error(Lot)"))
  ML_I_RI_sOLS <- summary(ML_I_RI_OLS <- aov(fm1, data = ds))
  ML_I_RI_OLS
  ML_I_RI_OLS_err <- sigma(summary.lm(ML_I_RI_OLS$Within))
  
  fm2 <- formula(paste(var,"~",1))
  
  ML_I_RI_sMLE <- summary(ML_I_RI_MLE <- lme(fm2, random = ~1|Lot, data = ds, method="ML"))
  
  ML_I_RI_sREML <- summary(ML_I_RI_REML <- lme(fm2,random = ~1|Lot, data = ds, method="REML"))
  
    l <- list(ML_I_RI_OLS$Within,ML_I_RI_sMLE,ML_I_RI_sREML,ML_I_RI_OLS,ML_I_RI_MLE,ML_I_RI_REML)
  names(l) <- c("OLS I RI","MLE I RI","REML I RI","OLS I RI base","MLE I RI base","REML I RI base")
  
  return(l)
} 

(Comp_I_RI_Ho1 <- Comp_I_RI_LM("Ho1",Av_DS))[1:3]

## $`OLS I RI`
## 
## Terms:
##                 Residuals
## Sum of Squares   12.90204
## Deg. of Freedom        36
## 
## Residual standard error: 0.5986569
## 
## $`MLE I RI`
## Linear mixed-effects model fit by maximum likelihood
##  Data: ds 
##        AIC      BIC    logLik
##   78.36824 83.58125 -36.18412
## 
## Random effects:
##  Formula: ~1 | Lot
##          (Intercept)  Residual
## StdDev: 1.136638e-05 0.5726893
## 
## Fixed effects: list(fm2) 
##                   Value  Std.Error DF    t-value p-value
## (Intercept) -0.03197749 0.08943905 36 -0.3575338  0.7228
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.65668158 -0.81181376 -0.06807865  0.88405597  1.73737481 
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`REML I RI`
## Linear mixed-effects model fit by REML
##  Data: ds 
##        AIC      BIC    logLik
##   81.37085 86.51157 -37.68543
## 
## Random effects:
##  Formula: ~1 | Lot
##          (Intercept)  Residual
## StdDev: 1.385267e-05 0.5796313
## 
## Fixed effects: list(fm2) 
##                   Value  Std.Error DF    t-value p-value
## (Intercept) -0.03197749 0.08943905 36 -0.3575338  0.7228
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.63684037 -0.80209109 -0.06726331  0.87346809  1.71656717 
## 
## Number of Observations: 42
## Number of Groups: 6

fmHo1 <- formula(paste("Ho1","~",1))
(eMLE_Ho1<-exactRLRT(lme(fmHo1, random = ~1|Lot, data = Av_DS, method="ML")))

## Using restricted likelihood evaluated at ML estimators.

## Refit with method="REML" for exact results.

## Warning in model.matrix.default(~m$groups[[n.levels - i + 1]] - 1,
## contrasts.arg = c("contr.treatment", : non-list contrasts argument ignored

## 
##  simulated finite sample distribution of RLRT.
##  
##  (p-value based on 10000 simulated values)
## 
## data:  
## RLRT = 0, p-value = 1

(eREML_Ho1<-exactRLRT(lme(fmHo1, random = ~1|Lot, data = Av_DS, method="REML")))

## Warning in model.matrix.default(~m$groups[[n.levels - i + 1]] - 1,
## contrasts.arg = c("contr.treatment", : non-list contrasts argument ignored

## 
##  simulated finite sample distribution of RLRT.
##  
##  (p-value based on 10000 simulated values)
## 
## data:  
## RLRT = 0, p-value = 1

(Comp_I_RI_Ho2 <- Comp_I_RI_LM("Ho2",Av_DS))[1:3]

## $`OLS I RI`
## 
## Terms:
##                 Residuals
## Sum of Squares  0.1450147
## Deg. of Freedom        36
## 
## Residual standard error: 0.06346799
## 
## $`MLE I RI`
## Linear mixed-effects model fit by maximum likelihood
##  Data: ds 
##         AIC       BIC   logLik
##   -80.06585 -74.85284 43.03293
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept)   Residual
## StdDev:   0.2142545 0.06346799
## 
## Fixed effects: list(fm2) 
##                  Value  Std.Error DF   t-value p-value
## (Intercept) -0.3835521 0.08908248 36 -4.305584   1e-04
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -1.6578389 -0.8302598  0.2248842  0.7899159  1.6254266 
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`REML I RI`
## Linear mixed-effects model fit by REML
##  Data: ds 
##         AIC       BIC   logLik
##   -77.13164 -71.99092 41.56582
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept)   Residual
## StdDev:   0.2349491 0.06346799
## 
## Fixed effects: list(fm2) 
##                  Value  Std.Error DF   t-value p-value
## (Intercept) -0.3835521 0.09641624 36 -3.978086   3e-04
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -1.6542401 -0.8289581  0.2228988  0.7935148  1.6267283 
## 
## Number of Observations: 42
## Number of Groups: 6

fmHo2 <- formula(paste("Ho2","~",1))
(eMLE_Ho2<-exactRLRT(lme(fmHo2, random = ~1|Lot, data = Av_DS, method="ML")))

## Using restricted likelihood evaluated at ML estimators.
## Refit with method="REML" for exact results.

## Warning in model.matrix.default(~m$groups[[n.levels - i + 1]] - 1,
## contrasts.arg = c("contr.treatment", : non-list contrasts argument ignored

## 
##  simulated finite sample distribution of RLRT.
##  
##  (p-value based on 10000 simulated values)
## 
## data:  
## RLRT = 81.254, p-value < 2.2e-16

(eREML_Ho2<-exactRLRT(lme(fmHo2, random = ~1|Lot, data = Av_DS, method="REML")))

## Warning in model.matrix.default(~m$groups[[n.levels - i + 1]] - 1,
## contrasts.arg = c("contr.treatment", : non-list contrasts argument ignored

## 
##  simulated finite sample distribution of RLRT.
##  
##  (p-value based on 10000 simulated values)
## 
## data:  
## RLRT = 81.33, p-value < 2.2e-16

Es podria considerar una mida mostral correcte la d’aquest model base amb 42 observacions, a part de no tenir més efectes fixos que la pròpia mitja. En aquest càlcul es comencen a veure les diferències dels datasets simulats:

• Ho1: Com s’espera no hi ha diferències molt grans entre el càlcul del model amb els 3 sistemes veient una desviació residual estàndard de (OLS Mètode ANOVA), 0.5727 (MLE) i 0.5796 (REML) ni tampoc comparat amb els errors estàndards dels models més simples sense efectes aleatoris.
• Ho2: Com en Ho1 no hi ha diferències significatives entre els 3 sistemes veient una desviació residual estàndard de (OLS Mètode ANOVA), 0.0635 (MLE) i 0.0635 (REML). A diferències de Ho1 sí que s’observa una disminució possiblement significativa de la desviació residual essent 0.2262 la del model sense efecte aleatori RI i NaN la mitja entre els errors estàndards dels 3 mètodes en els models amb l’efecte aleatori RI, a part que en aquest cas les proves de significació de l’efecte aleatori donen significatives pels dos mètodes de ML (p-valors de 0 i 0 respectivament).

Al existir un patró que modula la localització dels punts per cada lot en Ho2, encara que els lots tinguin una ubicació aleatòria, l’efecte aleatori sembla que permet ajustar aquesta característica per augmentar la certesa del model.

RESUM I DIAGNÒSTIC MODELS Ho

ML_Resum <- function(R,Key) {
  ML_table <- data.frame(Sigma = sapply(1:length(R),function(s){
    if (class(R[[s]])[1]=="aov" | class(R[[s]])[1]=="lmm"){
      round(sigma(summary.lm(R[[s]])),4)
    } else {
      round(sigma(R[[s]]),4)
    }
    })
    )
  ML_table$AIC <- sapply(1:length(R),function(s){
    round(AIC(R[[s]]),4)
    })
  ML_table$BIC <- sapply(1:length(R),function(s){
    round(BIC(R[[s]]),4)
  })
  rownames(ML_table) <- paste(Key,names(R),sep="--")
  
  R1 <- ML_table
  
  ML_table_p <- data.frame(Model = rep(rownames(ML_table),ncol(ML_table)))
  ML_table_p$Model <- factor(ML_table_p$Model, levels = rownames(ML_table))
  ML_table_p$ParQ <- unlist(ML_table)
  ML_table_p$Par <- rep(colnames(ML_table),each=nrow(ML_table))
  
  fact <- if(mean(abs(ML_table$AIC)) >= mean(abs(ML_table$Sigma))){ #Es fa un càlcul del factor més adient per graficar les 2 escales en el gràfic
    10^(round(log10(abs(mean(ML_table$AIC))/abs(mean(ML_table$Sigma)))))
  } else {
    1/(10^(round(log10(abs(mean(ML_table$Sigma))/abs(mean(ML_table$AIC))))))
  }
  
  if (fact==0){
    fact <-0.0001
  }
                
  ML_table_p$ParQ[ML_table_p$Par=="Sigma"] <- ML_table_p$ParQ[ML_table_p$Par=="Sigma"]*fact
  
  R2 <- ggplot(data=ML_table_p, aes(x=Model, y=ParQ, fill=Par)) + 
    geom_bar(stat="identity", position=position_dodge())+
    scale_fill_brewer(palette="Pastel1") + theme_minimal() +
    labs(title=Key,y = "AIC and BIC par.") +
    scale_y_continuous(expand = c(0, 0),sec.axis = sec_axis(~./fact, name = "Residual Std. Error (Sigma)")) + 
    theme(axis.text.x = element_text(angle=80,size = 10,hjust=0))
  
  return(list(R1,R2))
}

Ho1_R <- list(Comp_I_IS_Ho1[[4]],Comp_I_RI_Ho1[[1]],Comp_I_RI_Ho1[[2]],Comp_I_RI_Ho1[[3]])
names(Ho1_R) <- c(names(Comp_I_IS_Ho1)[4],names(Comp_I_RI_Ho1)[1:3])
ML_Resum(Ho1_R,"Ho1")

## [[1]]
##                  Sigma     AIC     BIC
## Ho1--OLS I base 0.5796 76.3682 79.8436
## Ho1--OLS I RI   0.5987 67.2228 68.8063
## Ho1--MLE I RI   0.5727 78.3682 83.5812
## Ho1--REML I RI  0.5796 81.3709 86.5116
## 
## [[2]]

Ho2_R <- list(Comp_I_IS_Ho2[[4]],Comp_I_RI_Ho2[[1]],Comp_I_RI_Ho2[[2]],Comp_I_RI_Ho2[[3]])
names(Ho2_R) <- c(names(Comp_I_IS_Ho2)[4],names(Comp_I_RI_Ho2)[1:3])
ML_Resum(Ho2_R,"Ho2")

## [[1]]
##                  Sigma      AIC      BIC
## Ho2--OLS I base 0.2262  -2.6861   0.7893
## Ho2--OLS I RI   0.0635 -94.3562 -92.7727
## Ho2--MLE I RI   0.0635 -80.0659 -74.8528
## Ho2--REML I RI  0.0635 -77.1316 -71.9909
## 
## [[2]]

Com s’ha comentat en l’anterior càlcul, es denoten aquí encara més les diferents característiques dels datasets:

• Ho1: S’observen resultats dels indicadors molt similars i que simplement a nivell de comparació podria semblar que el model ajustat per ANOVA clàssic seria. el millor, els errors estàndards residuals queden molt igualades. Al ser un model homogeni no es preveu necessari utilitzar models complexos com els de màxima versemblança essent correctes els models calculats de la manera clàssica.
• Ho2: Tant l’error estàndard residual com els indicadors d’ajust del model AIC i BIC donen a entendre les diferències entre el model bàsic i els models amb efectes aleatoris que semblen ser clau per millorar l’ajust del model. En aquest cas concret no hi ha grans diferències entr els 3 models amb efecte aleatori RI, tot i que a la vista dels gràfics podria semblar que el millor model seria l’ajust pels mètodes clàssics amb efecte aleatori.

S’extreu a partir dels models les estimacions de la matriu D per veure l’error estàndard dels efectes aleatoris i els residus, i la matriu R on es veu com es distribueix l’error estàndard residual en el model:

Ho1_B <- list(Comp_I_RI_Ho1[[5]])
names(Ho1_B) <- names(Comp_I_RI_Ho1)[5]

Ho2_B <- list(Comp_I_RI_Ho2[[5]])
names(Ho1_B) <- names(Comp_I_RI_Ho2)[5]

lapply(1:length(Ho1_B),function(B){
  VarCorr(Ho1_B[[B]])
})

## [[1]]
## Lot = pdLogChol(1) 
##             Variance     StdDev      
## (Intercept) 1.291945e-10 1.136638e-05
## Residual    3.279731e-01 5.726893e-01

lapply(1:length(Ho2_B),function(B){
  VarCorr(Ho2_B[[B]])
})

## [[1]]
## Lot = pdLogChol(1) 
##             Variance    StdDev    
## (Intercept) 0.045905006 0.21425454
## Residual    0.004028186 0.06346799

lapply(1:length(Ho1_B),function(B){
  getVarCov(Ho1_B[[B]],type="conditional")
})

## [[1]]
## Lot Lot 1 
## Conditional variance covariance matrix
##         1       2       3       4       5       6       7
## 1 0.32797 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## 2 0.00000 0.32797 0.00000 0.00000 0.00000 0.00000 0.00000
## 3 0.00000 0.00000 0.32797 0.00000 0.00000 0.00000 0.00000
## 4 0.00000 0.00000 0.00000 0.32797 0.00000 0.00000 0.00000
## 5 0.00000 0.00000 0.00000 0.00000 0.32797 0.00000 0.00000
## 6 0.00000 0.00000 0.00000 0.00000 0.00000 0.32797 0.00000
## 7 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.32797
##   Standard Deviations: 0.57269 0.57269 0.57269 0.57269 0.57269 0.57269 0.57269

lapply(1:length(Ho2_B),function(B){
  getVarCov(Ho2_B[[B]],type="conditional")
})

## [[1]]
## Lot Lot 1 
## Conditional variance covariance matrix
##           1         2         3         4         5         6         7
## 1 0.0040282 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 2 0.0000000 0.0040282 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 3 0.0000000 0.0000000 0.0040282 0.0000000 0.0000000 0.0000000 0.0000000
## 4 0.0000000 0.0000000 0.0000000 0.0040282 0.0000000 0.0000000 0.0000000
## 5 0.0000000 0.0000000 0.0000000 0.0000000 0.0040282 0.0000000 0.0000000
## 6 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0040282 0.0000000
## 7 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0040282
##   Standard Deviations: 0.063468 0.063468 0.063468 0.063468 0.063468 0.063468 0.063468

lapply(1:length(Ho1_B),function(B){
  a<-getVarCov(Ho1_B[[B]],type="marginal")
  b<-cov2cor(a[[1]])
  return(list(a,b))
})

## [[1]]
## [[1]][[1]]
## Lot Lot 1 
## Marginal variance covariance matrix
##            1          2          3          4          5          6
## 1 3.2797e-01 1.2919e-10 1.2919e-10 1.2919e-10 1.2919e-10 1.2919e-10
## 2 1.2919e-10 3.2797e-01 1.2919e-10 1.2919e-10 1.2919e-10 1.2919e-10
## 3 1.2919e-10 1.2919e-10 3.2797e-01 1.2919e-10 1.2919e-10 1.2919e-10
## 4 1.2919e-10 1.2919e-10 1.2919e-10 3.2797e-01 1.2919e-10 1.2919e-10
## 5 1.2919e-10 1.2919e-10 1.2919e-10 1.2919e-10 3.2797e-01 1.2919e-10
## 6 1.2919e-10 1.2919e-10 1.2919e-10 1.2919e-10 1.2919e-10 3.2797e-01
## 7 1.2919e-10 1.2919e-10 1.2919e-10 1.2919e-10 1.2919e-10 1.2919e-10
##            7
## 1 1.2919e-10
## 2 1.2919e-10
## 3 1.2919e-10
## 4 1.2919e-10
## 5 1.2919e-10
## 6 1.2919e-10
## 7 3.2797e-01
##   Standard Deviations: 0.57269 0.57269 0.57269 0.57269 0.57269 0.57269 0.57269 
## 
## [[1]][[2]]
##             1           2           3           4           5           6
## 1 1.00000e+00 3.93918e-10 3.93918e-10 3.93918e-10 3.93918e-10 3.93918e-10
## 2 3.93918e-10 1.00000e+00 3.93918e-10 3.93918e-10 3.93918e-10 3.93918e-10
## 3 3.93918e-10 3.93918e-10 1.00000e+00 3.93918e-10 3.93918e-10 3.93918e-10
## 4 3.93918e-10 3.93918e-10 3.93918e-10 1.00000e+00 3.93918e-10 3.93918e-10
## 5 3.93918e-10 3.93918e-10 3.93918e-10 3.93918e-10 1.00000e+00 3.93918e-10
## 6 3.93918e-10 3.93918e-10 3.93918e-10 3.93918e-10 3.93918e-10 1.00000e+00
## 7 3.93918e-10 3.93918e-10 3.93918e-10 3.93918e-10 3.93918e-10 3.93918e-10
##             7
## 1 3.93918e-10
## 2 3.93918e-10
## 3 3.93918e-10
## 4 3.93918e-10
## 5 3.93918e-10
## 6 3.93918e-10
## 7 1.00000e+00

lapply(1:length(Ho2_B),function(B){
  a<-getVarCov(Ho2_B[[B]],type="marginal")
  b<-cov2cor(a[[1]])
  return(list(a,b))
})

## [[1]]
## [[1]][[1]]
## Lot Lot 1 
## Marginal variance covariance matrix
##          1        2        3        4        5        6        7
## 1 0.049933 0.045905 0.045905 0.045905 0.045905 0.045905 0.045905
## 2 0.045905 0.049933 0.045905 0.045905 0.045905 0.045905 0.045905
## 3 0.045905 0.045905 0.049933 0.045905 0.045905 0.045905 0.045905
## 4 0.045905 0.045905 0.045905 0.049933 0.045905 0.045905 0.045905
## 5 0.045905 0.045905 0.045905 0.045905 0.049933 0.045905 0.045905
## 6 0.045905 0.045905 0.045905 0.045905 0.045905 0.049933 0.045905
## 7 0.045905 0.045905 0.045905 0.045905 0.045905 0.045905 0.049933
##   Standard Deviations: 0.22346 0.22346 0.22346 0.22346 0.22346 0.22346 0.22346 
## 
## [[1]][[2]]
##           1         2         3         4         5         6         7
## 1 1.0000000 0.9193285 0.9193285 0.9193285 0.9193285 0.9193285 0.9193285
## 2 0.9193285 1.0000000 0.9193285 0.9193285 0.9193285 0.9193285 0.9193285
## 3 0.9193285 0.9193285 1.0000000 0.9193285 0.9193285 0.9193285 0.9193285
## 4 0.9193285 0.9193285 0.9193285 1.0000000 0.9193285 0.9193285 0.9193285
## 5 0.9193285 0.9193285 0.9193285 0.9193285 1.0000000 0.9193285 0.9193285
## 6 0.9193285 0.9193285 0.9193285 0.9193285 0.9193285 1.0000000 0.9193285
## 7 0.9193285 0.9193285 0.9193285 0.9193285 0.9193285 0.9193285 1.0000000

Les matrius s’extreuen del model REML i es consideren que no variarien gaire en els altres models provats de efectes aleatoris:

• En els primers resultat s’observen a les matrius D altra vegada les diferències entre l’variància residual i l’efecte aleatori com canvia entre un conjunt i l’altra.
• Els resultats 3 i 4 on s’extreuen les matrius R que reflexen la matriu de variàncies-covariàncies residuals dels diferents temps que tenen tots els lots. S’està en el cas on es considera la variància homogènia i igual en tots els punts de temps.
• Els resultats 5 i 6 reflexen:
  • La matriu tipus marginal on es suma la variància i la variància degut a l’efecte aleatori en la diagonal, mentre que la que és degut a l’efecte aleatori forma part de la covariància en aquest cas. Això és degut a que en aquesta matriu es reflexa l’efecte de la variància total dins un dels individus (equivalent per tots els altres) que té la de l’efecte aleatori sumat a l’variància residual al centrar-nos en un punt, i es relaciona amb la resta de punts corresponents als altres temps amb el factor calculat aleatori de l’variància dins el mateix lot. És interessant notar també en aquesta matriu que en la diagonal on es sumen els errors estàndards resulta una variància global menor en el model que en teoria està més ben ajustat de Ho2.
  • Es reflecteix també a la segona matriu on es plasma la correlació de variàncies i on es veu que pel primer cas no hi ha correlació de variàncies entre els diferents temps i en canvi en el segon conjunt amb Ho2 sí que es pot veure un nivell clarament alt de correlació.

Es pot graficar el gràfic de distribució dels efectes aleatoris trobats, però com es veu a continuació té poc sentit al tenir un nombre de lots habitualment tant petit:

qqnorm(Comp_I_RI_Ho1[[5]], ~ranef(.))

qqnorm(Comp_I_RI_Ho2[[5]], ~ranef(.))

Finalment si es genera els gràfics bàsics de diagnòstic del model simple:

par(mfrow=c(1,2))
plot(Comp_I_IS_Ho1[[5]],which=c(1:2),main="OLS IS Ho1")

plot(Comp_I_IS_Ho2[[5]],which=c(1:2),main="OLS IS Ho2")

par(mfrow=c(1,1))

par(mfrow=c(1,2))
plot(Comp_I_IS_Ho1[[4]],which=c(1:2),main="OLS I Ho1")

plot(Comp_I_IS_Ho2[[4]],which=c(1:2),main="OLS I Ho2")

par(mfrow=c(1,1))

S’observa que la distribució dels residus és relativament homogènia pels dos casos quadrant amb el que es pretenia generar. En l’anàlisi d’aquest dataset no es considera per tant la inclusió de una funció per modular les variabilitat entre repeticions (o covariables).

PROVES DE VARIACIÓ DE LOTS/TEMPS Ho

Es desenvolupen unes fòrmules personalitzades dins el codi que ens permeti graficar heatmaps dels indicadors sigma, AIC i BIC passant pel rang definit de lots i temps possibles de l’estudi tant pel mètode clàssic ANOVA com pels mètodes de màxima versemblança (ML):

HeatMapRange_OLS <- function(Key,df,fm){

ML <- lapply(i,function(i_){
  lapply(j,function(j_){
    df <- data.frame(p = as.vector(unlist(df[c(1:j_),c(1:i_)])))
    colnames(df)[1] <- Key
    df$Lot <- rep(paste("Lot",c(1:i_),sep="-"),each=j_)
    aov(fm, data = df)$Within
  })
})

sigma <- sapply(1:length(i),function(i_l){
  sapply(1:length(j),function(j_l){
    sigma(summary.lm(ML[[i_l]][[j_l]]))
  })
})

sigma <- data.frame(sigma = as.vector(unlist(sigma)))
sigma$NLots <- rep(i,each=length(j))
sigma$NTemps <- rep(j,length(i))

AIC <- sapply(1:length(i),function(i_l){
  sapply(1:length(j),function(j_l){
    AIC(ML[[i_l]][[j_l]])
  })
})
AIC <- data.frame(AIC = as.vector(unlist(AIC)))
AIC$NLots <- rep(i,each=length(j))
AIC$NTemps <- rep(j,length(i))

BIC <- sapply(1:length(i),function(i_l){
  sapply(1:length(j),function(j_l){
    BIC(ML[[i_l]][[j_l]])
  })
})
BIC <- data.frame(BIC = as.vector(unlist(BIC)))
BIC$NLots <- rep(i,each=length(j))
BIC$NTemps <- rep(j,length(i))

p1<-ggplot(sigma, aes(x = NLots, y = NTemps, fill = sigma)) + geom_tile() + 
  labs(title="Sigma Random ef. OLS HeatM.")
p2<-ggplot(AIC, aes(x = NLots, y = NTemps, fill = AIC)) + geom_tile() + 
  labs(title="AIC Random ef. OLS HeatM.")
p3<-ggplot(BIC, aes(x = NLots, y = NTemps, fill = BIC)) + geom_tile() + 
  labs(title="BIC Random ef. OLS HeatM.")

grid.arrange(p1,p2,p3,ncol=2)
}

HeatMapRange_ML <- function(Key,df,fm,re,met,weight,weightF){ #Weight 0 per no tenir weights, 1 varIdent, 2 varExp, 3 varPower.

if (weight=="0"){
  ML <- lapply(i[1:(ncol(df)-i[1]+1)],function(i_){
  lapply(j[1:(nrow(df)-j[1]+1)],function(j_){
    df <- data.frame(p = as.vector(unlist(df[c(1:j_),c(1:i_)])))
    colnames(df)[1] <- Key
    df$Lot <- rep(paste("Lot",c(1:i_),sep="-"),each=j_)
    df$Temps <- rep(TA[1:j_],i_)
    summary(lme(fm, random = re, data = df, method=met))
  })
})
} else {
  ML <- lapply(i[1:(ncol(df)-i[1]+1)],function(i_){
  lapply(j[1:(nrow(df)-j[1]+1)],function(j_){
    df <- data.frame(p = as.vector(unlist(df[c(1:j_),c(1:i_)])))
    colnames(df)[1] <- Key
    df$Lot <- rep(paste("Lot",c(1:i_),sep="-"),each=j_)
    df$Temps <- rep(TA[1:j_],i_)
    summary(lme(fm, random = re, data = df, method=met, weights = 
                  if (weight=="1"){varIdent(form=weightF)}
                else if (weight=="2"){varExp(form=weightF)}
                else if (weight=="3"){
                  if(weightF=="0"){varPower()} else {varPower(form=weightF)}
                }
                    ))
  })
})
}

sigmaM <- sapply(1:length(i[1:(ncol(df)-i[1]+1)]),function(i_l){
  sapply(1:length(j[1:(nrow(df)-j[1]+1)]),function(j_l){
    ML[[i_l]][[j_l]]$sigma
  })
})

sigmaM <- data.frame(sigma = as.vector(unlist(sigmaM)))
sigmaM$NLots <- rep(i[1:(ncol(df)-i[1]+1)],each=length(j[1:(nrow(df)-j[1]+1)]))
sigmaM$NTemps <- rep(j[1:(nrow(df)-j[1]+1)],length(i[1:(ncol(df)-i[1]+1)]))

AICM <- sapply(1:length(i[1:(ncol(df)-i[1]+1)]),function(i_l){
  sapply(1:length(j[1:(nrow(df)-j[1]+1)]),function(j_l){
    ML[[i_l]][[j_l]]$AIC
  })
})
AICM <- data.frame(AIC = as.vector(unlist(AICM)))
AICM$NLots <- rep(i[1:(ncol(df)-i[1]+1)],each=length(j[1:(nrow(df)-j[1]+1)]))
AICM$NTemps <- rep(j[1:(nrow(df)-j[1]+1)],length(i[1:(ncol(df)-i[1]+1)]))

BICM <- sapply(1:length(i[1:(ncol(df)-i[1]+1)]),function(i_l){
  sapply(1:length(j[1:(nrow(df)-j[1]+1)]),function(j_l){
    ML[[i_l]][[j_l]]$BIC
  })
})
BICM <- data.frame(BIC = as.vector(unlist(BICM)))
BICM$NLots <- rep(i[1:(ncol(df)-i[1]+1)],each=length(j[1:(nrow(df)-j[1]+1)]))
BICM$NTemps <- rep(j[1:(nrow(df)-j[1]+1)],length(i[1:(ncol(df)-i[1]+1)]))

p1<-ggplot(sigmaM, aes(x = NLots, y = NTemps, fill = sigma)) + geom_tile() + 
  labs(title="Sigma Random ef. ML HeatM.")
p2<-ggplot(AICM, aes(x = NLots, y = NTemps, fill = AIC)) + geom_tile() + 
  labs(title="AIC Random ef. ML HeatM.")
p3<-ggplot(BICM, aes(x = NLots, y = NTemps, fill = BIC)) + geom_tile() + 
  labs(title="BIC Random ef. ML HeatM.")

grid.arrange(p1,p2,p3,ncol=2)
}

set.seed(7982)
i <- seq(i3,i10)
j <- seq(j4,j10)

DS_Ho1_df <- sapply(1:i[length(i)],function(a){
  runif(n=j[length(j)], min=-1, max=1)
}) #Es considera un Dataset amb el màxim de dades i després s'extreuen les que interessin a cada punt del rang a testar

fm_rOLS_Ho1 <- formula(Ho1 ~ 1 + Error(Lot))

HeatMapRange_OLS("Ho1",DS_Ho1_df,fm_rOLS_Ho1)

fm_ML_Ho1 <-  formula(Ho1 ~ 1)
re_RI <- formula(~1|Lot)
meMLE <- "ML"

HeatMapRange_ML("Ho1",DS_Ho1_df,fm_ML_Ho1,re_RI,meMLE,0)

meREML <- "REML"
HeatMapRange_ML("Ho1",DS_Ho1_df,fm_ML_Ho1,re_RI,meREML,0)

set.seed(7982)

DS_Ho2_df <- sapply(1:i[length(i)],function(a){
  m_lot <- runif(1, min=-1, max=1)
  runif(n=j[length(j)], min=m_lot-0.1, max=m_lot+0.1)
})

fm_rOLS_Ho2 <-  formula(Ho2 ~ 1 + Error(Lot))

HeatMapRange_OLS("Ho2",DS_Ho2_df,fm_rOLS_Ho2)

fm_ML_Ho2 <-  formula(Ho2 ~ 1)

HeatMapRange_ML("Ho2",DS_Ho2_df,fm_ML_Ho2,re_RI,meMLE,0)

HeatMapRange_ML("Ho2",DS_Ho2_df,fm_ML_Ho2,re_RI,meREML,0)

S’observa que per l’error estàndard residual sembla no haver un patró concret per tenir una variància residual més baixa o més alta possiblement pel fet de treballar amb una simulació de dades tipus núvol de punts en el cas Ho1 o els lots en posicions aleatòries en Ho2. Sí que s’observa clarament com s’inverteix el comportament pel cas dels AIC i BIC segons el tipus de datset que S’està tractant:

• Ho1: Al contrari del que habitualment s’esperaria, quan el model disminueix la mida tant en lots com en temps es calcula uns factors més baixos i per tant suposaria un millor ajust per les dades. Probablement es deu també al tipus de simulació de dataset utilitzat on un major nombre de punts simplement afegeix més variància al tractar-se d’un núvol de punts.
• Ho2: AIC i BIC disminueixen significativament quan s’augmenta el nombre de lots o de punts d’anàlisi possiblement degut a que un major nombre de mostres permet ajustar amb més precisió l’efecte aleatori provocat per la variació dels lots.

MODEL MIXT - DATASET HOMOGENI He

GENERACIÓ DATASET

S’aprofita la matriu base generada per generar les noves variables utilitzant sempre la mateixa llavor 7982.

GENERACIÓ DATASETS He1/He2/He3 (HETEROGENI DEPENENT DEL LOT)

En els dos datasets Ho s’assumeix la homogeneitat residual (una de les assumpcions habituals per poder ajustar els models lineals), pel que és interessant generar datasets pel cas que no es compleixi aquesta condició. Es simula un dels casos més habituals que és l’augment de l’error estàndard en funció del temps. Per fer això es considera un rang de variància variable entre medicions de cada lot. Es té en compte 3 graus de variació:

• Dataset He1 més simple on es depengui del temps multiplicat per 2.
• Dataset He2 on depengui del temps al quadrat.
• Dataset He3 més complex on intervingui una funció exponencial.

Per tots els casos es continua considerant una variació de les mitges de cada lot aleatòria entre un ordre superior de magnitud (-10 i 10). Per fer això es fa servir la funció runif combinada amb alguna de les funcions de la família apply per modificar el rang de cada lot i temps de manera seqüencial:

set.seed(7982)
Av_DS$He1 <- as.vector(sapply(1:av_i, function(i){
  m_lot <- runif(1, min=-10, max=10)
  sapply(1:av_j, function(j){
    runif(n=1, min=m_lot-(j*2), max=m_lot+(j*2))
  })
  }))
kable(head(Av_DS[,c(colnames(Av_DS)[1:3],"He1")]), booktabs = T) %>% kable_styling(latex_options =  c("striped","condensed"))

Lot	NTemps	Temps	He1
Lot 1	Temps 1	0	-1.2864362
Lot 1	Temps 2	3	0.7543398
Lot 1	Temps 3	6	-4.1725533
Lot 1	Temps 4	9	2.2019265
Lot 1	Temps 5	12	-2.7977619
Lot 1	Temps 6	18	-6.2111951

set.seed(7982)
Av_DS$He2 <- as.vector(sapply(1:av_i, function(i){
  m_lot <- runif(1, min=-10, max=10)
  sapply(1:av_j, function(j){
    runif(n=1, min=m_lot-(j^2), max=m_lot+(j^2))
  })
  }))
kable(head(Av_DS[,c(colnames(Av_DS)[1:3],"He2")]), booktabs = T) %>% kable_styling(latex_options =  c("striped","condensed"))

Lot	NTemps	Temps	He2
Lot 1	Temps 1	0	-1.4896161
Lot 1	Temps 2	3	0.7543398
Lot 1	Temps 3	6	-5.4124321
Lot 1	Temps 4	9	6.0966489
Lot 1	Temps 5	12	-4.4552110
Lot 1	Temps 6	18	-15.2479934

set.seed(7982)
Av_DS$He3 <- as.vector(sapply(1:av_i, function(i){
  m_lot <- runif(1, min=-10, max=10)
  sapply(1:av_j, function(j){
    runif(n=1, min=m_lot-(exp(j)), max=m_lot+(exp(j)))
  })
  }))
kable(head(Av_DS[,c(colnames(Av_DS)[1:3],"He3")]), booktabs = T) %>% kable_styling(latex_options =  c("striped","condensed"))

Lot	NTemps	Temps	He3
Lot 1	Temps 1	0	-1.140496
Lot 1	Temps 2	3	2.827710
Lot 1	Temps 3	6	-9.994006
Lot 1	Temps 4	9	24.887784
Lot 1	Temps 5	12	-18.091946
Lot 1	Temps 6	18	-153.597157

VISUALITZACIÓ DATASETS GENERATS He

Es visualitza gràficament les diferències entre els datasets generats. Per veure amb més claredat les diferències es visualitza primer mitjançant un gràfic tipus scatterplot on es veu l’error estàndard global i la diferenciació per lot, i posteriorment es visualitza mitjançant els diagrames de caixes per apreciar més fàcilment els errors estàndards en el temps:

SCATTERPLOTS He

es veu els datasets amb variació de variància en el temps:

labsHe_1 <- c("Dataset He1 - Factorx2","Dataset He2 - Factor quadràtic","Dataset He3 - Factor exponencial")
labsHe_2 <- c("Dataset He1 s/ lot","Dataset He2 s/ lot","Dataset He3 s/ lot")

HeVars <- c("He1","He2","He3")

Av_DS_plotHe_sc1 <- lapply(1:length(HeVars),function(graph){
  ggplot(Av_DS, aes(x=Temps, y=Av_DS[,HeVars[graph]])) + geom_point() + geom_smooth(linetype="dashed") + 
    labs(title=labsHe_1[graph], y = "Observacions")
  })

Av_DS_plotHe_sc2 <- lapply(1:length(HeVars),function(graph){
  ggplot(Av_DS, aes(x=Temps, y=Av_DS[,HeVars[graph]], color=Lot, group=Lot)) + geom_point() + geom_line() + 
    labs(title=labsHe_2[graph], y = "Observacions")
  })

ydens_DS_plotHe <- lapply(1:length(HeVars),function(graph){
  ggplot(Av_DS, aes(x=Av_DS[,HeVars[graph]])) + geom_density(color="darkblue", fill="lightblue") +
    geom_vline(aes(xintercept=mean(Av_DS[,HeVars[graph]])), color="blue", linetype="dashed", size=1)+ 
    coord_flip() + labs(x = "Observacions")
  })

grid.arrange(Av_DS_plotHe_sc1[[1]], ydens_DS_plotHe[[1]], Av_DS_plotHe_sc1[[2]], ydens_DS_plotHe[[2]],
             Av_DS_plotHe_sc1[[3]], ydens_DS_plotHe[[3]],ncol=2, widths=c(4, 1))

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

grid.arrange(Av_DS_plotHe_sc2[[1]], ydens_DS_plotHe[[1]], Av_DS_plotHe_sc2[[2]], ydens_DS_plotHe[[2]],
             Av_DS_plotHe_sc2[[3]], ydens_DS_plotHe[[3]],ncol=2, widths=c(4, 1))

Els gràfics de densitat ens mostren que al tenir més dispersió de dades hi ha una major concentració de punts en la mitja i que aquest efecte s’intensifica quan hi ha una progressió de variància més alta. Per la diferència d’escala i el tipus de gràfic no s’aprecia tant clar les diferències de variància, però sí que s’observa amb relativa claredat que la evolució en el temps és de diferent magnitud.

BOXPLOTS He

es veu els datasets amb variació de variància en el temps:

Av_DS_plotHe_bp1 <- lapply(1:length(HeVars),function(graph){
  ggplot(Av_DS, aes(x=NTemps, y=Av_DS[,HeVars[graph]])) + 
  geom_boxplot()+ geom_dotplot(binaxis='y', stackdir='center', dotsize=0.5)+
  labs(title=labsHe_1[graph], y = "Observacions")
  })

lapply(1:3, function(ga){
  grid.arrange(Av_DS_plotHe_bp1[[ga]], ydens_DS_plotHe[[ga]], 
               ncol=2, widths=c(4, 1))
})

## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

## [[1]]
## TableGrob (1 x 2) "arrange": 2 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 
## [[2]]
## TableGrob (1 x 2) "arrange": 2 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 
## [[3]]
## TableGrob (1 x 2) "arrange": 2 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]

S’observa clarament que en aquests 3 datasets hi ha una canvi de variància al augmentar el temps que és el que es buscava generar en la simulació.

MODELS LINEALS He

Com en el cas de l’anàlisi de Ho es pretén generar un model lineal simple i un model lineal d’efectes mixtos amb la intercepció aleatòria, afegint per aquest cas models amb variables per explicar el canvi de variància en el temps.

ML I/IS He

Es genera el model lineal simple i el model amb la variable temps per comparar-los i definir el model de base simple (s’aprofita la funció ja creada):

(Comp_I_IS_He1 <- Comp_I_IS_LM("He1",Av_DS))[1:3]

## $`OLS I`
## 
## Call:
## lm(formula = paste(var, "~", 1, sep = " "), data = ds)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17.2614  -1.9536   0.5654   3.5592  13.5827 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -4.3155     0.9737  -4.432 6.81e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## s: 6.31 on 41 degrees of freedom
## 
## 
## $`OLS IS`
## 
## Call:
## lm(formula = paste(var, "~", "Temps", sep = " "), data = ds)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.4919  -3.0986   0.0871   3.5572  15.1406 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  -2.2384     1.5739  -1.422    0.163
## Temps        -0.2019     0.1217  -1.659    0.105
## 
## s: 6.179 on 40 degrees of freedom
## Multiple R-squared: 0.06437,
## Adjusted R-squared: 0.04098 
## F-statistic: 2.752 on 1 and 40 DF,  p-value: 0.105 
## 
## 
## $`Anova comp. efecte S`
## Analysis of Variance Table
## 
## Model 1: He1 ~ Temps
## Model 2: He1 ~ 1
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     40 1527.4                           
## 2     41 1632.5 -1   -105.08 2.7518  0.105

(Comp_I_IS_He2 <- Comp_I_IS_LM("He2",Av_DS))[1:3]

## $`OLS I`
## 
## Call:
## lm(formula = paste(var, "~", 1, sep = " "), data = ds)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -49.843  -2.681   1.519   6.716  38.440 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   -6.060      2.639  -2.297   0.0268 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## s: 17.1 on 41 degrees of freedom
## 
## 
## $`OLS IS`
## 
## Call:
## lm(formula = paste(var, "~", "Temps", sep = " "), data = ds)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -41.040  -5.734  -1.675   8.464  43.391 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   0.5415     4.2100   0.129   0.8983  
## Temps        -0.6418     0.3256  -1.971   0.0557 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## s: 16.53 on 40 degrees of freedom
## Multiple R-squared: 0.08853,
## Adjusted R-squared: 0.06574 
## F-statistic: 3.885 on 1 and 40 DF,  p-value: 0.05567 
## 
## 
## $`Anova comp. efecte S`
## Analysis of Variance Table
## 
## Model 1: He2 ~ Temps
## Model 2: He2 ~ 1
##   Res.Df   RSS Df Sum of Sq     F  Pr(>F)  
## 1     40 10929                             
## 2     41 11990 -1   -1061.5 3.885 0.05567 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Comp_I_IS_He3 <- Comp_I_IS_LM("He3",Av_DS))[1:3]

## $`OLS I`
## 
## Call:
## lm(formula = paste(var, "~", 1, sep = " "), data = ds)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1019.79    37.40    58.94    66.83   449.78 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   -63.57      42.56  -1.493    0.143
## 
## s: 275.8 on 41 degrees of freedom
## 
## 
## $`OLS IS`
## 
## Call:
## lm(formula = paste(var, "~", "Temps", sep = " "), data = ds)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -827.75  -74.18   -8.00   71.26  572.91 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   80.467     65.116   1.236  0.22376   
## Temps        -14.003      5.037  -2.780  0.00824 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## s: 255.7 on 40 degrees of freedom
## Multiple R-squared: 0.162,
## Adjusted R-squared: 0.141 
## F-statistic:  7.73 on 1 and 40 DF,  p-value: 0.00824 
## 
## 
## $`Anova comp. efecte S`
## Analysis of Variance Table
## 
## Model 1: He3 ~ Temps
## Model 2: He3 ~ 1
##   Res.Df     RSS Df Sum of Sq      F  Pr(>F)   
## 1     40 2614521                               
## 2     41 3119769 -1   -505248 7.7299 0.00824 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Al haver simulat una variància que varia amb el temps al voltant de la mateixa mitja s’esperaria que el predictor Temps no fos significatiu per explicar la mitja de la resposta. Tot i així, al tenir un nombre de lots limitat, provoca que l’augment de variabilitat emmascari la mitja real sobre la qual es reparteixen i dona un fals efecte de tendència. És per aquest motiu probablement que s’observa que afegir la variable temps al model pren més significació al tenir un model amb un factor d’augment de variància més gran:

• Model He1 amb factor multiplicador dona un contrast no significatiu pel Temps, però ja més significatiu que el model anterior homogeni (p-valors: 0.105 (cas He1), 0.5612 (cas Ho1) i 0.5883 (cas Ho2)).
• Model He1 no significatiu, però al límit de la significació si es pren una significació límit de 0.05: p-valor de 0.0557.
• Model He2 significatiu: p-valor de 0.0082.

En aquest cas no seria tan senzill com descartar el model amb l’efecte fix de covariància del predictor Temps, però és interessant veure si aquest fals efecte de tendència en els propers models és necessari mantenir-lo o si es podrà explicar el model amb les variables d’efectes aleatoris i/o funcions de variància-covariància entre repeticions.

S’observa també com ha resultat l’error estimatiu i l’error estàndard residual del model segons el dataset :

Models I:

I_HE_res <- data.frame(Dataset = HeVars)
I_HE_res$"Error estimació" <- round(c(Comp_I_IS_He1[[1]]$coefficients[,2],
                                Comp_I_IS_He2[[1]]$coefficients[,2],Comp_I_IS_He3[[1]]$coefficients[,2]),4)
I_HE_res$"Error residual" <- round(c(Comp_I_IS_He1[[1]]$sigma,
                                Comp_I_IS_He2[[1]]$sigma,Comp_I_IS_He3[[1]]$sigma),4)

kable(I_HE_res, booktabs = T) %>% kable_styling(latex_options =  c("striped","condensed"))

Dataset	Error estimació	Error residual
He1	0.9737	6.3101
He2	2.6387	17.1010
He3	42.5642	275.8477

Models IS:

IS_HE_res <- data.frame(Dataset = HeVars)
IS_HE_res$"Error estimació intercepció" <- round(c(Comp_I_IS_He1[[2]]$coefficients[1,2],
                                Comp_I_IS_He2[[2]]$coefficients[1,2],Comp_I_IS_He3[[2]]$coefficients[1,2]),4)
IS_HE_res$"Error estimació predictor" <- round(c(Comp_I_IS_He1[[2]]$coefficients[2,2],
                                Comp_I_IS_He2[[2]]$coefficients[2,2],Comp_I_IS_He3[[2]]$coefficients[2,2]),4)
IS_HE_res$"Error residual" <- round(c(Comp_I_IS_He1[[2]]$sigma,
                                Comp_I_IS_He2[[2]]$sigma,Comp_I_IS_He3[[2]]$sigma),4)

kable(IS_HE_res, booktabs = T) %>% kable_styling(latex_options =  c("striped","condensed"))

Dataset	Error estimació intercepció	Error estimació predictor	Error residual
He1	1.5739	0.1217	6.1794
He2	4.2100	0.3256	16.5293
He3	65.1162	5.0367	255.6619

L’estimació i errors residuals augmenten a mesura que s’estudien models més dispersos ja que van lligats a la incertesa de predicció de les dades, però el més interessant en aquestes últimes taules és veure que tot i augmentar la significació del predictor, aquest no té la força suficient per provocar un canvi molt significatiu si es comparen els errors estàndards residuals dels mateixos datasets en models amb i sense la pendent S, cosa que ens donaria pistes en cas que fos un dataset desconegut de que hi ha algun factor addicional amagat en el model.

ML RI He

Avançant a partir del model simple amb intercepció ML I es fa l’ajust de models amb efecte aleatori (s’utilitza la funció generada). Es genera una funció addicional per fer el mateix a partir de la funció de model simple ML IS per veure també aquest ajust.

Comp_IS_RI_LM <- function(var,ds) {
  fm1 <- formula(paste(var,"~","Temps","+","Error(Lot)"))
  ML_IS_RI_sOLS <- summary(ML_IS_RI_OLS <- aov(fm1, data = ds))
  ML_IS_RI_OLS
  ML_IS_RI_OLS_err <- sigma(summary.lm(ML_IS_RI_OLS$Within))
  
  fm2 <- formula(paste(var,"~","Temps"))
  
  ML_IS_RI_sMLE <- summary(ML_IS_RI_MLE <- lme(fm2, random = ~1|Lot, data = ds, method="ML"))
  
  ML_IS_RI_sREML <- summary(ML_IS_RI_REML <- lme(fm2,random = ~1|Lot, data = ds, method="REML"))
  
    l <- list(ML_IS_RI_OLS$Within,ML_IS_RI_sMLE,ML_IS_RI_sREML,ML_IS_RI_OLS,ML_IS_RI_MLE,ML_IS_RI_REML)
  names(l) <- c("OLS IS RI","MLE IS RI","REML IS RI","OLS IS RI base","MLE IS RI base","REML IS RI base")
  
  return(l)
} 

(Comp_I_RI_He1 <- Comp_I_RI_LM("He1",Av_DS))[1:3]

## $`OLS I RI`
## 
## Terms:
##                 Residuals
## Sum of Squares   1175.276
## Deg. of Freedom        36
## 
## Residual standard error: 5.713716
## 
## $`MLE I RI`
## Linear mixed-effects model fit by maximum likelihood
##  Data: ds 
##       AIC     BIC   logLik
##   276.678 281.891 -135.339
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:    2.494505 5.713716
## 
## Fixed effects: list(fm2) 
##                Value Std.Error DF   t-value p-value
## (Intercept) -4.31553  1.363321 36 -3.165453  0.0031
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.4204112 -0.4077346  0.1686668  0.5467930  2.0070963 
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`REML I RI`
## Linear mixed-effects model fit by REML
##  Data: ds 
##        AIC      BIC   logLik
##   274.1559 279.2967 -134.078
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:    2.898365 5.713696
## 
## Fixed effects: list(fm2) 
##                Value Std.Error DF   t-value p-value
## (Intercept) -4.31553  1.475595 36 -2.924603  0.0059
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.3453677 -0.4085803  0.1651534  0.5794877  1.9608561 
## 
## Number of Observations: 42
## Number of Groups: 6

fmHe1 <- formula(paste("He1","~",1))
(eMLE_I_He1<-exactRLRT(lme(fmHe1, random = ~1|Lot, data = Av_DS, method="ML")))

## Using restricted likelihood evaluated at ML estimators.

## Refit with method="REML" for exact results.

## Warning in model.matrix.default(~m$groups[[n.levels - i + 1]] - 1,
## contrasts.arg = c("contr.treatment", : non-list contrasts argument ignored

## 
##  simulated finite sample distribution of RLRT.
##  
##  (p-value based on 10000 simulated values)
## 
## data:  
## RLRT = 2.9145, p-value = 0.033

(eREML_I_He1<-exactRLRT(lme(fmHe1, random = ~1|Lot, data = Av_DS, method="REML")))

## Warning in model.matrix.default(~m$groups[[n.levels - i + 1]] - 1,
## contrasts.arg = c("contr.treatment", : non-list contrasts argument ignored

## 
##  simulated finite sample distribution of RLRT.
##  
##  (p-value based on 10000 simulated values)
## 
## data:  
## RLRT = 2.9909, p-value = 0.0296

(Comp_I_RI_He2 <- Comp_I_RI_LM("He2",Av_DS))[1:3]

## $`OLS I RI`
## 
## Terms:
##                 Residuals
## Sum of Squares   10013.48
## Deg. of Freedom        36
## 
## Residual standard error: 16.6779
## 
## $`MLE I RI`
## Linear mixed-effects model fit by maximum likelihood
##  Data: ds 
##        AIC      BIC    logLik
##   362.5896 367.8026 -178.2948
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:    2.707311  16.6779
## 
## Fixed effects: list(fm2) 
##                 Value Std.Error DF   t-value p-value
## (Intercept) -6.060349  2.834711 36 -2.137907  0.0394
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.88054178 -0.18778591  0.09895573  0.33041501  2.23068218 
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`REML I RI`
## Linear mixed-effects model fit by REML
##  Data: ds 
##        AIC      BIC    logLik
##   358.6036 363.7443 -176.3018
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:    4.091778  16.6779
## 
## Fixed effects: list(fm2) 
##                 Value Std.Error DF   t-value p-value
## (Intercept) -6.060349   3.06808 36 -1.975291  0.0559
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.7829585 -0.2167656  0.1280920  0.3463577  2.1636963 
## 
## Number of Observations: 42
## Number of Groups: 6

fmHe2 <- formula(paste("He2","~",1))
(eMLE_I_He2<-exactRLRT(lme(fmHe2, random = ~1|Lot, data = Av_DS, method="ML")))

## Using restricted likelihood evaluated at ML estimators.
## Refit with method="REML" for exact results.

## Warning in model.matrix.default(~m$groups[[n.levels - i + 1]] - 1,
## contrasts.arg = c("contr.treatment", : non-list contrasts argument ignored

## 
##  simulated finite sample distribution of RLRT.
##  
##  (p-value based on 10000 simulated values)
## 
## data:  
## RLRT = 0.21998, p-value = 0.2656

(eREML_I_He2<-exactRLRT(lme(fmHe2, random = ~1|Lot, data = Av_DS, method="REML")))

## Warning in model.matrix.default(~m$groups[[n.levels - i + 1]] - 1,
## contrasts.arg = c("contr.treatment", : non-list contrasts argument ignored

## 
##  simulated finite sample distribution of RLRT.
##  
##  (p-value based on 10000 simulated values)
## 
## data:  
## RLRT = 0.29637, p-value = 0.2431

(Comp_I_RI_He3 <- Comp_I_RI_LM("He3",Av_DS))[1:3]

## $`OLS I RI`
## 
## Terms:
##                 Residuals
## Sum of Squares    2793880
## Deg. of Freedom        36
## 
## Residual standard error: 278.5817
## 
## $`MLE I RI`
## Linear mixed-effects model fit by maximum likelihood
##  Data: ds 
##       AIC     BIC   logLik
##   596.246 601.459 -295.123
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:  0.02002761  272.544
## 
## Fixed effects: list(fm2) 
##                 Value Std.Error DF  t-value p-value
## (Intercept) -63.56751  42.56422 36 -1.49345   0.144
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.7417541  0.1372290  0.2162638  0.2452097  1.6503027 
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`REML I RI`
## Linear mixed-effects model fit by REML
##  Data: ds 
##        AIC      BIC    logLik
##   586.9182 592.0589 -290.4591
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:  0.02857084 275.8477
## 
## Fixed effects: list(fm2) 
##                 Value Std.Error DF   t-value p-value
## (Intercept) -63.56751  42.56422 36 -1.493449   0.144
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.6969411  0.1355854  0.2136737  0.2422729  1.6305378 
## 
## Number of Observations: 42
## Number of Groups: 6

fmHe3 <- formula(paste("He3","~",1))
(eMLE_I_He3<-exactRLRT(lme(fmHe3, random = ~1|Lot, data = Av_DS, method="ML")))

## Using restricted likelihood evaluated at ML estimators.
## Refit with method="REML" for exact results.

## Warning in model.matrix.default(~m$groups[[n.levels - i + 1]] - 1,
## contrasts.arg = c("contr.treatment", : non-list contrasts argument ignored

## 
##  simulated finite sample distribution of RLRT.
##  
##  (p-value based on 10000 simulated values)
## 
## data:  
## RLRT = 0, p-value = 1

(eREML_I_He3<-exactRLRT(lme(fmHe3, random = ~1|Lot, data = Av_DS, method="REML")))

## Warning in model.matrix.default(~m$groups[[n.levels - i + 1]] - 1,
## contrasts.arg = c("contr.treatment", : non-list contrasts argument ignored

## 
##  simulated finite sample distribution of RLRT.
##  
##  (p-value based on 10000 simulated values)
## 
## data:  
## RLRT = 0, p-value = 1

(Comp_IS_RI_He1 <- Comp_IS_RI_LM("He1",Av_DS))[1:3]

## $`OLS IS RI`
## 
## Terms:
##                     Temps Residuals
## Sum of Squares   105.0784 1070.1974
## Deg. of Freedom         1        35
## 
## Residual standard error: 5.529654
## Estimated effects are balanced
## 
## $`MLE IS RI`
## Linear mixed-effects model fit by maximum likelihood
##  Data: ds 
##        AIC      BIC    logLik
##   275.3062 282.2569 -133.6531
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:    2.576779 5.452302
## 
## Fixed effects: list(fm2) 
##                  Value Std.Error DF   t-value p-value
## (Intercept) -2.2383685 1.7851673 35 -1.253871  0.2182
## Temps       -0.2019462 0.1100661 35 -1.834773  0.0751
##  Correlation: 
##       (Intr)
## Temps -0.634
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.98630953 -0.57575477  0.04054506  0.52365907  2.36305715 
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`REML IS RI`
## Linear mixed-effects model fit by REML
##  Data: ds 
##        AIC    BIC    logLik
##   275.3944 282.15 -133.6972
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:    2.948864 5.529644
## 
## Fixed effects: list(fm2) 
##                  Value Std.Error DF   t-value p-value
## (Intercept) -2.2383685 1.8527912 35 -1.208106  0.2351
## Temps       -0.2019462 0.1089372 35 -1.853786  0.0722
##  Correlation: 
##       (Intr)
## Temps -0.605
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.89801726 -0.56607915  0.03575724  0.52530932  2.29271896 
## 
## Number of Observations: 42
## Number of Groups: 6

fmHe1_S <- formula(paste("He1","~","Temps"))
(eMLE_IS_He1<-exactRLRT(lme(fmHe1_S, random = ~1|Lot, data = Av_DS, method="ML")))

## Using restricted likelihood evaluated at ML estimators.

## Refit with method="REML" for exact results.

## Warning in model.matrix.default(~m$groups[[n.levels - i + 1]] - 1,
## contrasts.arg = c("contr.treatment", : non-list contrasts argument ignored

## 
##  simulated finite sample distribution of RLRT.
##  
##  (p-value based on 10000 simulated values)
## 
## data:  
## RLRT = 3.3569, p-value = 0.0257

(eREML_IS_He1<-exactRLRT(lme(fmHe1_S, random = ~1|Lot, data = Av_DS, method="REML")))

## Warning in model.matrix.default(~m$groups[[n.levels - i + 1]] - 1,
## contrasts.arg = c("contr.treatment", : non-list contrasts argument ignored

## 
##  simulated finite sample distribution of RLRT.
##  
##  (p-value based on 10000 simulated values)
## 
## data:  
## RLRT = 3.4109, p-value = 0.0257

(Comp_IS_RI_He2 <- Comp_IS_RI_LM("He2",Av_DS))[1:3]

## $`OLS IS RI`
## 
## Terms:
##                    Temps Residuals
## Sum of Squares  1061.464  8952.018
## Deg. of Freedom        1        35
## 
## Residual standard error: 15.99287
## Estimated effects are balanced
## 
## $`MLE IS RI`
## Linear mixed-effects model fit by maximum likelihood
##  Data: ds 
##        AIC      BIC    logLik
##   360.5557 367.5064 -176.2779
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:    3.397311 15.76918
## 
## Fixed effects: list(fm2) 
##                  Value Std.Error DF    t-value p-value
## (Intercept)  0.5415064  4.354016 35  0.1243694  0.9017
## Temps       -0.6418471  0.318334 35 -2.0162705  0.0515
##  Correlation: 
##       (Intr)
## Temps -0.752
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.42268652 -0.39505387 -0.08976639  0.45194790  2.62815930 
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`REML IS RI`
## Linear mixed-effects model fit by REML
##  Data: ds 
##        AIC      BIC    logLik
##   357.0559 363.8114 -174.5279
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:    4.465404 15.99287
## 
## Fixed effects: list(fm2) 
##                  Value Std.Error DF    t-value p-value
## (Intercept)  0.5415064  4.462656 35  0.1213417  0.9041
## Temps       -0.6418471  0.315069 35 -2.0371650  0.0493
##  Correlation: 
##       (Intr)
## Temps -0.726
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.31082204 -0.40558024 -0.09537408  0.39209813  2.53787141 
## 
## Number of Observations: 42
## Number of Groups: 6

fmHe2_S <- formula(paste("He2","~","Temps"))
(eMLE_IS_He2<-exactRLRT(lme(fmHe2_S, random = ~1|Lot, data = Av_DS, method="ML")))

## Using restricted likelihood evaluated at ML estimators.
## Refit with method="REML" for exact results.

## Warning in model.matrix.default(~m$groups[[n.levels - i + 1]] - 1,
## contrasts.arg = c("contr.treatment", : non-list contrasts argument ignored

## 
##  simulated finite sample distribution of RLRT.
##  
##  (p-value based on 10000 simulated values)
## 
## data:  
## RLRT = 0.40809, p-value = 0.2143

(eREML_IS_He2<-exactRLRT(lme(fmHe2_S, random = ~1|Lot, data = Av_DS, method="REML")))

## Warning in model.matrix.default(~m$groups[[n.levels - i + 1]] - 1,
## contrasts.arg = c("contr.treatment", : non-list contrasts argument ignored

## 
##  simulated finite sample distribution of RLRT.
##  
##  (p-value based on 10000 simulated values)
## 
## data:  
## RLRT = 0.46211, p-value = 0.2002

(Comp_IS_RI_He3 <- Comp_IS_RI_LM("He3",Av_DS))[1:3]

## $`OLS IS RI`
## 
## Terms:
##                   Temps Residuals
## Sum of Squares   505248   2288632
## Deg. of Freedom       1        35
## 
## Residual standard error: 255.7137
## Estimated effects are balanced
## 
## $`MLE IS RI`
## Linear mixed-effects model fit by maximum likelihood
##  Data: ds 
##        AIC      BIC    logLik
##   590.8256 597.7762 -291.4128
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:  0.02502923 249.5005
## 
## Fixed effects: list(fm2) 
##                 Value Std.Error DF   t-value p-value
## (Intercept)  80.46670  65.11618 35  1.235740  0.2248
## Temps       -14.00333   5.03669 35 -2.780265  0.0087
##  Correlation: 
##       (Intr)
## Temps -0.796
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.31761615 -0.29731165 -0.03205673  0.28562712  2.29621302 
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`REML IS RI`
## Linear mixed-effects model fit by REML
##  Data: ds 
##        AIC     BIC    logLik
##   576.6155 583.371 -284.3077
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:   0.2328601 255.6619
## 
## Fixed effects: list(fm2) 
##                 Value Std.Error DF   t-value p-value
## (Intercept)  80.46670  65.11623 35  1.235740  0.2248
## Temps       -14.00333   5.03669 35 -2.780266  0.0087
##  Correlation: 
##       (Intr)
## Temps -0.796
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.23765998 -0.29014890 -0.03128171  0.27874419  2.24087309 
## 
## Number of Observations: 42
## Number of Groups: 6

fmHe3_S <- formula(paste("He3","~","Temps"))
(eMLE_IS_He3<-exactRLRT(lme(fmHe3_S, random = ~1|Lot, data = Av_DS, method="ML")))

## Using restricted likelihood evaluated at ML estimators.
## Refit with method="REML" for exact results.

## Warning in model.matrix.default(~m$groups[[n.levels - i + 1]] - 1,
## contrasts.arg = c("contr.treatment", : non-list contrasts argument ignored

## 
##  simulated finite sample distribution of RLRT.
##  
##  (p-value based on 10000 simulated values)
## 
## data:  
## RLRT = 0, p-value = 1

(eREML_IS_He3<-exactRLRT(lme(fmHe3_S, random = ~1|Lot, data = Av_DS, method="REML")))

## Warning in model.matrix.default(~m$groups[[n.levels - i + 1]] - 1,
## contrasts.arg = c("contr.treatment", : non-list contrasts argument ignored

## 
##  simulated finite sample distribution of RLRT.
##  
##  (p-value based on 10000 simulated values)
## 
## data:  
## RLRT = 0, p-value = 1

Es genera una taula resum dels resultats obtinguts del test RLRT per començar:

I_IS_RI_HE_RLRT <- data.frame(Dataset = paste(rep(paste(rep(HeVars,each=2),c("I RI","IS RI")),each=3),c("Clàssic OLS","MLE","REML")))

Comp_RI_He_2 <- list("-",eMLE_I_He1,eREML_I_He1,"-",eMLE_IS_He1,eREML_IS_He1,
                     "-",eMLE_I_He2,eREML_I_He2,"-",eMLE_IS_He2,eREML_IS_He2,
                     "-",eMLE_I_He3,eREML_I_He3,"-",eMLE_IS_He3,eREML_IS_He3)
I_IS_RI_HE_RLRT$"p-valor RLRT" <- sapply(1:length(Comp_RI_He_2),function(l1){
  if (class(Comp_RI_He_2[[l1]])!="htest"){
    return("-")
    } else {
      round(Comp_RI_He_2[[l1]]$p.value,4)
    }
  })

kable(I_IS_RI_HE_RLRT, booktabs = T) %>% kable_styling(latex_options =  c("striped","condensed"))

Dataset	p-valor RLRT
He1 I RI Clàssic OLS
He1 I RI MLE	0.033
He1 I RI REML	0.0296
He1 IS RI Clàssic OLS
He1 IS RI MLE	0.0257
He1 IS RI REML	0.0257
He2 I RI Clàssic OLS
He2 I RI MLE	0.2656
He2 I RI REML	0.2431
He2 IS RI Clàssic OLS
He2 IS RI MLE	0.2143
He2 IS RI REML	0.2002
He3 I RI Clàssic OLS
He3 I RI MLE	1
He3 I RI REML	1
He3 IS RI Clàssic OLS
He3 IS RI MLE	1
He3 IS RI REML	1

Es podria considerar una mida mostral correcte la d’aquest model base amb 42 observacions, i entre 1 i 2 efectes fixos (si es considera la mitja com un efecte fix). En aquests resultats com a comentari general es pot apreciar que:

• No sembla haver-hi diferències significatives entre mètodes (MLE; REML).
• En el cas de la utilització de la covariable temps semblaria que al afegir-la disminueix el p-valor calculat de significació de l’efecte aleatori (és a dir, es tornaria més significatiu).
• Entre els tres conjunts de dades donaria la impressió que com més disperses es troben les dades menys significatiu es torna l’efecte aleatori.

ML RI Weights He

#Es genera una funció amb la simple instrucció de quan es facin llistats dintre llistats converteixi tot en una llista normal sense subnivells per facilitar els resums
ConvList <- function(List){
  ListC <- list()
  for(l1 in 1:length(List)) {
    if(class(List[[l1]])=="list"){
    for (l2 in 1:length(List[[l1]])) {
      ListC[[length(ListC)+1]] <- List[[l1]][[l2]]
      names(ListC)[length(ListC)] <- names(List[[l1]])[l2]
      }
    } else {
      ListC[[length(ListC)+1]] <- List[[l1]]
      names(ListC)[length(ListC)] <- names(List)[l1]
    }
  }
    
return(ListC)
}

S’observa els gràfics de diagnòstic de residus dels models amb la covariable temps generats pels conjunts He:

Diag_Hom_He1 <- list(Comp_I_IS_He1[5],Comp_IS_RI_He1[5:6])

Diag_Hom_lHe1 <- ConvList(Diag_Hom_He1)
lapply(1:length(Diag_Hom_lHe1),function(l){
  plot(Diag_Hom_lHe1[[l]],which=1,main=paste("He1",names(Diag_Hom_lHe1)[l],sep="--"))
})

## [[1]]
## NULL
## 
## [[2]]

## 
## [[3]]

Diag_Hom_He2 <- list(Comp_I_IS_He2[5],Comp_IS_RI_He2[5:6])

Diag_Hom_lHe2 <- ConvList(Diag_Hom_He2)
lapply(1:length(Diag_Hom_lHe2),function(l){
  plot(Diag_Hom_lHe2[[l]],which=1,main=paste("He2",names(Diag_Hom_lHe2)[l],sep="--"))
})

## [[1]]
## NULL
## 
## [[2]]

## 
## [[3]]

Diag_Hom_He3 <- list(Comp_I_IS_He3[5],Comp_IS_RI_He3[5:6])

Diag_Hom_lHe3 <- ConvList(Diag_Hom_He3)
lapply(1:length(Diag_Hom_lHe3),function(l){
  plot(Diag_Hom_lHe3[[l]],which=1,main=paste("He3",names(Diag_Hom_lHe3)[l],sep="--"))
})

## [[1]]
## NULL
## 
## [[2]]

## 
## [[3]]

Com estava previst en la simulació s’observa heterogeneitat de residus (amb més claredat com més dispers és el conjunt de dades).

En aquest apartat s’avança amb el conjunts He al intentar aplicar i veure com es comporta el model amb les funcions per modular la variància de l’error en els casos on ja no es considera homogeneitat residual. A nivell de programació és complicat afegir aquests paràmetres de pes als objectes aov pel que tot i que el model ajustat per ANOVA dona més bons resultats que els ML, es veurà les funcions per Heteroscedasticitat aplicat als models ajustats per màxima versemblança més preparats per l’aplicació d’aquestes funcions a nivell d’algoritmes computacionals:

Comp_W_LM <- function(fML,fRE,ds,WVI,WVE,WVP) # WVI;WVE;WVP Per decidir quines funcions s'apliquen essent-> 0:No aplica, 1:Aplicar MLE, 2:Aplicar REML, 3:Aplicar les 2
  {
  
  baseML <- lme(fML, random = fRE, control=ctrl, data = ds,method="ML")
  baseREML <- lme(fML, random = fRE, control=ctrl, data = ds,method="REML")
  l <- list()
  
  if (paste(fML[3])=="1"){
    fMLk <- "I"
  } else if (paste(fML[3])=="Temps"){
    fMLk <- "IS"
  } else {
    fMLk <- "desconegut"
  }
  
  if (paste(fRE[2])=="1 | Lot"){
    fREk <- "RI"
  } else if (paste(fRE[2])=="1 + Temps | Lot"){
    fREk <- "RIS"
  } else if (class(fRE)=="list"){
    fREk <- "DRIS"
  } else {
    fREk <- "desconegut"
  }
  
  
    if(WVI==1 | WVI==3){
      ML_I_RI_WVI_sMLE <- summary(ML_I_RI_WVI_MLE <- update(baseML, weights = varIdent(form = ~Temps)))
      le <- length(l)+1
      l[[le]] <- ML_I_RI_WVI_sMLE
      names(l)[le] <- paste(fMLk,"MLE",fREk,"VarIdent")
      le <- length(l)+1
      l[[le]] <- ML_I_RI_WVI_MLE
      names(l)[le] <- paste(fMLk,"MLE",fREk,"VarIdent Base")
    } else {
      le <- length(l)+1
      l[[le]] <- "Pendent ML_I_RI_WVI_sMLE"
      names(l)[le] <- paste("NE Model",fMLk,"MLE",fREk,"VarIdent")
      le <- length(l)+1
      l[[le]] <- "Pendent ML_I_RI_WVI_MLE"
      names(l)[le] <- paste("NE Model",fMLk,"MLE",fREk,"VarIdent Base")
      }
    if(WVI==2 | WVI==3){
      ML_I_RI_WVI_sREML <- summary(ML_I_RI_WVI_REML <- update(baseREML, weights = varIdent(form = ~Temps)))
      le <- length(l)+1
      l[[le]] <- ML_I_RI_WVI_sREML
      names(l)[le] <- paste(fMLk,"REML",fREk,"VarIdent")
      le <- length(l)+1
      l[[le]] <- ML_I_RI_WVI_REML
      names(l)[le] <- paste(fMLk,"REML",fREk,"VarIdent Base")
    } else {
      le <- length(l)+1
      l[[le]] <- "Pendent ML_I_RI_WVI_sREML"
      names(l)[le] <- paste("NE Model",fMLk,"REML",fREk,"VarIdent")
      le <- length(l)+1
      l[[le]] <- "Pendent ML_I_RI_WVI_REML"
      names(l)[le] <- paste("NE Model",fMLk,"REML",fREk,"VarIdent Base")
    }

    if(WVE==1 | WVE==3){
      ML_I_RI_WVE_sMLE <- summary(ML_I_RI_WVE_MLE <- update(baseML, weights = varExp(form = ~Temps)))
      le <- length(l)+1
      l[[le]] <- ML_I_RI_WVE_sMLE
      names(l)[le] <- paste(fMLk,"MLE",fREk,"VarExp")
      le <- length(l)+1
      l[[le]] <- ML_I_RI_WVE_MLE
      names(l)[le] <- paste(fMLk,"MLE",fREk,"VarExp Base")
    } else{
      le <- length(l)+1
      l[[le]] <- "Pendent ML_I_RI_WVE_sMLE"
      names(l)[le] <- paste("NE Model",fMLk,"MLE",fREk,"VarExp")
      le <- length(l)+1
      l[[le]] <- "Pendent ML_I_RI_WVE_MLE"
      names(l)[le] <- paste("NE Model",fMLk,"MLE",fREk,"VarExp Base")
    }
    if(WVE==2 | WVE==3){
      ML_I_RI_WVE_sREML <- summary(ML_I_RI_WVE_REML <- update(baseREML, weights = varExp(form = ~Temps)))
      le <- length(l)+1
      l[[le]] <- ML_I_RI_WVE_sREML
      names(l)[le] <- paste(fMLk,"REML",fREk,"VarExp")
      le <- length(l)+1
      l[[le]] <- ML_I_RI_WVE_REML
      names(l)[le] <- paste(fMLk,"REML",fREk,"VarExp Base")
    } else{
      le <- length(l)+1
      l[[le]] <- "Pendent ML_I_RI_WVE_sREML"
      names(l)[le] <- paste("NE Model",fMLk,"REML",fREk,"VarExp")
      le <- length(l)+1
      l[[le]] <- "Pendent ML_I_RI_WVE_REML"
      names(l)[le] <- paste("NE Model",fMLk,"REML",fREk,"VarExp Base")
    }

  
  #La funció varPower s'ha hagut de deixar sense formules ja que al fer-la dependre de la covariable temps no convergia a un resultat real en els algoritmes de nlme pels conjunts donats

    if(WVP==1 | WVP==3){
      ML_I_RI_WVP_sMLE <- summary(ML_I_RI_WVP_MLE <- update(baseML, weights = varPower()))
      le <- length(l)+1
      l[[le]] <- ML_I_RI_WVP_sMLE
      names(l)[le] <- paste(fMLk,"MLE",fREk,"VarPower")
      le <- length(l)+1
      l[[le]] <- ML_I_RI_WVP_MLE
      names(l)[le] <- paste(fMLk,"MLE",fREk,"VarPower Base")
    } else{
      le <- length(l)+1
      l[[le]] <- "Pendent ML_I_RI_WVP_sMLE"
      names(l)[le] <- paste("NE Model",fMLk,"MLE",fREk,"VarPower")
      le <- length(l)+1
      l[[le]] <- "Pendent ML_I_RI_WVP_MLE"
      names(l)[le] <- paste("NE Model",fMLk,"MLE",fREk,"VarPower Base")
    }
    if(WVP==2 | WVP==3){
      ML_I_RI_WVP_sREML <- summary(ML_I_RI_WVP_REML <- update(baseREML, weights = varPower()))
      le <- length(l)+1
      l[[le]] <- ML_I_RI_WVP_sREML
      names(l)[le] <- paste(fMLk,"REML",fREk,"VarPower")
      le <- length(l)+1
      l[[le]] <- ML_I_RI_WVP_REML
      names(l)[le] <- paste(fMLk,"REML",fREk,"VarPower Base")
    } else{
      le <- length(l)+1
      l[[le]] <- "Pendent ML_I_RI_WVP_sREML"
      names(l)[le] <- paste("NE Model",fMLk,"REML",fREk,"VarPower")
      le <- length(l)+1
      l[[le]] <- "Pendent ML_I_RI_WVP_REML"
      names(l)[le] <- paste("NE Model",fMLk,"REML",fREk,"VarPower Base")
    }

  return(l)
} 

Pel moment només s’aplica la funció sense la covariant en la fòrmula base del model per avançar d’una manera lògica en l’anàlisi de la simulació.

fML_I_He1 <- formula(He1 ~ 1)
fRE_RI <- formula(~1|Lot)
(Comp_I_RI_W_He1 <- Comp_W_LM(fML_I_He1,fRE_RI,Av_DS,3,3,3))[1:3]

## $`I MLE RI VarIdent`
## Linear mixed-effects model fit by maximum likelihood
##  Data: ds 
##       AIC     BIC   logLik
##   276.678 281.891 -135.339
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:    2.494569 5.713703
## 
## Fixed effects: list(fML) 
##                Value Std.Error DF  t-value p-value
## (Intercept) -4.31553   1.36334 36 -3.16541  0.0031
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.4204024 -0.4077357  0.1686653  0.5467960  2.0070921 
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`I MLE RI VarIdent Base`
## Linear mixed-effects model fit by maximum likelihood
##   Data: ds 
##   Log-likelihood: -135.339
##   Fixed: fML 
## (Intercept) 
##    -4.31553 
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:    2.494569 5.713703
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`I REML RI VarIdent`
## Linear mixed-effects model fit by REML
##  Data: ds 
##        AIC      BIC   logLik
##   274.1559 279.2967 -134.078
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:    2.898278 5.713711
## 
## Fixed effects: list(fML) 
##                Value Std.Error DF   t-value p-value
## (Intercept) -4.31553  1.475567 36 -2.924657  0.0059
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.3453774 -0.4085791  0.1651512  0.5794780  1.9608608 
## 
## Number of Observations: 42
## Number of Groups: 6

fML_I_He2 <- formula(He2 ~ 1)
(Comp_I_RI_W_He2 <- Comp_W_LM(fML_I_He2,fRE_RI,Av_DS,3,3,2))[1:3]

## $`I MLE RI VarIdent`
## Linear mixed-effects model fit by maximum likelihood
##  Data: ds 
##        AIC      BIC    logLik
##   362.5896 367.8027 -178.2948
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:     2.73105 16.67456
## 
## Fixed effects: list(fML) 
##                 Value Std.Error DF   t-value p-value
## (Intercept) -6.060349  2.838118 36 -2.135341  0.0396
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.87947887 -0.18837291  0.09860122  0.32935606  2.23000320 
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`I MLE RI VarIdent Base`
## Linear mixed-effects model fit by maximum likelihood
##   Data: ds 
##   Log-likelihood: -178.2948
##   Fixed: fML 
## (Intercept) 
##   -6.060349 
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:     2.73105 16.67456
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`I REML RI VarIdent`
## Linear mixed-effects model fit by REML
##  Data: ds 
##        AIC      BIC    logLik
##   358.6036 363.7443 -176.3018
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:    4.091809 16.67789
## 
## Fixed effects: list(fML) 
##                 Value Std.Error DF   t-value p-value
## (Intercept) -6.060349  3.068086 36 -1.975286  0.0559
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.7829570 -0.2167661  0.1280906  0.3463597  2.1636954 
## 
## Number of Observations: 42
## Number of Groups: 6

fML_I_He3 <- formula(He3 ~ 1)
(Comp_I_RI_W_He3 <- Comp_W_LM(fML_I_He3,fRE_RI,Av_DS,3,3,1))[1:3]

## $`I MLE RI VarIdent`
## Linear mixed-effects model fit by maximum likelihood
##  Data: ds 
##        AIC      BIC    logLik
##   596.2473 601.4603 -295.1236
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:    2.872539 272.5329
## 
## Fixed effects: list(fML) 
##                 Value Std.Error DF  t-value p-value
## (Intercept) -63.56751  42.57904 36 -1.49293  0.1442
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.7415220  0.1370678  0.2161018  0.2450529  1.6500662 
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`I MLE RI VarIdent Base`
## Linear mixed-effects model fit by maximum likelihood
##   Data: ds 
##   Log-likelihood: -295.1236
##   Fixed: fML 
## (Intercept) 
##   -63.56751 
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:    2.872539 272.5329
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`I REML RI VarIdent`
## Linear mixed-effects model fit by REML
##  Data: ds 
##        AIC      BIC    logLik
##   586.9194 592.0602 -290.4597
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:    4.271676 275.8235
## 
## Fixed effects: list(fML) 
##                 Value Std.Error DF   t-value p-value
## (Intercept) -63.56751  42.59621 36 -1.492328  0.1443
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.6964464  0.1352420  0.2134523  0.2419388  1.6300339 
## 
## Number of Observations: 42
## Number of Groups: 6

fML_IS_He1 <- formula(He1 ~ Temps)
fRE_RI <- formula(~1|Lot)
(Comp_IS_RI_W_He1 <- Comp_W_LM(fML_IS_He1,fRE_RI,Av_DS,3,3,1))[1:3]

## $`IS MLE RI VarIdent`
## Linear mixed-effects model fit by maximum likelihood
##  Data: ds 
##        AIC      BIC    logLik
##   275.3062 282.2569 -133.6531
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:     2.57674  5.45231
## 
## Fixed effects: list(fML) 
##                  Value Std.Error DF   t-value p-value
## (Intercept) -2.2383685 1.7851591 35 -1.253876  0.2182
## Temps       -0.2019462 0.1100662 35 -1.834770  0.0751
##  Correlation: 
##       (Intr)
## Temps -0.634
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.98631530 -0.57575418  0.04054785  0.52365704  2.36305907 
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`IS MLE RI VarIdent Base`
## Linear mixed-effects model fit by maximum likelihood
##   Data: ds 
##   Log-likelihood: -133.6531
##   Fixed: fML 
## (Intercept)       Temps 
##  -2.2383685  -0.2019462 
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:     2.57674  5.45231
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`IS REML RI VarIdent`
## Linear mixed-effects model fit by REML
##  Data: ds 
##        AIC    BIC    logLik
##   275.3944 282.15 -133.6972
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:     2.94882 5.529652
## 
## Fixed effects: list(fML) 
##                  Value Std.Error DF   t-value p-value
## (Intercept) -2.2383685 1.8527810 35 -1.208113  0.2351
## Temps       -0.2019462 0.1089373 35 -1.853784  0.0722
##  Correlation: 
##       (Intr)
## Temps -0.605
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.89802256 -0.56607859  0.03575378  0.52530744  2.29272070 
## 
## Number of Observations: 42
## Number of Groups: 6

fML_IS_He2 <- formula(He2 ~ Temps)
(Comp_IS_RI_W_He2 <- Comp_W_LM(fML_IS_He2,fRE_RI,Av_DS,3,3,3))[1:3]

## $`IS MLE RI VarIdent`
## Linear mixed-effects model fit by maximum likelihood
##  Data: ds 
##        AIC      BIC    logLik
##   360.5557 367.5064 -176.2779
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:    3.397521 15.76915
## 
## Fixed effects: list(fML) 
##                  Value Std.Error DF    t-value p-value
## (Intercept)  0.5415064  4.354036 35  0.1243688  0.9017
## Temps       -0.6418471  0.318333 35 -2.0162751  0.0515
##  Correlation: 
##       (Intr)
## Temps -0.752
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.42267454 -0.39506060 -0.08976259  0.45193695  2.62815324 
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`IS MLE RI VarIdent Base`
## Linear mixed-effects model fit by maximum likelihood
##   Data: ds 
##   Log-likelihood: -176.2779
##   Fixed: fML 
## (Intercept)       Temps 
##   0.5415064  -0.6418471 
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:    3.397521 15.76915
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`IS REML RI VarIdent`
## Linear mixed-effects model fit by REML
##  Data: ds 
##        AIC      BIC    logLik
##   357.0559 363.8114 -174.5279
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:    4.465284 15.99289
## 
## Fixed effects: list(fML) 
##                  Value Std.Error DF    t-value p-value
## (Intercept)  0.5415064  4.462641 35  0.1213421  0.9041
## Temps       -0.6418471  0.315069 35 -2.0371625  0.0493
##  Correlation: 
##       (Intr)
## Temps -0.726
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.31082851 -0.40557597 -0.09537236  0.39210406  2.53787466 
## 
## Number of Observations: 42
## Number of Groups: 6

fML_IS_He3 <- formula(He3 ~ Temps)
(Comp_IS_RI_W_He3 <- Comp_W_LM(fML_IS_He3,fRE_RI,Av_DS,3,3,1))[1:3]

## $`IS MLE RI VarIdent`
## Linear mixed-effects model fit by maximum likelihood
##  Data: ds 
##        AIC      BIC    logLik
##   590.8267 597.7774 -291.4134
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:    3.686905 249.4768
## 
## Fixed effects: list(fML) 
##                 Value Std.Error DF   t-value p-value
## (Intercept)  80.46670  65.12825 35  1.235511  0.2249
## Temps       -14.00333   5.03621 35 -2.780530  0.0087
##  Correlation: 
##       (Intr)
## Temps -0.795
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.31710779 -0.29797584 -0.03138667  0.28580675  2.29584125 
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`IS MLE RI VarIdent Base`
## Linear mixed-effects model fit by maximum likelihood
##   Data: ds 
##   Log-likelihood: -291.4134
##   Fixed: fML 
## (Intercept)       Temps 
##    80.46670   -14.00333 
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:    3.686905 249.4768
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`IS REML RI VarIdent`
## Linear mixed-effects model fit by REML
##  Data: ds 
##        AIC      BIC    logLik
##   576.6177 583.3732 -284.3088
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:    16.54275 255.2063
## 
## Fixed effects: list(fML) 
##                 Value Std.Error DF   t-value p-value
## (Intercept)  80.46670  65.35003 35  1.231318  0.2264
## Temps       -14.00333   5.02771 35 -2.785229  0.0086
##  Correlation: 
##       (Intr)
## Temps -0.791
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.22836471 -0.28918194 -0.02014172  0.28163899  2.23407518 
## 
## Number of Observations: 42
## Number of Groups: 6

fML_IS_He1 <- formula(He1 ~ Temps)
fRE_RIS <- formula(~1 + Temps|Lot)
(Comp_IS_RIS_W_He1 <- Comp_W_LM(fML_IS_He1,fRE_RIS,Av_DS,3,3,0))[1:3]

## $`IS MLE RIS VarIdent`
## Linear mixed-effects model fit by maximum likelihood
##  Data: ds 
##        AIC      BIC   logLik
##   275.5219 285.9479 -131.761
## 
## Random effects:
##  Formula: ~1 + Temps | Lot
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev    Corr  
## (Intercept) 0.3575705 (Intr)
## Temps       0.2451125 0.403 
## Residual    5.0469114       
## 
## Fixed effects: list(fML) 
##                  Value Std.Error DF   t-value p-value
## (Intercept) -2.2383685 1.3256403 35 -1.688519  0.1002
## Temps       -0.2019462 0.1445477 35 -1.397090  0.1712
##  Correlation: 
##       (Intr)
## Temps -0.525
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.82242826 -0.73863965  0.06046065  0.59009730  2.20606476 
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`IS MLE RIS VarIdent Base`
## Linear mixed-effects model fit by maximum likelihood
##   Data: ds 
##   Log-likelihood: -131.761
##   Fixed: fML 
## (Intercept)       Temps 
##  -2.2383685  -0.2019462 
## 
## Random effects:
##  Formula: ~1 + Temps | Lot
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev    Corr  
## (Intercept) 0.3575705 (Intr)
## Temps       0.2451125 0.403 
## Residual    5.0469114       
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`IS REML RIS VarIdent`
## Linear mixed-effects model fit by REML
##  Data: ds 
##       AIC      BIC   logLik
##   275.476 285.6092 -131.738
## 
## Random effects:
##  Formula: ~1 + Temps | Lot
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev    Corr  
## (Intercept) 0.4223037 (Intr)
## Temps       0.2747036 0.381 
## Residual    5.1178359       
## 
## Fixed effects: list(fML) 
##                  Value Std.Error DF  t-value p-value
## (Intercept) -2.2383685 1.3148464 35 -1.70238  0.0976
## Temps       -0.2019462 0.1508063 35 -1.33911  0.1892
##  Correlation: 
##       (Intr)
## Temps -0.49 
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.81310975 -0.72359743  0.04793833  0.57533392  2.16363940 
## 
## Number of Observations: 42
## Number of Groups: 6

fML_IS_He2 <- formula(He2 ~ Temps)
(Comp_IS_RIS_W_He2 <- Comp_W_LM(fML_IS_He2,fRE_RIS,Av_DS,3,3,0))[1:3]

## $`IS MLE RIS VarIdent`
## Linear mixed-effects model fit by maximum likelihood
##  Data: ds 
##        AIC      BIC    logLik
##   360.0596 370.4856 -174.0298
## 
## Random effects:
##  Formula: ~1 + Temps | Lot
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev     Corr  
## (Intercept)  1.4623857 (Intr)
## Temps        0.6883651 -0.781
## Residual    13.9851505       
## 
## Fixed effects: list(fML) 
##                  Value Std.Error DF    t-value p-value
## (Intercept)  0.5415064  3.700844 35  0.1463197  0.8845
## Temps       -0.6418471  0.403271 35 -1.5916030  0.1205
##  Correlation: 
##       (Intr)
## Temps -0.642
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.12232367 -0.50804753 -0.09390181  0.40876992  2.59165264 
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`IS MLE RIS VarIdent Base`
## Linear mixed-effects model fit by maximum likelihood
##   Data: ds 
##   Log-likelihood: -174.0298
##   Fixed: fML 
## (Intercept)       Temps 
##   0.5415064  -0.6418471 
## 
## Random effects:
##  Formula: ~1 + Temps | Lot
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev     Corr  
## (Intercept)  1.4623857 (Intr)
## Temps        0.6883651 -0.781
## Residual    13.9851505       
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`IS REML RIS VarIdent`
## Linear mixed-effects model fit by REML
##  Data: ds 
##        AIC      BIC    logLik
##   356.1293 366.2626 -172.0647
## 
## Random effects:
##  Formula: ~1 + Temps | Lot
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev     Corr  
## (Intercept)  1.6296805 (Intr)
## Temps        0.7737087 -0.77 
## Residual    14.1868088       
## 
## Fixed effects: list(fML) 
##                  Value Std.Error DF   t-value p-value
## (Intercept)  0.5415064  3.674070 35  0.147386  0.8837
## Temps       -0.6418471  0.421764 35 -1.521817  0.1370
##  Correlation: 
##       (Intr)
## Temps -0.623
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.11730124 -0.50782927 -0.09011089  0.39414295  2.53608344 
## 
## Number of Observations: 42
## Number of Groups: 6

fML_IS_He3 <- formula(He3 ~ Temps)
(Comp_IS_RIS_W_He3 <- Comp_W_LM(fML_IS_He3,fRE_RIS,Av_DS,3,3,3))[1:3]

## Warning in logLik.reStruct(object, conLin): Singular precision matrix in
## level -1, block 6

## Warning in logLik.reStruct(object, conLin): Singular precision matrix in
## level -1, block 1

## Warning in logLik.reStruct(object, conLin): Singular precision matrix in
## level -1, block 1

## Warning in logLik.reStruct(object, conLin): Singular precision matrix in
## level -1, block 1

## $`IS MLE RIS VarIdent`
## Linear mixed-effects model fit by maximum likelihood
##  Data: ds 
##        AIC      BIC    logLik
##   587.8952 598.3212 -287.9476
## 
## Random effects:
##  Formula: ~1 + Temps | Lot
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev    Corr  
## (Intercept)  74.32588 (Intr)
## Temps        14.66258 -0.97 
## Residual    206.88256       
## 
## Fixed effects: list(fML) 
##                 Value Std.Error DF   t-value p-value
## (Intercept)  80.46670  62.30614 35  1.291473  0.2050
## Temps       -14.00333   7.42061 35 -1.887086  0.0675
##  Correlation: 
##       (Intr)
## Temps -0.788
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.38691287 -0.31561560 -0.04199045  0.49771014  2.70225065 
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`IS MLE RIS VarIdent Base`
## Linear mixed-effects model fit by maximum likelihood
##   Data: ds 
##   Log-likelihood: -287.9476
##   Fixed: fML 
## (Intercept)       Temps 
##    80.46670   -14.00333 
## 
## Random effects:
##  Formula: ~1 + Temps | Lot
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev    Corr  
## (Intercept)  74.32588 (Intr)
## Temps        14.66258 -0.97 
## Residual    206.88256       
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`IS REML RIS VarIdent`
## Linear mixed-effects model fit by REML
##  Data: ds 
##       AIC      BIC   logLik
##   572.962 583.0952 -280.481
## 
## Random effects:
##  Formula: ~1 + Temps | Lot
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev    Corr  
## (Intercept)  82.99108 (Intr)
## Temps        16.41286 -0.969
## Residual    209.79908       
## 
## Fixed effects: list(fML) 
##                 Value Std.Error DF   t-value p-value
## (Intercept)  80.46670  63.27106 35  1.271777  0.2118
## Temps       -14.00333   7.87274 35 -1.778711  0.0840
##  Correlation: 
##       (Intr)
## Temps -0.794
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.35238381 -0.30889313 -0.03774284  0.47470858  2.66513171 
## 
## Number of Observations: 42
## Number of Groups: 6

RESUM I DIAGNÒSTIC MODELS He

En el codi resum que es mostra a continuació s’han realitzat varies generaciones de llistats i de gràfics segons els resultats que s’han anat obtenint tenint en compte els tres indicadors de qualitat. Es comenten al final ja que les conclusions són similars pels 3 conjunts.

#He1 sense funció de variància
He1_R0 <- list(Comp_I_IS_He1[4:5],Comp_I_RI_He1[1:3],Comp_IS_RI_He1[1:3])
He1_R <- ConvList(He1_R0)
ML_Resum(He1_R,"He1")

## [[1]]
##                   Sigma      AIC      BIC
## He1--OLS I base  6.3101 276.9192 280.3946
## He1--OLS IS base 6.1794 276.1249 281.3379
## He1--OLS I RI    5.7137 229.6502 231.2337
## He1--MLE I RI    5.7137 276.6780 281.8910
## He1--REML I RI   5.7137 274.1559 279.2967
## He1--OLS IS RI   5.5297 228.2784 231.4455
## He1--MLE IS RI   5.4523 275.3062 282.2569
## He1--REML IS RI  5.5296 275.3944 282.1500
## 
## [[2]]

#He1 funcions var1
He1_RWa0 <- list(Comp_I_RI_W_He1[c(1,3,5,7,9,11)])
He1_RWa <- ConvList(He1_RWa0)
ML_Resum(He1_RWa,"He1")

## [[1]]
##                          Sigma      AIC      BIC
## He1--I MLE RI VarIdent  5.7137 276.6780 281.8910
## He1--I REML RI VarIdent 5.7137 274.1559 279.2967
## He1--I MLE RI VarExp    2.0089 258.3204 265.2711
## He1--I REML RI VarExp   2.0215 256.5681 263.4224
## He1--I MLE RI VarPower  4.0571 277.7818 284.7325
## He1--I REML RI VarPower 4.3576 275.3532 282.2075
## 
## [[2]]

#He1 funcions var2
He1_RWb0 <- list(Comp_IS_RI_W_He1[c(1,3,5,7,9)])
He1_RWb <- ConvList(He1_RWb0)
ML_Resum(He1_RWb,"He1")

## [[1]]
##                           Sigma      AIC      BIC
## He1--IS MLE RI VarIdent  5.4523 275.3062 282.2569
## He1--IS REML RI VarIdent 5.5297 275.3944 282.1500
## He1--IS MLE RI VarExp    2.0204 259.7077 268.3960
## He1--IS REML RI VarExp   2.0582 260.5776 269.0220
## He1--IS MLE RI VarPower  0.1512 263.0674 271.7558
## 
## [[2]]

#He1 funcions var3
He1_RWc0 <- list(Comp_IS_RIS_W_He1[c(1,3,5,7)])
He1_RWc <- ConvList(He1_RWc0)
ML_Resum(He1_RWc,"He1")

## [[1]]
##                            Sigma      AIC      BIC
## He1--IS MLE RIS VarIdent  5.0469 275.5219 285.9479
## He1--IS REML RIS VarIdent 5.1178 275.4760 285.6092
## He1--IS MLE RIS VarExp    2.0843 263.2686 275.4323
## He1--IS REML RIS VarExp   2.1276 264.0129 275.8350
## 
## [[2]]

He1_RGW0 <- list(He1_R[c(3,6)],He1_RWa[c(4,6)],He1_RWb[c(4,5)],He1_RWc[c(3,4)])
He1_RGW <- ConvList(He1_RGW0)
ML_Resum(He1_RGW,"He1")

## [[1]]
##                          Sigma      AIC      BIC
## He1--OLS I RI           5.7137 229.6502 231.2337
## He1--OLS IS RI          5.5297 228.2784 231.4455
## He1--I REML RI VarExp   2.0215 256.5681 263.4224
## He1--I REML RI VarPower 4.3576 275.3532 282.2075
## He1--IS REML RI VarExp  2.0582 260.5776 269.0220
## He1--IS MLE RI VarPower 0.1512 263.0674 271.7558
## He1--IS MLE RIS VarExp  2.0843 263.2686 275.4323
## He1--IS REML RIS VarExp 2.1276 264.0129 275.8350
## 
## [[2]]

#He2 sense funció de variància
He2_R0 <- list(Comp_I_IS_He2[4:5],Comp_I_RI_He2[1:3],Comp_IS_RI_He2[1:3])
He2_R <- ConvList(He2_R0)
ML_Resum(He2_R,"He2")

## [[1]]
##                    Sigma      AIC      BIC
## He2--OLS I base  17.1010 360.6663 364.1417
## He2--OLS IS base 16.5293 358.7732 363.9862
## He2--OLS I RI    16.6779 306.7776 308.3612
## He2--MLE I RI    16.6779 362.5896 367.8026
## He2--REML I RI   16.6779 358.6036 363.7443
## He2--OLS IS RI   15.9929 304.7437 307.9108
## He2--MLE IS RI   15.7692 360.5557 367.5064
## He2--REML IS RI  15.9929 357.0559 363.8114
## 
## [[2]]

#He2 funcions var
He2_RWa0 <- list(Comp_I_RI_W_He2[c(1,3,5,7,11)])
He2_RWa <- ConvList(He2_RWa0)
ML_Resum(He2_RWa,"He2")

## [[1]]
##                           Sigma      AIC      BIC
## He2--I MLE RI VarIdent  16.6746 362.5896 367.8027
## He2--I REML RI VarIdent 16.6779 358.6036 363.7443
## He2--I MLE RI VarExp     2.1259 309.9968 316.9474
## He2--I REML RI VarExp    2.1342 308.2841 315.1383
## He2--I REML RI VarPower 10.3850 359.9306 366.7849
## 
## [[2]]

#He2 funcions var
He2_RWb0 <- list(Comp_IS_RI_W_He2[c(1,3,5,7,9,11)])
He2_RWb <- ConvList(He2_RWb0)
ML_Resum(He2_RWb,"He2")

## [[1]]
##                            Sigma      AIC      BIC
## He2--IS MLE RI VarIdent  15.7691 360.5557 367.5064
## He2--IS REML RI VarIdent 15.9929 357.0559 363.8114
## He2--IS MLE RI VarExp     2.1252 311.9960 320.6844
## He2--IS REML RI VarExp    2.1970 311.6457 320.0901
## He2--IS MLE RI VarPower   0.0000 312.6716 321.3600
## He2--IS REML RI VarPower  0.0000 312.5798 321.0242
## 
## [[2]]

He2_RWbb0 <- list(Comp_IS_RI_W_He2[c(9,11)])
He2_RWbb <- ConvList(He2_RWbb0)
ML_Resum(He2_RWbb,"He2")

## [[1]]
##                          Sigma      AIC      BIC
## He2--IS MLE RI VarPower      0 312.6716 321.3600
## He2--IS REML RI VarPower     0 312.5798 321.0242
## 
## [[2]]

## Warning: Removed 2 rows containing missing values (geom_bar).

#He2 funcions var3
He2_RWc0 <- list(Comp_IS_RIS_W_He2[c(1,3,5,7)])
He2_RWc <- ConvList(He2_RWc0)
ML_Resum(He2_RWc,"He2")

## [[1]]
##                             Sigma      AIC      BIC
## He2--IS MLE RIS VarIdent  13.9852 360.0596 370.4856
## He2--IS REML RIS VarIdent 14.1868 356.1293 366.2626
## He2--IS MLE RIS VarExp     2.1253 315.9962 328.1599
## He2--IS REML RIS VarExp    2.5671 316.2920 328.1141
## 
## [[2]]

He2_RGW0 <- list(He2_R[c(3,6)],He2_RWa[c(4,5)],He2_RWb[c(4,6)],He2_RWc[c(3,4)])
He2_RGW <- ConvList(He2_RGW0)
ML_Resum(He2_RGW,"He2")

## [[1]]
##                            Sigma      AIC      BIC
## He2--OLS I RI            16.6779 306.7776 308.3612
## He2--OLS IS RI           15.9929 304.7437 307.9108
## He2--I REML RI VarExp     2.1342 308.2841 315.1383
## He2--I REML RI VarPower  10.3850 359.9306 366.7849
## He2--IS REML RI VarExp    2.1970 311.6457 320.0901
## He2--IS REML RI VarPower  0.0000 312.5798 321.0242
## He2--IS MLE RIS VarExp    2.1253 315.9962 328.1599
## He2--IS REML RIS VarExp   2.5671 316.2920 328.1141
## 
## [[2]]

#He3 sense funció de variància
He3_R0 <- list(Comp_I_IS_He3[4:5],Comp_I_RI_He3[1:3],Comp_IS_RI_He3[1:3])
He3_R <- ConvList(He3_R0)
ML_Resum(He3_R,"He3")

## [[1]]
##                     Sigma      AIC      BIC
## He3--OLS I base  275.8477 594.2460 597.7214
## He3--OLS IS base 255.6619 588.8256 594.0386
## He3--OLS I RI    278.5817 509.5028 511.0863
## He3--MLE I RI    272.5440 596.2460 601.4590
## He3--REML I RI   275.8477 586.9182 592.0589
## He3--OLS IS RI   255.7137 504.3216 507.4887
## He3--MLE IS RI   249.5005 590.8256 597.7762
## He3--REML IS RI  255.6619 576.6155 583.3710
## 
## [[2]]

#He3 funcions var
He3_RWa0 <- list(Comp_I_RI_W_He3[c(1,3,5,7,9)])
He3_RWa <- ConvList(He3_RWa0)
ML_Resum(He3_RWa,"He3")

## [[1]]
##                            Sigma      AIC      BIC
## He3--I MLE RI VarIdent  272.5329 596.2473 601.4603
## He3--I REML RI VarIdent 275.8235 586.9194 592.0602
## He3--I MLE RI VarExp      2.6404 429.2263 436.1770
## He3--I REML RI VarExp     2.7460 427.4230 434.2773
## He3--I MLE RI VarPower    0.5964 598.2460 605.1967
## 
## [[2]]

He3_RWaa0 <- list(Comp_I_RI_W_He3[c(5,7,9)])
He3_RWaa <- ConvList(He3_RWaa0)
ML_Resum(He3_RWaa,"He3")

## [[1]]
##                         Sigma      AIC      BIC
## He3--I MLE RI VarExp   2.6404 429.2263 436.1770
## He3--I REML RI VarExp  2.7460 427.4230 434.2773
## He3--I MLE RI VarPower 0.5964 598.2460 605.1967
## 
## [[2]]

#He3 funcions var b
He3_RWb0 <- list(Comp_IS_RI_W_He3[c(1,3,5,7,9)])
He3_RWb <- ConvList(He3_RWb0)
ML_Resum(He3_RWb,"He3")

## [[1]]
##                             Sigma      AIC      BIC
## He3--IS MLE RI VarIdent  249.4768 590.8267 597.7774
## He3--IS REML RI VarIdent 255.2063 576.6177 583.3732
## He3--IS MLE RI VarExp      2.6227 431.1668 439.8552
## He3--IS REML RI VarExp     2.8108 428.6916 437.1360
## He3--IS MLE RI VarPower    0.0103 429.3102 437.9985
## 
## [[2]]

He3_RWbb0 <- list(Comp_IS_RI_W_He3[c(5,7,9)])
He3_RWbb <- ConvList(He3_RWbb0)
ML_Resum(He3_RWbb,"He3")

## [[1]]
##                          Sigma      AIC      BIC
## He3--IS MLE RI VarExp   2.6227 431.1668 439.8552
## He3--IS REML RI VarExp  2.8108 428.6916 437.1360
## He3--IS MLE RI VarPower 0.0103 429.3102 437.9985
## 
## [[2]]

#He3 funcions var3
He3_RWc0 <- list(Comp_IS_RIS_W_He3[c(1,3,5,7,9,11)])
He3_RWc <- ConvList(He3_RWc0)
ML_Resum(He3_RWc,"He3")

## [[1]]
##                              Sigma      AIC      BIC
## He3--IS MLE RIS VarIdent  206.8826 587.8952 598.3212
## He3--IS REML RIS VarIdent 209.7991 572.9620 583.0952
## He3--IS MLE RIS VarExp      2.6347 435.1639 447.3276
## He3--IS REML RIS VarExp     2.8263 432.6939 444.5161
## He3--IS MLE RIS VarPower    0.0103 433.3102 445.4739
## He3--IS REML RIS VarPower   0.0112 431.3242 443.1463
## 
## [[2]]

He3_RWcc0 <- list(Comp_IS_RIS_W_He3[c(5,7,9,11)])
He3_RWcc <- ConvList(He3_RWcc0)
ML_Resum(He3_RWcc,"He3")

## [[1]]
##                            Sigma      AIC      BIC
## He3--IS MLE RIS VarExp    2.6347 435.1639 447.3276
## He3--IS REML RIS VarExp   2.8263 432.6939 444.5161
## He3--IS MLE RIS VarPower  0.0103 433.3102 445.4739
## He3--IS REML RIS VarPower 0.0112 431.3242 443.1463
## 
## [[2]]

He3_RGW0 <- list(He3_R[c(3,6)],He3_RWa[c(4,5)],He3_RWb[c(4,5)],He3_RWc[c(4,6)])
He3_RGW <- ConvList(He3_RGW0)
ML_Resum(He3_RGW,"He3")

## [[1]]
##                              Sigma      AIC      BIC
## He3--OLS I RI             278.5817 509.5028 511.0863
## He3--OLS IS RI            255.7137 504.3216 507.4887
## He3--I REML RI VarExp       2.7460 427.4230 434.2773
## He3--I MLE RI VarPower      0.5964 598.2460 605.1967
## He3--IS REML RI VarExp      2.8108 428.6916 437.1360
## He3--IS MLE RI VarPower     0.0103 429.3102 437.9985
## He3--IS REML RIS VarExp     2.8263 432.6939 444.5161
## He3--IS REML RIS VarPower   0.0112 431.3242 443.1463
## 
## [[2]]

He3_RGWb0 <- list(He3_RWa[c(4,5)],He3_RWb[c(4,5)],He3_RWc[c(4,6)])
He3_RGWb <- ConvList(He3_RGWb0)
ML_Resum(He3_RGWb,"He3")

## [[1]]
##                            Sigma      AIC      BIC
## He3--I REML RI VarExp     2.7460 427.4230 434.2773
## He3--I MLE RI VarPower    0.5964 598.2460 605.1967
## He3--IS REML RI VarExp    2.8108 428.6916 437.1360
## He3--IS MLE RI VarPower   0.0103 429.3102 437.9985
## He3--IS REML RIS VarExp   2.8263 432.6939 444.5161
## He3--IS REML RIS VarPower 0.0112 431.3242 443.1463
## 
## [[2]]

Els models sense la funció de l’error estàndard i amb efectes aleatoris de l’intercepció són similars entre ells. Al afegir al model la funció de modificació de l’error estàndard s’observa que es redueix significativament l’error estàndard residual, sobretot quan s’apliquen la funció varPower juntament amb el model fixe que conté la covariable temps que sembla el model on es redueix més l’error estàndard residual (concretament als models amb la funció varPower es nota la diferència al afegir la variable temps a la part d’efectes fixos en relació al model on només compté la intercepció com efecte fix).

Pel que fa als indicadors AIC/BIC no hi ha realment una millora substancial al anar augmentant la complexitat del model encara que sí que denota una lleugera reducció (i per tant millora en la bondat d’ajust) en general al afegir les funcions de modulació de variància.

Segons el conjunt no s’ha pogut aplicar la varPower a totes les variants del model ja que en algun cas s’arribava a resultats inconsistents (cal recordar que al aplicar algoritmes d’aquesta complexitat matemàtica segons com es distribueixen les nostres dades es pot no arribar a una solució matemàticament lògica), tot i així no es denota una millora significativa al modificar els efectes aleatoris de la intercepció per efectes aleatoris de pendent i intercepció en aquests casos, que ja quadraria amb els conjunts simulats. S’entén que la importància d’afegir la pendent en els efectes fixos ve com en el conjunt Ho2 pel fet que la pròpia dispersió al haver un nombre reduït de lots es falseja en forma de pendent segons com s’hagi produit la dispersió.

En aquest tipus de conjunts de dades es pot extreure que es considera important aplicar la funció de modular l’error estàndard i en cas necessari per compensar la mida dels lots aplicar el factor de la pendent, tot i que en aquest especte s’ha d’anar amb compte si es vol utilitzar com un model de predicció per futurs lots.

Prenent el model més efectiu com el Model IS MLE/REML RI VarPower s’extreuen les matrius internes d’interès per analitzar-lo amb més detall:

He_B <- list(Comp_IS_RI_W_He1[[10]],Comp_IS_RI_W_He2[[10]],Comp_IS_RI_W_He3[[10]])
names(He_B) <- c(names(Comp_IS_RI_W_He1)[10],names(Comp_IS_RI_W_He2)[10],names(Comp_IS_RI_W_He3)[10])

lapply(1:length(He_B),function(B){
  VarCorr(He_B[[B]])
})

## [[1]]
## Lot = pdLogChol(1) 
##             Variance     StdDev      
## (Intercept) 1.522746e-07 0.0003902238
## Residual    2.286472e-02 0.1512108315
## 
## [[2]]
## Lot = pdLogChol(1) 
##             Variance     StdDev      
## (Intercept) 8.002497e-49 8.945668e-25
## Residual    1.592017e-43 3.990008e-22
## 
## [[3]]
## Lot = pdLogChol(1) 
##             Variance     StdDev      
## (Intercept) 4.884495e-09 6.988916e-05
## Residual    1.063421e-04 1.031223e-02

lapply(1:length(He_B),function(B){
  getVarCov(He_B[[B]],type="conditional")
})

## [[1]]
## Lot Lot 1 
## Conditional variance covariance matrix
##        1      2      3      4      5      6      7
## 1 5.7349 0.0000  0.000  0.000  0.000  0.000   0.00
## 2 0.0000 9.4577  0.000  0.000  0.000  0.000   0.00
## 3 0.0000 0.0000 14.904  0.000  0.000  0.000   0.00
## 4 0.0000 0.0000  0.000 22.614  0.000  0.000   0.00
## 5 0.0000 0.0000  0.000  0.000 33.228  0.000   0.00
## 6 0.0000 0.0000  0.000  0.000  0.000 66.311   0.00
## 7 0.0000 0.0000  0.000  0.000  0.000  0.000 121.68
##   Standard Deviations: 2.3948 3.0753 3.8606 4.7554 5.7644 8.1431 11.031 
## 
## [[2]]
## Lot Lot 1 
## Conditional variance covariance matrix
##        1      2     3      4      5      6      7
## 1 5.9403  0.000  0.00  0.000   0.00   0.00    0.0
## 2 0.0000 12.833  0.00  0.000   0.00   0.00    0.0
## 3 0.0000  0.000 27.52  0.000   0.00   0.00    0.0
## 4 0.0000  0.000  0.00 58.593   0.00   0.00    0.0
## 5 0.0000  0.000  0.00  0.000 123.87   0.00    0.0
## 6 0.0000  0.000  0.00  0.000   0.00 542.33    0.0
## 7 0.0000  0.000  0.00  0.000   0.00   0.00 2311.7
##   Standard Deviations: 2.4373 3.5823 5.2459 7.6546 11.13 23.288 48.08 
## 
## [[3]]
## Lot Lot 1 
## Conditional variance covariance matrix
##        1      2      3      4      5     6      7
## 1 2.9705  0.000   0.00    0.0    0.0     0      0
## 2 0.0000 33.714   0.00    0.0    0.0     0      0
## 3 0.0000  0.000 245.19    0.0    0.0     0      0
## 4 0.0000  0.000   0.00 1311.9    0.0     0      0
## 5 0.0000  0.000   0.00    0.0 5607.2     0      0
## 6 0.0000  0.000   0.00    0.0    0.0 63498      0
## 7 0.0000  0.000   0.00    0.0    0.0     0 461120
##   Standard Deviations: 1.7235 5.8064 15.659 36.22 74.881 251.99 679.06

lapply(1:length(He_B),function(B){
  a<-getVarCov(He_B[[B]],type="marginal")
  b<-cov2cor(a[[1]])
  return(list(a,b))
})

## [[1]]
## [[1]][[1]]
## Lot Lot 1 
## Marginal variance covariance matrix
##            1          2          3          4          5          6
## 1 5.7349e+00 1.5227e-07 1.5227e-07 1.5227e-07 1.5227e-07 1.5227e-07
## 2 1.5227e-07 9.4577e+00 1.5227e-07 1.5227e-07 1.5227e-07 1.5227e-07
## 3 1.5227e-07 1.5227e-07 1.4904e+01 1.5227e-07 1.5227e-07 1.5227e-07
## 4 1.5227e-07 1.5227e-07 1.5227e-07 2.2614e+01 1.5227e-07 1.5227e-07
## 5 1.5227e-07 1.5227e-07 1.5227e-07 1.5227e-07 3.3228e+01 1.5227e-07
## 6 1.5227e-07 1.5227e-07 1.5227e-07 1.5227e-07 1.5227e-07 6.6311e+01
## 7 1.5227e-07 1.5227e-07 1.5227e-07 1.5227e-07 1.5227e-07 1.5227e-07
##            7
## 1 1.5227e-07
## 2 1.5227e-07
## 3 1.5227e-07
## 4 1.5227e-07
## 5 1.5227e-07
## 6 1.5227e-07
## 7 1.2168e+02
##   Standard Deviations: 2.3948 3.0753 3.8606 4.7554 5.7644 8.1431 11.031 
## 
## [[1]][[2]]
##              1            2            3            4            5
## 1 1.000000e+00 2.067628e-08 1.647076e-08 1.337155e-08 1.103090e-08
## 2 2.067628e-08 1.000000e+00 1.282573e-08 1.041238e-08 8.589726e-09
## 3 1.647076e-08 1.282573e-08 1.000000e+00 8.294522e-09 6.842590e-09
## 4 1.337155e-08 1.041238e-08 8.294522e-09 1.000000e+00 5.555059e-09
## 5 1.103090e-08 8.589726e-09 6.842590e-09 5.555059e-09 1.000000e+00
## 6 7.808606e-09 6.080536e-09 4.843765e-09 3.932341e-09 3.243997e-09
## 7 5.764336e-09 4.488669e-09 3.575681e-09 2.902866e-09 2.394728e-09
##              6            7
## 1 7.808606e-09 5.764336e-09
## 2 6.080536e-09 4.488669e-09
## 3 4.843765e-09 3.575681e-09
## 4 3.932341e-09 2.902866e-09
## 5 3.243997e-09 2.394728e-09
## 6 1.000000e+00 1.695191e-09
## 7 1.695191e-09 1.000000e+00
## 
## 
## [[2]]
## [[2]][[1]]
## Lot Lot 1 
## Marginal variance covariance matrix
##            1          2          3          4          5          6
## 1 5.9403e+00 8.0025e-49 8.0025e-49 8.0025e-49 8.0025e-49 8.0025e-49
## 2 8.0025e-49 1.2833e+01 8.0025e-49 8.0025e-49 8.0025e-49 8.0025e-49
## 3 8.0025e-49 8.0025e-49 2.7520e+01 8.0025e-49 8.0025e-49 8.0025e-49
## 4 8.0025e-49 8.0025e-49 8.0025e-49 5.8593e+01 8.0025e-49 8.0025e-49
## 5 8.0025e-49 8.0025e-49 8.0025e-49 8.0025e-49 1.2387e+02 8.0025e-49
## 6 8.0025e-49 8.0025e-49 8.0025e-49 8.0025e-49 8.0025e-49 5.4233e+02
## 7 8.0025e-49 8.0025e-49 8.0025e-49 8.0025e-49 8.0025e-49 8.0025e-49
##            7
## 1 8.0025e-49
## 2 8.0025e-49
## 3 8.0025e-49
## 4 8.0025e-49
## 5 8.0025e-49
## 6 8.0025e-49
## 7 2.3117e+03
##   Standard Deviations: 2.4373 3.5823 5.2459 7.6546 11.13 23.288 48.08 
## 
## [[2]][[2]]
##              1            2            3            4            5
## 1 1.000000e+00 9.165659e-50 6.258933e-50 4.289438e-50 2.950086e-50
## 2 9.165659e-50 1.000000e+00 4.258369e-50 2.918390e-50 2.007140e-50
## 3 6.258933e-50 4.258369e-50 1.000000e+00 1.992874e-50 1.370611e-50
## 4 4.289438e-50 2.918390e-50 1.992874e-50 1.000000e+00 9.393217e-51
## 5 2.950086e-50 2.007140e-50 1.370611e-50 9.393217e-51 1.000000e+00
## 6 1.409916e-50 9.592599e-51 6.550476e-51 4.489241e-51 3.087502e-51
## 7 6.829012e-51 4.646231e-51 3.172761e-51 2.174390e-51 1.495450e-51
##              6            7
## 1 1.409916e-50 6.829012e-51
## 2 9.592599e-51 4.646231e-51
## 3 6.550476e-51 3.172761e-51
## 4 4.489241e-51 2.174390e-51
## 5 3.087502e-51 1.495450e-51
## 6 1.000000e+00 7.147109e-52
## 7 7.147109e-52 1.000000e+00
## 
## 
## [[3]]
## [[3]][[1]]
## Lot Lot 1 
## Marginal variance covariance matrix
##            1          2          3          4          5          6
## 1 2.9705e+00 4.8845e-09 4.8845e-09 4.8845e-09 4.8845e-09 4.8845e-09
## 2 4.8845e-09 3.3714e+01 4.8845e-09 4.8845e-09 4.8845e-09 4.8845e-09
## 3 4.8845e-09 4.8845e-09 2.4519e+02 4.8845e-09 4.8845e-09 4.8845e-09
## 4 4.8845e-09 4.8845e-09 4.8845e-09 1.3119e+03 4.8845e-09 4.8845e-09
## 5 4.8845e-09 4.8845e-09 4.8845e-09 4.8845e-09 5.6072e+03 4.8845e-09
## 6 4.8845e-09 4.8845e-09 4.8845e-09 4.8845e-09 4.8845e-09 6.3498e+04
## 7 4.8845e-09 4.8845e-09 4.8845e-09 4.8845e-09 4.8845e-09 4.8845e-09
##            7
## 1 4.8845e-09
## 2 4.8845e-09
## 3 4.8845e-09
## 4 4.8845e-09
## 5 4.8845e-09
## 6 4.8845e-09
## 7 4.6112e+05
##   Standard Deviations: 1.7235 5.8064 15.659 36.22 74.881 251.99 679.06 
## 
## [[3]][[2]]
##              1            2            3            4            5
## 1 1.000000e+00 4.880886e-10 1.809884e-10 7.824538e-11 3.784728e-11
## 2 4.880886e-10 1.000000e+00 5.372265e-11 2.322551e-11 1.123418e-11
## 3 1.809884e-10 5.372265e-11 1.000000e+00 8.612266e-12 4.165752e-12
## 4 7.824538e-11 2.322551e-11 8.612266e-12 1.000000e+00 1.800949e-12
## 5 3.784728e-11 1.123418e-11 4.165752e-12 1.800949e-12 1.000000e+00
## 6 1.124670e-11 3.338348e-12 1.237895e-12 5.351699e-13 2.588616e-13
## 7 4.173490e-12 1.238814e-12 4.593652e-13 1.985939e-13 9.605987e-14
##              6            7
## 1 1.124670e-11 4.173490e-12
## 2 3.338348e-12 1.238814e-12
## 3 1.237895e-12 4.593652e-13
## 4 5.351699e-13 1.985939e-13
## 5 2.588616e-13 9.605987e-14
## 6 1.000000e+00 2.854515e-14
## 7 2.854515e-14 1.000000e+00

Tot i obtenir resultats diferents per cada conjunt He les conclusions són semblants:

• En els primers resultats s’observen a les matrius D donant poca informació, simplement es veu la baixa escala obtinguda per l’variància residual.
• La segona sèria de matrius són les matrius R que reflexen la matriu de variàncies-covariàncies residuals dels diferents temps que tenen tots els lots. S’està en el cas on es considera l’variància residual heterogènia i per tant es veu com la funció addicional recalcula de manera que a mesura que s’avança en els punts de temps augmenta l’variància pròpia del temps (i si el conjunt és més dispers també més variància es calcula).
• La tercera sèria de matrius mostra:
  • La matriu primera de cada punt de la sèria és de tipus marginal on es suma l’variància residual i l’variància degut a l’efecte aleatori en la diagonal, mentre que l’variància degut a l’efecte aleatori forma part de la covariància en aquest cas. Això és degut a que en aquesta matriu es reflexa l’efecte de l’variància total dins un dels individus (equivalent per tots els altres) que té l’variància de l’efecte aleatori sumat a l’variància residual al centrar-nos en un punt, i es relaciona amb la resta de punts corresponents als altres temps amb el factor calculat aleatori de l’variància dins el mateix lot.
  • Es reflecteix també a la segona matriu on es plasma la correlació de variàncies i on es veu que s’ha calculat una correlació variant segons quin temps es relaciona amb quin, però baixes en general (és a dir que per aquest cas no s’ha aconseguit trobar un model que expliqui bé aquesta correlació, o simplement al haver una dispersió en funció del temps baixa la correlació).

Es visualitzen els gràfics de diagnòstic del model:

He_Diag0 <- list(Comp_I_IS_He1[5],Comp_I_IS_He2[5],Comp_I_IS_He3[5],
                Comp_IS_RI_W_He1[10],Comp_IS_RI_W_He2[10],Comp_IS_RI_W_He3[10])

He_Diag <- ConvList(He_Diag0)

names(He_Diag) <- paste(HeVars,names(He_Diag))

lapply(1:length(He_Diag),function(QQ){
  if (class(He_Diag[[QQ]])=="lme") {
    list(qqnorm(He_Diag[[QQ]],abline = c(0,1),main=names(He_Diag)[QQ]),plot(He_Diag[[QQ]]))
  } else {
    list(plot(He_Diag[[QQ]],which=1,main=names(He_Diag)[QQ]),plot(He_Diag[[QQ]],which=2,main=names(He_Diag)[QQ]))
  }
  })

## Warning in if (class(He_Diag[[QQ]]) == "lme") {: the condition has length >
## 1 and only the first element will be used

## Warning in if (class(He_Diag[[QQ]]) == "lme") {: the condition has length >
## 1 and only the first element will be used

## Warning in if (class(He_Diag[[QQ]]) == "lme") {: the condition has length >
## 1 and only the first element will be used

## [[1]]
## [[1]][[1]]
## NULL
## 
## [[1]][[2]]
## NULL
## 
## 
## [[2]]
## [[2]][[1]]
## NULL
## 
## [[2]][[2]]
## NULL
## 
## 
## [[3]]
## [[3]][[1]]
## NULL
## 
## [[3]][[2]]
## NULL
## 
## 
## [[4]]
## [[4]][[1]]

## 
## [[4]][[2]]

## 
## 
## [[5]]
## [[5]][[1]]

## 
## [[5]][[2]]

## 
## 
## [[6]]
## [[6]][[1]]

## 
## [[6]][[2]]

Es denota com a mesura que augmenta la dispersió també hi ha una deviació més pronunciada de la normalitat que comença semblant només petit efecte de cua llarga (long tails) i que es va acusant a mesura que es pren conjunts més dispersos de dades. Probablement és un dels motius de l’anterior indicador de pitjor ajust del model MLE en comparació al clàssic OLS.

També s’observa com al aplicar el model amb la funció moduladora de l’error estàndard sembla que es corregeix bastant l’efecte d’heterogeneitat que es veia en els residus així com la desviació de la normalitat.

PROVES DE VARIACIÓ DE LOTS/TEMPS He

S’utilitzen les fòrmules ja generades per graficar el heatmap de variació de variants sobre el model trobat com a òptim dins les proves realitzades:

set.seed(7982)
i <- seq(i3,i10)
j <- seq(j4,j10)

DS_He1_df <- sapply(1:i[length(i)],function(a){
  m_lot <- runif(1, min=-10, max=10)
  sapply(1:j[length(j)], function(j){
    runif(n=1, min=m_lot-(j*2), max=m_lot+(j*2))
  })
  })
  
  #Es considera un Dataset amb el màxim de dades i després s'extreuen les que interessin a cada punt del rang a testar

fm_ML_He1 <-  formula(He1 ~ 1 + Temps)
re_RI <- formula(~1|Lot)
fm_W <- formula(~Temps)

meREML <- "REML"
HeatMapRange_ML("He1",DS_He1_df,fm_ML_He1,re_RI,meREML,2,fm_W)

set.seed(7982)
DS_He2_df <- sapply(1:i[length(i)],function(a){
  m_lot <- runif(1, min=-10, max=10)
  sapply(1:j[length(j)], function(j){
    runif(n=1, min=m_lot-(j^2), max=m_lot+(j^2))
  })
  })
  
  #Es considera un Dataset amb el màxim de dades i després s'extreuen les que interessin a cada punt del rang a testar

fm_ML_He2 <-  formula(He2 ~ 1 + Temps)
fm_W <- formula(~Temps)

meREML <- "REML"
HeatMapRange_ML("He2",DS_He2_df,fm_ML_He2,re_RI,meREML,2,fm_W)

set.seed(7982)
DS_He3_df <- sapply(1:i[length(i)],function(a){
  m_lot <- runif(1, min=-10, max=10)
  sapply(1:j[length(j)], function(j){
    runif(n=1, min=m_lot-(exp(j)), max=m_lot+(exp(j)))
  })
  })
  
  #Es considera un Dataset amb el màxim de dades i després s'extreuen les que interessin a cada punt del rang a testar

fm_ML_He3 <-  formula(He3 ~ 1 + Temps)
fm_W <- formula(~Temps)

meREML <- "REML"
HeatMapRange_ML("He3",DS_He3_df,fm_ML_He3,re_RI,meREML,2,fm_W)

Per generar els heatmaps de comparació entre la variació de lots i temps d’anàlisis s’ha utilitzat el segon model més eficient que s’ha trobat en els anteriors càlculs degut a que en el cas de la funció varPower, al calcular un nombre tan gran de datasets, en alguns d’ells no troba convergència en les iteracions per generar un resultat de la funció de variància. La funció utilitzada ha estat la de varExp que a efectes de comprovar els efectes de variar aquests 2 paràmetres ens pot servir de símil correctament.

A part de les diferències ja vistes entre els 3 datasets que provoquen l’augment de dispersió de dades en aquest cas s’observa que:

• Pels paràmetres indicatius de la qualitat del model AIC i BIC sembla que al tenir més lots/Temps augmenta aquest paràmetre disminuïnt l’eficàcia del model. Possiblement l’efecte de dispersió augmenta al afegir més lots i temps i per tant fa més dificil l’ajust d’un model eficient amb els recursos que s’han utilitzat fins ara.
• El paràmetre de l’error estàndard residual semblaria que segueix una patró similar a AIC i BIC, però és curiòs com sembla que pren més importància el paràmetre del nombre de temps. És a dir al variar el nombre de lots sembla no haver un rang molt gran de canvi mentre que en el nombre de temps sembla haver-hi dos zones molt diferenciades amb sigmes alts i baixos.

Tot i que en aquest apartat es podria treure la conclusió que és millor tenir un nombre més limitat de lots i temps, en la realitat sempre és millor disposar del màxim de dades possibles ja que habitualment les dades d’un mateix lot no es comportaran amb una dispersió tan exagerada com s’ha simulat en aquests models. Si es donés aquest cas amb aquest diagnòstic, potser s’hauria de començar a sospitar que hi ha algun factor que no s’ha tingut en compte o algun error en les dades.

MODEL MIXT - DATASET AMB TENDÈNCIA

GENERACIÓ DATASET

Per la generació de datasets s’utilitza sempre la mateixa llavor essent aquesta 7982 i el dataset base Av_DS com a referencia sobre el qual es van generant les noves variables.

GENERACIÓ DATASET So1 (TENDÈNCIA AMB RECTES PARALELES CENTRADES)

Es generen punts seguint una tendència concreta amb els punts centrats en una recta concreta amb poca dispersió (a efectes d’escalat al voltant de 1 d’una recta amb pendent -1) mitjançant la funció runif:

set.seed(7982)
Av_DS$So1 <- as.vector(sapply(1:av_i,function(i){
  sapply(1:av_j, function(j){
    runif(n=1, min=(TA[length(TA)]-TA[j])-1, max=(TA[length(TA)]-TA[j])+1) #Observacions constant de rang petit entre -1 i 1 per simular que no hi Ha tendències ni diferències
})
})
)

kable(head(Av_DS[,c(colnames(Av_DS)[1:3],"So1")]), booktabs = T) %>% kable_styling(latex_options =  c("striped","condensed"))

Lot	NTemps	Temps	So1
Lot 1	Temps 1	0	41.83072
Lot 1	Temps 2	3	39.20318
Lot 1	Temps 3	6	36.61178
Lot 1	Temps 4	9	32.58671
Lot 1	Temps 5	12	30.48684
Lot 1	Temps 6	18	23.88950

GENERACIÓ DATASET So2 (TENDÈNCIA AMB RECTES PARALELES DISPERSES)

El mateix cas que l’anterior, però suposant que hi ha dispersió més gran entre els lots (aquesta simulació seria un dels casos més habituals en situació ideal que es dona en els estudis d’estabilitat si el disseny és correcte):

set.seed(7982)
Av_DS$So2 <- as.vector(sapply(1:av_i, function(i){
  m_lot <- runif(1, min=TA[length(TA)]-30, max=TA[length(TA)]+30)
  sapply(1:av_j,function(j){
    runif(n=1, min=(m_lot-TA[j])-1, max=(m_lot-TA[j])+1)
  })
  }) 
  )

kable(head(Av_DS[,c(colnames(Av_DS)[1:3],"So2")]), booktabs = T) %>% kable_styling(latex_options =  c("striped","condensed"))

Lot	NTemps	Temps	So2
Lot 1	Temps 1	0	37.12479
Lot 1	Temps 2	3	34.53340
Lot 1	Temps 3	6	30.50832
Lot 1	Temps 4	9	28.40845
Lot 1	Temps 5	12	24.81112
Lot 1	Temps 6	18	18.54508

GENERACIÓ DATASET So3 (TENDÈNCIA AMB RECTES NO PARALELES DISPERSES)

En aquest cas es considera certa dispersió en la intercepció dels lots i una variabilitat en la pendent segons el lot dins un rang definit (aquesta simulació seria un dels altres casos habituals en els estudis d’estabilitat en les fases de testatge de possible combinació de lots):

set.seed(7982)
Av_DS$So3 <- as.vector(sapply(1:av_i, function(i){
  m_lot <- runif(1, min=TA[length(TA)]-10, max=TA[length(TA)]+10)
  s_lot <- runif(1, min=0.1, max=2)
  sapply(1:av_j,function(j){
    runif(n=1, min=(m_lot-(s_lot*TA[j]))-1, max=(m_lot-(s_lot*TA[j]))+1)
  })
  })
  )

kable(head(Av_DS[,c(colnames(Av_DS)[1:3],"So3")]), booktabs = T) %>% kable_styling(latex_options =  c("striped","condensed"))

Lot	NTemps	Temps	So3
Lot 1	Temps 1	0	40.91899
Lot 1	Temps 2	3	36.16485
Lot 1	Temps 3	6	33.33592
Lot 1	Temps 4	9	29.00952
Lot 1	Temps 5	12	25.01442
Lot 1	Temps 6	18	17.12526

VISUALITZACIÓ DATASETS GENERATS So

Es visualitza gràficament les diferències entre els datasets generats. Per veure amb més claredat les diferències es visualitza primer mitjançant un gràfic tipus scatterplot on es veu l’error estàndard global i la diferenciació per lot, i posteriorment es visualitza mitjançant els diagrames de caixes per apreciar més fàcilment per on es mouen les dades:

SCATTERPLOTS So

es veu els datasets amb variància homogència.

labsSo_1 <- c("Dataset So1 - Rectes paraleles centrades","Dataset So2 - Rectes paraleles disperses","Dataset Soe - Rectes no paraleles disperses")
labsSo_2 <- c("Dataset So1 s/ lot","Dataset So2 s/ lot","Dataset So3 s/ lot")

SoVars <- c("So1","So2","So3")

Av_DS_plotSo_sc1 <- lapply(1:length(SoVars),function(graph){
  ggplot(Av_DS, aes(x=Temps, y=Av_DS[,SoVars[graph]])) + geom_point() + geom_smooth(linetype="dashed") + 
    labs(title=labsSo_1[graph], y = "Observacions")
  })

Av_DS_plotSo_sc2 <- lapply(1:length(SoVars),function(graph){
  ggplot(Av_DS, aes(x=Temps, y=Av_DS[,SoVars[graph]], color=Lot, group=Lot)) + geom_point() + geom_line() + 
    labs(title=labsSo_2[graph], y = "Observacions")
  })

ydens_DS_plotSo <- lapply(1:length(SoVars),function(graph){
  ggplot(Av_DS, aes(x=Av_DS[,SoVars[graph]])) + geom_density(color="darkblue", fill="lightblue") +
    geom_vline(aes(xintercept=mean(Av_DS[,SoVars[graph]])), color="blue", linetype="dashed", size=1)+ 
    coord_flip() + labs(x = "Observacions")
  })

grid.arrange(Av_DS_plotSo_sc1[[1]], ydens_DS_plotSo[[1]], Av_DS_plotSo_sc1[[2]], ydens_DS_plotSo[[2]],
             Av_DS_plotSo_sc1[[3]], ydens_DS_plotSo[[3]],ncol=2, widths=c(4, 1))

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

grid.arrange(Av_DS_plotSo_sc2[[1]], ydens_DS_plotSo[[1]], Av_DS_plotSo_sc2[[2]], ydens_DS_plotSo[[2]],
             Av_DS_plotSo_sc2[[3]], ydens_DS_plotSo[[3]],ncol=2, widths=c(4, 1))

S’ha intentat diferenciar entre els lots centrats amb certa dispersió (mínima en comparació a la magnitud de la tendència), el conjunt on realment hi ha diferències en la intercepció, però amb mateixa pendent i l’últim conjunt on hi ha diferències de pendents.

BOXPLOTS So

es veu boxplots dels datasets.

Av_DS_plotSo_bp1 <- lapply(1:length(SoVars),function(graph){
  ggplot(Av_DS, aes(x=NTemps, y=Av_DS[,SoVars[graph]])) + 
  geom_boxplot()+ geom_dotplot(binaxis='y', stackdir='center', dotsize=1)+
  labs(title=labsSo_1[graph], y = "Observacions")
  })

grid.arrange(Av_DS_plotSo_bp1[[1]], ydens_DS_plotSo[[1]], Av_DS_plotSo_bp1[[2]], ydens_DS_plotSo[[2]],
             Av_DS_plotSo_bp1[[3]], ydens_DS_plotSo[[3]],
             ncol=2, widths=c(4, 1))

## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

Es reflexa amb claredat en quin cas els lots estan centrats i en quin no.

MODELS LINEALS So

Per l’ajust dels models lineals dels datasets en aquest cas es consideren diversos models diferents per comparar entre ells. Es pretén generar un model lineal simple, un model lineal d’efectes mixtos amb la intercepció aleatòria i també el d’efectes aleatoris de pendent i intercepció. És important tenir en compte que els datasets generats s’han de considerar independents entre ells i per tant no són directament comparables al ajustar els models, pel que l’anàlisi que es segueix és veure les diferències entre l’ajust dels diferents models a un mateix dataset i posteriorment es compara els diferents datasets, però en termes de la bondat d’ajust dels models aplicats:

ML I/IS So

Es genera el model lineal simple i el model amb la variable temps per comparar-los i definir el model de base simple amb la funció creada per aquest fi:

(Comp_I_IS_So1 <- Comp_I_IS_LM("So1",Av_DS))[1:3]

## $`OLS I`
## 
## Call:
## lm(formula = paste(var, "~", 1, sep = " "), data = ds)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -14.313  -7.566   1.096   6.684  11.281 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   31.682      1.235   25.66   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## s: 8.002 on 41 degrees of freedom
## 
## 
## $`OLS IS`
## 
## Call:
## lm(formula = paste(var, "~", "Temps", sep = " "), data = ds)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.93932 -0.41288 -0.05226  0.47529  0.92559 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 42.03741    0.14883  282.46   <2e-16 ***
## Temps       -1.00675    0.01151  -87.45   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## s: 0.5843 on 40 degrees of freedom
## Multiple R-squared: 0.9948,
## Adjusted R-squared: 0.9947 
## F-statistic:  7648 on 1 and 40 DF,  p-value: < 2.2e-16 
## 
## 
## $`Anova comp. efecte S`
## Analysis of Variance Table
## 
## Model 1: So1 ~ Temps
## Model 2: So1 ~ 1
##   Res.Df     RSS Df Sum of Sq      F    Pr(>F)    
## 1     40   13.66                                  
## 2     41 2625.11 -1   -2611.4 7648.3 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Comp_I_IS_So2 <- Comp_I_IS_LM("So2",Av_DS))[1:3]

## $`OLS I`
## 
## Call:
## lm(formula = paste(var, "~", 1, sep = " "), data = ds)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -26.6589  -7.8191   0.9073   7.7518  16.9864 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   20.138      1.602   12.57 1.18e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## s: 10.38 on 41 degrees of freedom
## 
## 
## $`OLS IS`
## 
## Call:
## lm(formula = paste(var, "~", "Temps", sep = " "), data = ds)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.861  -3.343   1.893   5.205   6.953 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  30.6764     1.6665  18.408  < 2e-16 ***
## Temps        -1.0245     0.1289  -7.948  9.3e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## s: 6.543 on 40 degrees of freedom
## Multiple R-squared: 0.6123,
## Adjusted R-squared: 0.6026 
## F-statistic: 63.17 on 1 and 40 DF,  p-value: 9.299e-10 
## 
## 
## $`Anova comp. efecte S`
## Analysis of Variance Table
## 
## Model 1: So2 ~ Temps
## Model 2: So2 ~ 1
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1     40 1712.5                                  
## 2     41 4417.0 -1   -2704.5 63.171 9.299e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Comp_I_IS_So3 <- Comp_I_IS_LM("So3",Av_DS))[1:3]

## $`OLS I`
## 
## Call:
## lm(formula = paste(var, "~", 1, sep = " "), data = ds)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -30.551 -10.069   1.123  10.964  21.600 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   28.221      2.066   13.66   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## s: 13.39 on 41 degrees of freedom
## 
## 
## $`OLS IS`
## 
## Call:
## lm(formula = paste(var, "~", "Temps", sep = " "), data = ds)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.3509  -5.1996  -0.9243   4.7701  24.4059 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  41.8705     2.1350  19.612  < 2e-16 ***
## Temps        -1.3271     0.1651  -8.036 7.06e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## s: 8.382 on 40 degrees of freedom
## Multiple R-squared: 0.6175,
## Adjusted R-squared: 0.608 
## F-statistic: 64.58 on 1 and 40 DF,  p-value: 7.063e-10 
## 
## 
## $`Anova comp. efecte S`
## Analysis of Variance Table
## 
## Model 1: So3 ~ Temps
## Model 2: So3 ~ 1
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1     40 2810.6                                  
## 2     41 7348.2 -1   -4537.6 64.579 7.063e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Per començar s’observa que, com és d’esperar, afegir la variable temps al model és una diferència significativa i una suposada millora al model per explicar la resposta ja que el resultat de l’Anova corresponent entre models dona un p-valor de 0 (cas So1), 0 (cas So2) i 0 (cas So3), pel que per regla general sempre és preferible utilitzar el model més simple possible que no afegeix tants factors redundants que poden fer augmentar la variabilitat.

En quant a l’error estàndard residual s’observa a primera vista que és més baix en el cas de So1 ja que és el conjunt amb menys variabilitat i més fàcil de definir en general:

• So1: 0.5843 (Error residual)
• So2: 6.5431 (Error residual)
• So3: 6.5431 (Error residual)

ML LotxTemps So

En el cas concret de tenir gràfiques amb pendent i possibles combinacions de lot, una de les proves habituals també és la comprovació de si en un model simple l’efecte pendent amb lot és significatiu. És a dir si hi ha equivalència de lots en la pendent o no:

(anova(Creuat_IS_So1 <- lm(So1 ~ 1 + Temps*Lot,data=Av_DS)))

(anova(Creuat_IS_So2 <- lm(So2 ~ 1 + Temps*Lot,data=Av_DS)))

(anova(Creuat_IS_So3 <- lm(So3 ~ 1 + Temps*Lot,data=Av_DS)))

Com s’observa en els 2 primers casos de línies paraleles aquest terme no té significació. Pel contrari sí que en té al model So3 on sí que hi ha diferències de pendents. Aquest diagnòstic ens serviria per veure s’ha d’aplicar algun terme que ajusti correctament el model com podria ser els efectes aleatoris de pendents.

ML RI So

Avançant a partir del model simple amb intercepció ML I es fa l’ajust de models amb efecte aleatori generant la funció generada:

(Comp_IS_RI_So1 <- Comp_IS_RI_LM("So1",Av_DS))[1:3]

## $`OLS IS RI`
## 
## Terms:
##                     Temps Residuals
## Sum of Squares  2611.4501   12.7848
## Deg. of Freedom         1        35
## 
## Residual standard error: 0.604384
## Estimated effects are balanced
## 
## $`MLE IS RI`
## Linear mixed-effects model fit by maximum likelihood
##  Data: ds 
##        AIC     BIC    logLik
##   80.00922 86.9599 -36.00461
## 
## Random effects:
##  Formula: ~1 | Lot
##          (Intercept)  Residual
## StdDev: 1.137151e-05 0.5702469
## 
## Fixed effects: list(fm2) 
##                Value  Std.Error DF   t-value p-value
## (Intercept) 42.03741 0.14882654 35 282.45907       0
## Temps       -1.00675 0.01151162 35 -87.45473       0
##  Correlation: 
##       (Intr)
## Temps -0.796
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.64721254 -0.72403556 -0.09165138  0.83348276  1.62314245 
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`REML IS RI`
## Linear mixed-effects model fit by REML
##  Data: ds 
##       AIC      BIC    logLik
##   90.1237 96.87922 -41.06185
## 
## Random effects:
##  Formula: ~1 | Lot
##          (Intercept)  Residual
## StdDev: 1.376499e-05 0.5843292
## 
## Fixed effects: list(fm2) 
##                Value  Std.Error DF   t-value p-value
## (Intercept) 42.03741 0.14882654 35 282.45907       0
## Temps       -1.00675 0.01151162 35 -87.45473       0
##  Correlation: 
##       (Intr)
## Temps -0.796
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.60751484 -0.70658635 -0.08944259  0.81339589  1.58402484 
## 
## Number of Observations: 42
## Number of Groups: 6

fmSo1_S <- formula(paste("So1","~","Temps"))
(eMLE_IS_So1<-exactRLRT(lme(fmSo1_S, random = ~1|Lot, data = Av_DS, method="ML")))

## Using restricted likelihood evaluated at ML estimators.

## Refit with method="REML" for exact results.

## Warning in model.matrix.default(~m$groups[[n.levels - i + 1]] - 1,
## contrasts.arg = c("contr.treatment", : non-list contrasts argument ignored

## 
##  simulated finite sample distribution of RLRT.
##  
##  (p-value based on 10000 simulated values)
## 
## data:  
## RLRT = 0, p-value = 1

(eREML_IS_So1<-exactRLRT(lme(fmSo1_S, random = ~1|Lot, data = Av_DS, method="REML")))

## Warning in model.matrix.default(~m$groups[[n.levels - i + 1]] - 1,
## contrasts.arg = c("contr.treatment", : non-list contrasts argument ignored

## 
##  simulated finite sample distribution of RLRT.
##  
##  (p-value based on 10000 simulated values)
## 
## data:  
## RLRT = 0, p-value = 1

(Comp_IS_RI_So2 <- Comp_IS_RI_LM("So2",Av_DS))[1:3]

## $`OLS IS RI`
## 
## Terms:
##                     Temps Residuals
## Sum of Squares  2704.4996   12.9518
## Deg. of Freedom         1        35
## 
## Residual standard error: 0.6083181
## Estimated effects are balanced
## 
## $`MLE IS RI`
## Linear mixed-effects model fit by maximum likelihood
##  Data: ds 
##        AIC      BIC    logLik
##   124.2667 131.2174 -58.13336
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept)  Residual
## StdDev:    6.357206 0.5998098
## 
## Fixed effects: list(fm2) 
##                 Value Std.Error DF   t-value p-value
## (Intercept) 30.676361 2.6640135 35  11.51509       0
## Temps       -1.024525 0.0121084 35 -84.61265       0
##  Correlation: 
##       (Intr)
## Temps -0.047
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.03640421 -0.82543522  0.01705584  0.82694033  1.80929780 
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`REML IS RI`
## Linear mixed-effects model fit by REML
##  Data: ds 
##        AIC      BIC    logLik
##   127.4564 134.2119 -59.72818
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept)  Residual
## StdDev:    6.964603 0.6083181
## 
## Fixed effects: list(fm2) 
##                 Value Std.Error DF   t-value p-value
## (Intercept) 30.676361 2.8475056 35  10.77306       0
## Temps       -1.024525 0.0119842 35 -85.48949       0
##  Correlation: 
##       (Intr)
## Temps -0.043
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.00691559 -0.81421877  0.01511302  0.81418512  1.78499773 
## 
## Number of Observations: 42
## Number of Groups: 6

fmSo2_S <- formula(paste("So2","~","Temps"))
(eMLE_IS_So2<-exactRLRT(lme(fmSo2_S, random = ~1|Lot, data = Av_DS, method="ML")))

## Using restricted likelihood evaluated at ML estimators.
## Refit with method="REML" for exact results.

## Warning in model.matrix.default(~m$groups[[n.levels - i + 1]] - 1,
## contrasts.arg = c("contr.treatment", : non-list contrasts argument ignored

## 
##  simulated finite sample distribution of RLRT.
##  
##  (p-value based on 10000 simulated values)
## 
## data:  
## RLRT = 155.87, p-value < 2.2e-16

(eREML_IS_So2<-exactRLRT(lme(fmSo2_S, random = ~1|Lot, data = Av_DS, method="REML")))

## Warning in model.matrix.default(~m$groups[[n.levels - i + 1]] - 1,
## contrasts.arg = c("contr.treatment", : non-list contrasts argument ignored

## 
##  simulated finite sample distribution of RLRT.
##  
##  (p-value based on 10000 simulated values)
## 
## data:  
## RLRT = 155.92, p-value < 2.2e-16

(Comp_IS_RI_So3 <- Comp_IS_RI_LM("So3",Av_DS))[1:3]

## $`OLS IS RI`
## 
## Terms:
##                    Temps Residuals
## Sum of Squares  4537.636   591.197
## Deg. of Freedom        1        35
## 
## Residual standard error: 4.109909
## Estimated effects are balanced
## 
## $`MLE IS RI`
## Linear mixed-effects model fit by maximum likelihood
##  Data: ds 
##        AIC      BIC    logLik
##   263.4209 270.3716 -127.7105
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:    7.106096 4.052425
## 
## Fixed effects: list(fm2) 
##                Value Std.Error DF   t-value p-value
## (Intercept) 41.87051 3.1552307 35  13.27019       0
## Temps       -1.32707 0.0818066 35 -16.22203       0
##  Correlation: 
##       (Intr)
## Temps -0.267
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.04303851 -0.35964158 -0.02602916  0.40132348  2.94667999 
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`REML IS RI`
## Linear mixed-effects model fit by REML
##  Data: ds 
##        AIC      BIC    logLik
##   262.5228 269.2783 -127.2614
## 
## Random effects:
##  Formula: ~1 | Lot
##         (Intercept) Residual
## StdDev:    7.810128 4.109909
## 
## Fixed effects: list(fm2) 
##                Value Std.Error DF   t-value p-value
## (Intercept) 41.87051  3.355905 35  12.47667       0
## Temps       -1.32707  0.080968 35 -16.39014       0
##  Correlation: 
##       (Intr)
## Temps -0.248
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.02648445 -0.33795748 -0.02440195  0.38665664  2.88533664 
## 
## Number of Observations: 42
## Number of Groups: 6

fmSo3_S <- formula(paste("So3","~","Temps"))
(eMLE_IS_So3<-exactRLRT(lme(fmSo3_S, random = ~1|Lot, data = Av_DS, method="ML")))

## Using restricted likelihood evaluated at ML estimators.
## Refit with method="REML" for exact results.

## Warning in model.matrix.default(~m$groups[[n.levels - i + 1]] - 1,
## contrasts.arg = c("contr.treatment", : non-list contrasts argument ignored

## 
##  simulated finite sample distribution of RLRT.
##  
##  (p-value based on 10000 simulated values)
## 
## data:  
## RLRT = 40.621, p-value < 2.2e-16

(eREML_IS_So3<-exactRLRT(lme(fmSo3_S, random = ~1|Lot, data = Av_DS, method="REML")))

## Warning in model.matrix.default(~m$groups[[n.levels - i + 1]] - 1,
## contrasts.arg = c("contr.treatment", : non-list contrasts argument ignored

## 
##  simulated finite sample distribution of RLRT.
##  
##  (p-value based on 10000 simulated values)
## 
## data:  
## RLRT = 40.675, p-value < 2.2e-16

Es genera una taula resum dels resultats obtinguts del test RLRT per començar:

I_IS_RI_So_RLRT <- data.frame(Dataset = paste(rep(paste(SoVars,c("I RI","IS RI")),each=3),c("Clàssic OLS","MLE","REML")))

Comp_RI_So_2 <- list("-",eMLE_IS_So1,eREML_IS_So1,
                     "-",eMLE_IS_So2,eREML_IS_So2,
                     "-",eMLE_IS_So3,eREML_IS_So3)
I_IS_RI_So_RLRT$"p-valor RLRT" <- sapply(1:length(Comp_RI_So_2),function(l1){
  if (class(Comp_RI_So_2[[l1]])!="htest"){
    return("-")
    } else {
      round(Comp_RI_So_2[[l1]]$p.value,4)
    }
  })

kable(I_IS_RI_So_RLRT, booktabs = T) %>% kable_styling(latex_options =  c("striped","condensed"))

Dataset	p-valor RLRT
So1 I RI Clàssic OLS
So1 I RI MLE	1
So1 I RI REML	1
So2 IS RI Clàssic OLS
So2 IS RI MLE	0
So2 IS RI REML	0
So3 I RI Clàssic OLS
So3 I RI MLE	0
So3 I RI REML	0

En aquests resultats com a comentari general es pot apreciar que:

• No sembla haver-hi diferències significatives entre mètodes (MLE; REML).
• Entre els tres conjunts de dades sembla que la dispersió de dades és la que activa la signifació de l’efecte aleatori de la intercepció com és de lògica.

ML RIS So

Es modifica l’efecte aleatori afegint la variable temps:

Comp_IS_RIS_LM <- function(var,ds) {
  fm1 <- formula(paste(var,"~",1," + Temps","+ Error(Lot/Temps)"))
  ML_IS_RIS_sOLS <- summary(ML_IS_RIS_OLS <- aov(fm1, data = ds))
  ML_IS_RIS_OLS
  ML_IS_RIS_OLS_err <- sigma(summary.lm(ML_IS_RIS_OLS$Within))
  
  fm2 <- formula(paste(var,"~","Temps"))
  
  ML_IS_RIS_sMLE <- summary(ML_IS_RIS_MLE <- lme(fm2, random = ~1+Temps|Lot, data = ds, control = ctrl, method="ML"))
  
  ML_IS_RIS_sREML <- summary(ML_IS_RIS_REML <- lme(fm2,random = ~1+Temps|Lot, data = ds, control = ctrl, method="REML"))
  
    l <- list(ML_IS_RIS_OLS$Within,ML_IS_RIS_sMLE,ML_IS_RIS_sREML,ML_IS_RIS_OLS,ML_IS_RIS_MLE,ML_IS_RIS_REML)
  names(l) <- c("OLS IS RIS","MLE IS RIS","REML IS RIS","OLS IS RIS base","MLE IS RIS base","REML IS RIS base")
  
  return(l)
} 

#No es pot aplicar la funció exactRLRT quan hi ha més d'un efecte aleatori, pel que es farà la comparació amb el model anterior on només es tenia l'efecte RI mitjançant un ANOVA comparatiu.
  
(Comp_IS_RIS_So1 <- Comp_IS_RIS_LM("So1",Av_DS))[1:3]

## $`OLS IS RIS`
## 
## Terms:
##                 Residuals
## Sum of Squares    10.8835
## Deg. of Freedom        30
## 
## Residual standard error: 0.6023149
## 
## $`MLE IS RIS`
## Linear mixed-effects model fit by maximum likelihood
##  Data: ds 
##        AIC      BIC    logLik
##   84.00964 94.43566 -36.00482
## 
## Random effects:
##  Formula: ~1 + Temps | Lot
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev       Corr  
## (Intercept) 0.0005893073 (Intr)
## Temps       0.0002344358 -0.041
## Residual    0.5702414856       
## 
## Fixed effects: list(fm2) 
##                Value  Std.Error DF   t-value p-value
## (Intercept) 42.03741 0.14882534 35 282.46136       0
## Temps       -1.00675 0.01151193 35 -87.45238       0
##  Correlation: 
##       (Intr)
## Temps -0.796
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.64726584 -0.72395941 -0.09165991  0.83352426  1.62315553 
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`REML IS RIS`
## Linear mixed-effects model fit by REML
##  Data: ds 
##        AIC      BIC    logLik
##   94.12393 104.2572 -41.06197
## 
## Random effects:
##  Formula: ~1 + Temps | Lot
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev       Corr  
## (Intercept) 0.0008666710 (Intr)
## Temps       0.0002152367 -0.031
## Residual    0.5843245686       
## 
## Fixed effects: list(fm2) 
##                Value  Std.Error DF   t-value p-value
## (Intercept) 42.03741 0.14882579 35 282.46050       0
## Temps       -1.00675 0.01151187 35 -87.45287       0
##  Correlation: 
##       (Intr)
## Temps -0.796
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.60755693 -0.70652387 -0.08945138  0.81343297  1.58403234 
## 
## Number of Observations: 42
## Number of Groups: 6

(Comp_IS_RIS_So2 <- Comp_IS_RIS_LM("So2",Av_DS))[1:3]

## $`OLS IS RIS`
## 
## Terms:
##                 Residuals
## Sum of Squares   10.77873
## Deg. of Freedom        30
## 
## Residual standard error: 0.5994089
## 
## $`MLE IS RIS`
## Linear mixed-effects model fit by maximum likelihood
##  Data: ds 
##        AIC      BIC    logLik
##   127.1721 137.5981 -57.58606
## 
## Random effects:
##  Formula: ~1 + Temps | Lot
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev     Corr  
## (Intercept) 6.21011469 (Intr)
## Temps       0.01558689 0.792 
## Residual    0.58599457       
## 
## Fixed effects: list(fm2) 
##                 Value Std.Error DF   t-value p-value
## (Intercept) 30.676361 2.6023751 35  11.78783       0
## Temps       -1.024525 0.0135076 35 -75.84827       0
##  Correlation: 
##       (Intr)
## Temps 0.341 
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.1207227 -0.6926417  0.1027614  0.6612646  1.8739939 
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`REML IS RIS`
## Linear mixed-effects model fit by REML
##  Data: ds 
##        AIC      BIC    logLik
##   130.3496 140.4829 -59.17481
## 
## Random effects:
##  Formula: ~1 + Temps | Lot
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev     Corr  
## (Intercept) 6.80187843 (Intr)
## Temps       0.01828676 0.733 
## Residual    0.59193641       
## 
## Fixed effects: list(fm2) 
##                 Value Std.Error DF   t-value p-value
## (Intercept) 30.676361 2.7809450 35  11.03091       0
## Temps       -1.024525 0.0138465 35 -73.99178       0
##  Correlation: 
##       (Intr)
## Temps 0.358 
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.0668230 -0.7105783  0.1075472  0.6535442  1.8239090 
## 
## Number of Observations: 42
## Number of Groups: 6

(Comp_IS_RIS_So3 <- Comp_IS_RIS_LM("So3",Av_DS))[1:3]

## $`OLS IS RIS`
## 
## Terms:
##                 Residuals
## Sum of Squares   8.074923
## Deg. of Freedom        30
## 
## Residual standard error: 0.5188103
## 
## $`MLE IS RIS`
## Linear mixed-effects model fit by maximum likelihood
##  Data: ds 
##        AIC      BIC    logLik
##   153.0026 163.4287 -70.50132
## 
## Random effects:
##  Formula: ~1 + Temps | Lot
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev    Corr  
## (Intercept) 6.3922707 (Intr)
## Temps       0.4750688 -0.191
## Residual    0.5188103       
## 
## Fixed effects: list(fm2) 
##                Value Std.Error DF   t-value p-value
## (Intercept) 41.87051 2.6775045 35 15.637885       0
## Temps       -1.32707 0.1990113 35 -6.668313       0
##  Correlation: 
##       (Intr)
## Temps -0.193
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -1.559095887 -0.575570678 -0.006151924  0.555801370  1.903070973 
## 
## Number of Observations: 42
## Number of Groups: 6 
## 
## $`REML IS RIS`
## Linear mixed-effects model fit by REML
##  Data: ds 
##        AIC      BIC    logLik
##   150.5445 160.6778 -69.27226
## 
## Random effects:
##  Formula: ~1 + Temps | Lot
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev    Corr  
## (Intercept) 7.0038777 (Intr)
## Temps       0.5205323 -0.191
## Residual    0.5188103       
## 
## Fixed effects: list(fm2) 
##                Value Std.Error DF   t-value p-value
## (Intercept) 41.87051 2.8623728 35 14.627902       0
## Temps       -1.32707 0.2127521 35 -6.237635       0
##  Correlation: 
##       (Intr)
## Temps -0.193
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -1.558916061 -0.581202613 -0.004035162  0.554407903  1.897107378 
## 
## Number of Observations: 42
## Number of Groups: 6

#NO ES PODEN FER ANOVAS DE RESULTATS DE FUNCIONS JA QUE NECESSITA LES FORMULES ORIGINALS, ES FA INDIVIDUALMENT AMB ELS MODELS AJUSTATS PER REML

So1_RI <- lme(So1 ~ 1 + Temps, random = ~1|Lot,control=ctrl, method="REML",data=Av_DS)
So1_RIS <- lme(So1 ~ 1 + Temps, random = ~1+Temps|Lot,control=ctrl, method="REML",data=Av_DS)
(So1_RI_RIS_an <- anova(So1_RI,So1_RIS))

So2_RI <- lme(So2 ~ 1 + Temps, random = ~1|Lot,control=ctrl, method="REML",data=Av_DS)
So2_RIS <- lme(So2 ~ 1 + Temps, random = ~1+Temps|Lot,control=ctrl, method="REML",data=Av_DS)
(So2_RI_RIS_an <- anova(So2_RI,So2_RIS))

So3_RI <- lme(So3 ~ 1 + Temps, random = ~1|Lot,control=ctrl, method="REML",data=Av_DS)
So3_RIS <- lme(So3 ~ 1 + Temps, random = ~1+Temps|Lot,control=ctrl, method="REML",data=Av_DS)
(So3_RI_RIS_an <- anova(So3_RI,So3_RIS))

Aquest possiblemenet seria un dels punts claus per entendre la utilitat de la utilització dels efectes aleatoris amb la covariant inclosa:

• So1: El model centrat no hi ha significació entre els models RI i RIS ja que amb l’efecte aleatori dels lots o inclús sense aquest es podria generar el model fàcilment al voltant de les dades.
• So2: Tot i ser no significatiu ja comença a baixar el p-valor donat que comença a haver més dispersió, però al ser encara rectes paraleles pot ser modelat correctament per l’efecte aleatori de la intercepció del lot.
• So3: Al tenir diferents pendents en els diferents lots pren un efecte molt significatiu el fet d’afegir la variable de l’efecte aleatori relacionat amb la pendent o el que és el mateix amb la relació amb la covariant temps.

A continuació es fa la modificació afegint als models RIS les funcions de modulació de variància per veure si serien fets diferencials per millorar els models:

#NO ES PODEN FER ANOVAS DE RESULTATS DE FUNCIONS JA QUE NECESSITA LES FORMULES ORIGINALS, ES FA INDIVIDUALMENT AMB ELS MODELS AJUSTATS PER REML

So1_RI <- lme(So1 ~ 1 + Temps, random = ~1|Lot,control=ctrl, method="REML",data=Av_DS)
So1_RIS <- lme(So1 ~ 1 + Temps, random = ~1+Temps|Lot,control=ctrl, method="REML",data=Av_DS)
So1_RIS_w1 <- update(So1_RIS,control=ctrl,weights= varExp(form = ~Temps))

## Warning in logLik.reStruct(object, conLin): Singular precision matrix in
## level -1, block 1

## Warning in logLik.reStruct(object, conLin): Singular precision matrix in
## level -1, block 1

(So1_RI_RIS_an <- anova(So1_RIS,So1_RIS_w1))

So2_RI <- lme(So2 ~ 1 + Temps, random = ~1|Lot,control=ctrl, method="REML",data=Av_DS)
So2_RIS <- lme(So2 ~ 1 + Temps, random = ~1+Temps|Lot,control=ctrl, method="REML",data=Av_DS)
So2_RIS_w1 <- update(So2_RIS,control=ctrl,weights= varExp(form = ~Temps))

## Warning in logLik.reStruct(object, conLin): Singular precision matrix in
## level -1, block 1

(So2_RI_RIS_an <- anova(So2_RIS,So2_RIS_w1))

So3_RI <- lme(So3 ~ 1 + Temps, random = ~1|Lot,control=ctrl, method="REML",data=Av_DS)
So3_RIS <- lme(So3 ~ 1 + Temps, random = ~1+Temps|Lot,control=ctrl, method="REML",data=Av_DS)
So3_RIS_w1 <- update(So3_RIS,control=ctrl,weights=varPower())

(So3_RI_RIS_an <- anova(So3_RIS,So3_RIS_w1))

So3_RIS_w2 <- update(So3_RIS,control=ctrl,weights=varExp(form = ~Temps))

(So3_RI_RIS_an2 <- anova(So3_RIS,So3_RIS_w2))

En els casos on ha estat possible s’ha aplicat la funció varPower i els casos on no s’ha arribat a cap convergència en les iteracions s’ha optat per l’aplicació de varExp en funció del temps.

Les conclusions no són molt clares en aquest punt, però semblaria que en l’únic cas on s’ha arribat a un model amb varPower dona un resultat no significatiu, però amb un p-valor bastant baix (NA, 0.2194) al afegir aquest terme al model. La resta provat amb varExp resulten no significatius.

De fet, de lògica, no seria necessari en aquest tipus de conjunt simulats la funció de modular l’error estàndard ja que s’ha dissenyat dins de cada lot una dispersió relativament petita, pel que habitualment no es veuria necessari d’aplicar.

RESUM I DIAGNÒSTIC MODELS So

En el codi resum que es mostra a continuació es comparen els difeents models generats per veure els indicadors AIC i BIC, així com l’error estàndard residual:

#So1
So1_R0 <- list(Comp_I_IS_So1[5],Comp_IS_RI_So1[c(1,5:6)],Comp_IS_RIS_So1[c(1,5:6)])
So1_R <- ConvList(So1_R0)
So1_R[[length(So1_R)+1]] <- So1_RIS_w1
names(So1_R)[length(So1_R)] <- "So1 Model IS REML RIS W"
ML_Resum(So1_R,"So1")

## [[1]]
##                               Sigma     AIC      BIC
## So1--OLS IS base             0.5843 78.0092  83.2222
## So1--OLS IS RI               0.6044 68.8941  72.0612
## So1--MLE IS RI base          0.5702 80.0092  86.9599
## So1--REML IS RI base         0.5843 90.1237  96.8792
## So1--OLS IS RIS              0.6023 56.7178  58.1190
## So1--MLE IS RIS base         0.5702 84.0096  94.4357
## So1--REML IS RIS base        0.5843 94.1239 104.2572
## So1--So1 Model IS REML RIS W 0.6139 96.1129 107.9351
## 
## [[2]]

#So2
So2_R0 <- list(Comp_I_IS_So2[5],Comp_IS_RI_So2[c(1,5:6)],Comp_IS_RIS_So2[c(1,5:6)])
So2_R <- ConvList(So2_R0)
So2_R[[length(So2_R)+1]] <- So2_RIS_w1
names(So2_R)[length(So2_R)] <- "So2 Model IS REML RIS W"
ML_Resum(So2_R,"So2")

## [[1]]
##                               Sigma      AIC      BIC
## So2--OLS IS base             6.5431 280.9286 286.1416
## So2--OLS IS RI               0.6083  69.3613  72.5283
## So2--MLE IS RI base          0.5998 124.2667 131.2174
## So2--REML IS RI base         0.6083 127.4564 134.2119
## So2--OLS IS RIS              0.5994  56.4276  57.8288
## So2--MLE IS RIS base         0.5860 127.1721 137.5981
## So2--REML IS RIS base        0.5919 130.3496 140.4829
## So2--So2 Model IS REML RIS W 0.5980 132.3462 144.1684
## 
## [[2]]

#So3
So3_R0 <- list(Comp_I_IS_So3[5],Comp_IS_RI_So3[c(1,5:6)],Comp_IS_RIS_So3[c(1,5:6)])
So3_R <- ConvList(So3_R0)
So3_R[[length(So3_R)+1]] <- So3_RIS_w1
names(So3_R)[length(So3_R)] <- "So3 Model IS REML RIS W"
ML_Resum(So3_R,"So3")

## [[1]]
##                               Sigma      AIC      BIC
## So3--OLS IS base             8.3824 301.7370 306.9500
## So3--OLS IS RI               4.1099 206.9143 210.0813
## So3--MLE IS RI base          4.0524 263.4209 270.3716
## So3--REML IS RI base         4.1099 262.5228 269.2783
## So3--OLS IS RIS              0.5188  47.7633  49.1645
## So3--MLE IS RIS base         0.5188 153.0026 163.4287
## So3--REML IS RIS base        0.5188 150.5445 160.6778
## So3--So3 Model IS REML RIS W 0.1601 151.0362 162.8584
## 
## [[2]]

Els indicadors donen prova del que s’ha anat veient en els diferents models:

• So1: Els models tenen poca diferència de qualitat en variància residual. Possiblement les diferències en AIC / BIC són pel tema de la distribució normal que es veurà també en els diagnòstics següents.
• So2: S’aprecia clarament tant en els indicadors AIC/BIC com l’error estàndard residual la reducció al afegir el factor aleatori. Dins els models amb factor aleatori en canvi no s’aprecia diferència entre si és l’efecte d’intercepció aleatòria o intercepció/pendents aleatoris. Finalment les diferències amb AIC/BIC entre el tipus de càlcul de model s’observa possiblement pel compliment o no de normalitat.
• So3: Disminució significativa al afegir els models amb efecte RI i encara més significativa al afegir els RIS. Quedaria bastant clar quin és el camí per enfocar aquest tipus de models amb efectes aleatoris. Tot i la aparent millora que afegeix el terme varPower s’estudiarà cada cas concret si realment aquesta funció dona una millora a l’ajust del model.

Es visualitzen els gràfics de diagnòstic dels models base i REML:

So1_Diag0 <- list(So1_R[c(1,4,7)])

So1_Diag <- ConvList(So1_Diag0)

names(So1_Diag) <- paste(rep("So1",each=length(So1_R[c(1,4,7)])),names(So1_Diag))

lapply(1:length(So1_Diag),function(QQ){
  if (class(So1_Diag[[QQ]])=="lme") {
    list(qqnorm(So1_Diag[[QQ]],abline = c(0,1),main=names(So1_Diag)[QQ]),plot(So1_Diag[[QQ]],main=names(So1_Diag)[QQ]))
  } else {
    list(plot(So1_Diag[[QQ]],which=1,main=names(So1_Diag)[QQ]),plot(So1_Diag[[QQ]],which=2,main=names(So1_Diag)[QQ]))
  }
  })

## Warning in if (class(So1_Diag[[QQ]]) == "lme") {: the condition has length
## > 1 and only the first element will be used

## [[1]]
## [[1]][[1]]
## NULL
## 
## [[1]][[2]]
## NULL
## 
## 
## [[2]]
## [[2]][[1]]

## 
## [[2]][[2]]

## 
## 
## [[3]]
## [[3]][[1]]

## 
## [[3]][[2]]

So2_Diag0 <- list(So2_R[c(1,4,7)])

So2_Diag <- ConvList(So2_Diag0)

names(So2_Diag) <- paste(rep("So2",each=length(So2_R[c(1,4,7)])),names(So2_Diag))

lapply(1:length(So2_Diag),function(QQ){
  if (class(So2_Diag[[QQ]])=="lme") {
    list(qqnorm(So2_Diag[[QQ]],abline = c(0,1),main=names(So2_Diag)[QQ]),plot(So2_Diag[[QQ]],main=names(So2_Diag)[QQ]))
  } else {
    list(plot(So2_Diag[[QQ]],which=1,main=names(So2_Diag)[QQ]),plot(So2_Diag[[QQ]],which=2,main=names(So2_Diag)[QQ]))
  }
  })

## Warning in if (class(So2_Diag[[QQ]]) == "lme") {: the condition has length
## > 1 and only the first element will be used

## [[1]]
## [[1]][[1]]
## NULL
## 
## [[1]][[2]]
## NULL
## 
## 
## [[2]]
## [[2]][[1]]

## 
## [[2]][[2]]

## 
## 
## [[3]]
## [[3]][[1]]

## 
## [[3]][[2]]

So3_Diag0 <- list(So3_R[c(1,4,7)])

So3_Diag <- ConvList(So3_Diag0)

names(So3_Diag) <- paste(rep("So3",each=length(So3_R[c(1,4,7)])),names(So3_Diag))

lapply(1:length(So3_Diag),function(QQ){
  if (class(So3_Diag[[QQ]])=="lme") {
    list(qqnorm(So3_Diag[[QQ]],abline = c(0,1),main=names(So3_Diag)[QQ]),plot(So3_Diag[[QQ]],main=names(So3_Diag)[QQ]))
  } else {
    list(plot(So3_Diag[[QQ]],which=1,main=names(So3_Diag)[QQ]),plot(So3_Diag[[QQ]],which=2,main=names(So3_Diag)[QQ]))
  }
  })

## Warning in if (class(So3_Diag[[QQ]]) == "lme") {: the condition has length
## > 1 and only the first element will be used

## [[1]]
## [[1]][[1]]
## NULL
## 
## [[1]][[2]]
## NULL
## 
## 
## [[2]]
## [[2]][[1]]

## 
## [[2]][[2]]

## 
## 
## [[3]]
## [[3]][[1]]

## 
## [[3]][[2]]

Sobretot en els models So2 i So3 s’observa que al anar afegint paràmetres de millora al model els residus tenen menys desviació respecte a la distribució normal, així com la repartició en el gràfic de residus per homogeneitat de variància ja que d’alguna manera, el model aconsegueix corregir els comportaments addicionals que el model simple no era capaç.

Prenent els models més efectius de la tècnica de màxima versemblança s’extreuen les matrius internes d’interès per analitzar més el detall:

So_B <- list(Comp_IS_RI_So1[[6]],Comp_IS_RI_So2[[6]],Comp_IS_RIS_So3[[6]])
names(So_B) <- c(names(Comp_IS_RI_So1)[6],names(Comp_IS_RI_So2)[6],names(Comp_IS_RIS_So3)[6])

lapply(1:length(So_B),function(B){
  VarCorr(So_B[[B]])
})

## [[1]]
## Lot = pdLogChol(1) 
##             Variance     StdDev      
## (Intercept) 1.894749e-10 1.376499e-05
## Residual    3.414406e-01 5.843292e-01
## 
## [[2]]
## Lot = pdLogChol(1) 
##             Variance   StdDev   
## (Intercept) 48.5056974 6.9646032
## Residual     0.3700509 0.6083181
## 
## [[3]]
## Lot = pdLogChol(1 + Temps) 
##             Variance   StdDev    Corr  
## (Intercept) 49.0543031 7.0038777 (Intr)
## Temps        0.2709538 0.5205323 -0.191
## Residual     0.2691641 0.5188103

lapply(1:length(So_B),function(B){
  getVarCov(So_B[[B]],type="conditional")
})

## [[1]]
## Lot Lot 1 
## Conditional variance covariance matrix
##         1       2       3       4       5       6       7
## 1 0.34144 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## 2 0.00000 0.34144 0.00000 0.00000 0.00000 0.00000 0.00000
## 3 0.00000 0.00000 0.34144 0.00000 0.00000 0.00000 0.00000
## 4 0.00000 0.00000 0.00000 0.34144 0.00000 0.00000 0.00000
## 5 0.00000 0.00000 0.00000 0.00000 0.34144 0.00000 0.00000
## 6 0.00000 0.00000 0.00000 0.00000 0.00000 0.34144 0.00000
## 7 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.34144
##   Standard Deviations: 0.58433 0.58433 0.58433 0.58433 0.58433 0.58433 0.58433 
## 
## [[2]]
## Lot Lot 1 
## Conditional variance covariance matrix
##         1       2       3       4       5       6       7
## 1 0.37005 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## 2 0.00000 0.37005 0.00000 0.00000 0.00000 0.00000 0.00000
## 3 0.00000 0.00000 0.37005 0.00000 0.00000 0.00000 0.00000
## 4 0.00000 0.00000 0.00000 0.37005 0.00000 0.00000 0.00000
## 5 0.00000 0.00000 0.00000 0.00000 0.37005 0.00000 0.00000
## 6 0.00000 0.00000 0.00000 0.00000 0.00000 0.37005 0.00000
## 7 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.37005
##   Standard Deviations: 0.60832 0.60832 0.60832 0.60832 0.60832 0.60832 0.60832 
## 
## [[3]]
## Lot Lot 1 
## Conditional variance covariance matrix
##         1       2       3       4       5       6       7
## 1 0.26916 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
## 2 0.00000 0.26916 0.00000 0.00000 0.00000 0.00000 0.00000
## 3 0.00000 0.00000 0.26916 0.00000 0.00000 0.00000 0.00000
## 4 0.00000 0.00000 0.00000 0.26916 0.00000 0.00000 0.00000
## 5 0.00000 0.00000 0.00000 0.00000 0.26916 0.00000 0.00000
## 6 0.00000 0.00000 0.00000 0.00000 0.00000 0.26916 0.00000
## 7 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.26916
##   Standard Deviations: 0.51881 0.51881 0.51881 0.51881 0.51881 0.51881 0.51881

lapply(1:length(So_B),function(B){
  a<-getVarCov(So_B[[B]],type="marginal")
  b<-cov2cor(a[[1]])
  return(list(a,b))
})

## [[1]]
## [[1]][[1]]
## Lot Lot 1 
## Marginal variance covariance matrix
##            1          2          3          4          5          6
## 1 3.4144e-01 1.8947e-10 1.8947e-10 1.8947e-10 1.8947e-10 1.8947e-10
## 2 1.8947e-10 3.4144e-01 1.8947e-10 1.8947e-10 1.8947e-10 1.8947e-10
## 3 1.8947e-10 1.8947e-10 3.4144e-01 1.8947e-10 1.8947e-10 1.8947e-10
## 4 1.8947e-10 1.8947e-10 1.8947e-10 3.4144e-01 1.8947e-10 1.8947e-10
## 5 1.8947e-10 1.8947e-10 1.8947e-10 1.8947e-10 3.4144e-01 1.8947e-10
## 6 1.8947e-10 1.8947e-10 1.8947e-10 1.8947e-10 1.8947e-10 3.4144e-01
## 7 1.8947e-10 1.8947e-10 1.8947e-10 1.8947e-10 1.8947e-10 1.8947e-10
##            7
## 1 1.8947e-10
## 2 1.8947e-10
## 3 1.8947e-10
## 4 1.8947e-10
## 5 1.8947e-10
## 6 1.8947e-10
## 7 3.4144e-01
##   Standard Deviations: 0.58433 0.58433 0.58433 0.58433 0.58433 0.58433 0.58433 
## 
## [[1]][[2]]
##              1            2            3            4            5
## 1 1.000000e+00 5.549279e-10 5.549279e-10 5.549279e-10 5.549279e-10
## 2 5.549279e-10 1.000000e+00 5.549279e-10 5.549279e-10 5.549279e-10
## 3 5.549279e-10 5.549279e-10 1.000000e+00 5.549279e-10 5.549279e-10
## 4 5.549279e-10 5.549279e-10 5.549279e-10 1.000000e+00 5.549279e-10
## 5 5.549279e-10 5.549279e-10 5.549279e-10 5.549279e-10 1.000000e+00
## 6 5.549279e-10 5.549279e-10 5.549279e-10 5.549279e-10 5.549279e-10
## 7 5.549279e-10 5.549279e-10 5.549279e-10 5.549279e-10 5.549279e-10
##              6            7
## 1 5.549279e-10 5.549279e-10
## 2 5.549279e-10 5.549279e-10
## 3 5.549279e-10 5.549279e-10
## 4 5.549279e-10 5.549279e-10
## 5 5.549279e-10 5.549279e-10
## 6 1.000000e+00 5.549279e-10
## 7 5.549279e-10 1.000000e+00
## 
## 
## [[2]]
## [[2]][[1]]
## Lot Lot 1 
## Marginal variance covariance matrix
##        1      2      3      4      5      6      7
## 1 48.876 48.506 48.506 48.506 48.506 48.506 48.506
## 2 48.506 48.876 48.506 48.506 48.506 48.506 48.506
## 3 48.506 48.506 48.876 48.506 48.506 48.506 48.506
## 4 48.506 48.506 48.506 48.876 48.506 48.506 48.506
## 5 48.506 48.506 48.506 48.506 48.876 48.506 48.506
## 6 48.506 48.506 48.506 48.506 48.506 48.876 48.506
## 7 48.506 48.506 48.506 48.506 48.506 48.506 48.876
##   Standard Deviations: 6.9911 6.9911 6.9911 6.9911 6.9911 6.9911 6.9911 
## 
## [[2]][[2]]
##           1         2         3         4         5         6         7
## 1 1.0000000 0.9924287 0.9924287 0.9924287 0.9924287 0.9924287 0.9924287
## 2 0.9924287 1.0000000 0.9924287 0.9924287 0.9924287 0.9924287 0.9924287
## 3 0.9924287 0.9924287 1.0000000 0.9924287 0.9924287 0.9924287 0.9924287
## 4 0.9924287 0.9924287 0.9924287 1.0000000 0.9924287 0.9924287 0.9924287
## 5 0.9924287 0.9924287 0.9924287 0.9924287 1.0000000 0.9924287 0.9924287
## 6 0.9924287 0.9924287 0.9924287 0.9924287 0.9924287 1.0000000 0.9924287
## 7 0.9924287 0.9924287 0.9924287 0.9924287 0.9924287 0.9924287 1.0000000
## 
## 
## [[3]]
## [[3]][[1]]
## Lot Lot 1 
## Marginal variance covariance matrix
##        1      2      3      4       5       6       7
## 1 49.323 46.962 44.870 42.777  40.685  36.500  32.316
## 2 46.962 47.577 47.654 48.001  48.347  49.039  49.732
## 3 44.870 47.654 50.708 53.224  56.009  61.579  67.148
## 4 42.777 48.001 53.224 58.717  63.671  74.118  84.565
## 5 40.685 48.347 56.009 63.671  71.602  86.657 101.980
## 6 36.500 49.039 61.579 74.118  86.657 112.000 136.810
## 7 32.316 49.732 67.148 84.565 101.980 136.810 171.920
##   Standard Deviations: 7.0231 6.8976 7.121 7.6627 8.4618 10.583 13.112 
## 
## [[3]][[2]]
##           1         2         3         4         5         6         7
## 1 1.0000000 0.9694362 0.8971922 0.7948876 0.6846115 0.4910798 0.3509365
## 2 0.9694362 1.0000000 0.9702041 0.9081688 0.8283350 0.6717812 0.5498928
## 3 0.8971922 0.9702041 1.0000000 0.9754107 0.9295106 0.8170943 0.7191797
## 4 0.7948876 0.9081688 0.9754107 1.0000000 0.9819686 0.9139533 0.8416867
## 5 0.6846115 0.8283350 0.9295106 0.9819686 1.0000000 0.9676601 0.9191761
## 6 0.4910798 0.6717812 0.8170943 0.9139533 0.9676601 1.0000000 0.9859485
## 7 0.3509365 0.5498928 0.7191797 0.8416867 0.9191761 0.9859485 1.0000000

Comentaris sobre les matriux extretes dels conjunts So:

• En els primers resultats s’observen a les matrius D donant poca informació, simplement es veu l’variància residual obtinguda en cada conjunt. Com a detall interessant veure que el tercer conjunt que podria ser el més complex s’ha aconseguit reduïr l’variància residual lleugerament per sota dels dos primers conjunts que en teoria semblarien més fàcils de tractar.
• La segona sèria de matrius són les matrius R que reflexen la matriu de variàncies-covariàncies residuals dels diferents temps que tenen tots els lots. S’està en el cas on es considera l’variància residual homogènia i per tant es veu com és la mateix per tots els temps i correspon a la diagonal.
• La tercera sèria de matrius mostra:
  • La matriu primera de cada punt de la sèria és de tipus marginal on es suma l’variància residual i l’variància degut a l’efecte aleatori en la diagonal, mentre que l’variància degut a l’efecte aleatori forma part de la covariància en aquest cas. Això és degut a que en aquesta matriu es reflexa l’efecte de l’variància total dins un dels individus (equivalent per tots els altres) que té l’variància de l’efecte aleatori sumat a l’variància residual al centrar-nos en un punt, i es relaciona amb la resta de punts corresponents als altres temps amb el factor calculat aleatori de l’variància dins el mateix lot.
  • A la segona matriu es plasma la correlació de variàncies i es reflexa la correlació entre els diferents temps del model per un individu (equivalent als altres).
  • Alguns comentaris d’aquestes 2 matrius:
    • En el primer conjunt So1 sembla tenir poca significació la covariant corresponent a l’efecte aleatori i de fet la correlació entre els temps és molt baixa propera a zero. Contràriament en els conjunts So2 i So3 augmenta considerablement l’efecte aleatori i també es reflexa en la correlació entre els temps que és relativament gran. Aquest fet en principi seria causat pel fet que al primer conjunt possiblement no és necessari l’efecte aleatori per explicar el model al tenir un model molt centrat i paralel, cosa que no passaria en els models So2 i So3 on (de diferent manera) hi ha una major dispersió que requeriria de l’efecte aleatori per ser explicada.
    • Tot i que en el conjunt So1 i So2 no varia la covariància i correlació entre temps, sí que ho fa en el conjunt So3, pel que en aquest cas sí que la funció ha aconseguit trobar certa correlació per poder reflectir en aquesta matriu. En els conjunts So1 i So2 la variació s’entén que és de intercepció i per tant no hi hauria diferències segons el temps. El conjunt So3 on s’aplica també l’efecte aleatori de la pendent dona per fet que els diferents temps estaran correlacionats per poder construir la pendent que toqui en cada moment (de fet ja es veu l’estructura típica on els temps més propers estan més correlacionats que els temps més llunyans i que conicidiria amb un comportament habitual dels models en els estudis d’estabilitat).

PROVES DE VARIACIÓ DE LOTS/TEMPS So

S’utilitzen les fòrmules ja generades per graficar el heatmap de variació de variants sobre el model trobat com a òptim dins les proves realitzades amb el càlcul de màxima versemblança:

set.seed(7982)
i <- seq(i3,i10)
j <- seq(j4,j10)

DS_So1_df <- sapply(1:i[length(i)],function(i){
  sapply(1:j[length(j)], function(j){
    runif(n=1, min=(TA[length(TA)]-TA[j])-1, max=(TA[length(TA)]-TA[j])+1) #Observacions constant de rang petit entre -1 i 1 per simular que no hi Ha tendències ni diferències
  })
  })
  
  #Es considera un Dataset amb el màxim de dades i després s'extreuen les que interessin a cada punt del rang a testar

fm_ML_So1 <-  formula(So1 ~ 1 + Temps)
re_RI <- formula(~1|Lot)

meREML <- "REML"
HeatMapRange_ML("So1",DS_So1_df,fm_ML_So1,re_RI,meREML,0,0)

El cas més perfecte dels lots centrats en la pendent de So1 reflexa com l’eficàcia del model depèn totalment de la variació del nombre de temps i no depèn del nombre de lots utilitzats. Possiblement es dona això degut a que la mitja dels lots sempre es situa allà mateix i únicament queda més o menys definida segons si conté més o menys repeticions o temps d’anàlisi.

set.seed(7982)
DS_So2_df <- sapply(1:i[length(i)], function(i){
  m_lot <- runif(1, min=TA[length(TA)]-30, max=TA[length(TA)]+30)
  sapply(1:j[length(j)],function(j){
    runif(n=1, min=(m_lot-TA[j])-1, max=(m_lot-TA[j])+1)
  })
  }) 

#Es considera un Dataset amb el màxim de dades i després s'extreuen les que interessin a cada punt del rang a testar

fm_ML_So2 <-  formula(So2 ~ 1 + Temps)

meREML <- "REML"
HeatMapRange_ML("So2",DS_So2_df,fm_ML_So2,re_RI,meREML,0,0)

set.seed(7982)
DS_So3_df <- sapply(1:i[length(i)], function(i){
  m_lot <- runif(1, min=TA[length(TA)]-10, max=TA[length(TA)]+10)
  s_lot <- runif(1, min=0.1, max=2)
  sapply(1:j[length(j)],function(j){
    runif(n=1, min=(m_lot-(s_lot*TA[j]))-1, max=(m_lot-(s_lot*TA[j]))+1)
  })
  })
  
#Es considera un Dataset amb el màxim de dades i després s'extreuen les que interessin a cada punt del rang a testar

fm_ML_So3 <-  formula(So3 ~ 1 + Temps)
re_RIS <- formula(~Temps|Lot)

meREML <- "REML"
HeatMapRange_ML("So3",DS_So3_df,fm_ML_So3,re_RIS,meREML,0,0)

A part de les diferències ja vistes entre els 3 datasets en aquest cas s’observa que:

• Pels paràmetres indicatius de la qualitat del model AIC i BIC sembla que al tenir més lots/Temps augmenta aquest paràmetre disminuïnt l’eficàcia del model. Possiblement com en altres conjunts l’efecte de dispersió augmenta al afegir més lots i temps i per tant fa més dificil l’ajust d’un model eficient amb els recursos que s’han utilitzat fins ara.
• El paràmetre de l’error estàndard residual no sembla seguir un patró concret, probablament degut a que tenim una mida mostral no molt gran i el càlcul d’aquests varia molt segons quines dades es posen o treuen.

Com en casos anteriors tot i que en aquest apartat es podria treure la conclusió que és millor tenir un nombre més limitat de lots i temps, en la realitat sempre és millor disposar del màxim de dades possibles ja que habitualment les dades d’un mateix lot no es comportaran amb una dispersió tan exagerada com s’ha simulat en aquests models. En tot cas és millor intentar millorar el model per adaptar-se a aquest augment de dades.