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
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
.
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 |
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 |
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:
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.
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.
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:
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:
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.
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:
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.
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:
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:
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).
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:
S’aprofita la matriu base generada per generar les noves variables utilitzant sempre la mateixa llavor 7982
.
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ó:
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 |
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:
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.
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ó.
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.
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:
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.
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:
#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
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:
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.
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:
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.
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.
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 |
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 |
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 |
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:
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.
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.
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:
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:
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.
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:
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:
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.
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:
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:
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:
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.