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")
usePackage("haven")
usePackage("summarytools")
usePackage("psych")
usePackage("stringr")
usePackage("randomcoloR")
usePackage("MASS")
usePackage("stats")
usePackage("lmtest")
usePackage("faraway")
usePackage("piecewiseSEM")
ctrl <- lmeControl(opt='optim') # Mesura de control per l'aplicació dels paquets lme per evitar errors per incorrecte optimització d'iteracions
Petita funció que ha servir d’utilitat ja creada en els anteriors codis per reorganitzar objectes del tipus llista:
#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’importa el dataset a analitzar i es visualitza una part de la taula:
workingDir <-getwd()
setwd(workingDir)
LAKE0 <- read_sav(paste(workingDir,"/BdD Clinical Trial/LAKE/lake.sav",sep=""))
head(LAKE0[,1:10])
Es mostra el tipus de variables que formen el model:
dim(LAKE0)
## [1] 116 219
table(sapply(LAKE0,class))
##
## character Date haven_labelled numeric
## 3 8 30 178
which(sapply(LAKE0,class)=="character")
## nusuario npac especificar
## 2 3 8
which(sapply(LAKE0,class)=="Date")
## fecha_nac fecha_ini_lake fecha_vih Fecha_0 Fecha_12
## 5 11 19 24 61
## Fecha_24 Fecha_36 Fecha_48
## 98 138 175
La quantitat d’individus inicial és gran i suficient per analitzar de manera correcta semblaria (116 individus). HI ha també una quantitat enorme de variables (219). És d’interès pel model que es vol testar modelitzar amb poques variables d’inici ja que l’objectiu no és l’anàlisi multivariable en sí.
Respecte a la identificació del tipus de variables hi ha:
Grupo
que marca els dos tractaments de l’estudi i sexo
que podria aportar informació interessant al tipus de model a aplicar.Si s’observa quines són les variables numèriques:
which(sapply(LAKE0,class)=="numeric")
## nvisita a19 a28
## 4 9 16
## estadio_VIH_31 a32 Estado
## 17 18 20
## edad week_0 CargaViral_0
## 21 23 25
## CD4A_0 CD4P_0 CD8A_0
## 26 27 28
## CD8P_0 Hematocrito_0 Hemoglobina_0
## 29 30 31
## Plaquetas_0 Leucocitos_0 LinfosTotales_0
## 32 33 34
## Glucosa_mg_0 Urea_mg_0 Creatinina_mumol_0
## 35 36 37
## Sodio_0 Potasio_0 Cloro_0
## 38 39 40
## Calcio_0 Bilirrubina_mumol_0 GPT_0
## 41 42 43
## GOT_0 GGT_0 ProteinasTotales_0
## 44 45 46
## Albumina_0 Colesterol_mg_0 LDL_mg_0
## 47 48 49
## HDL_mg_0 Trigliceridos_mg_0 Amilasa_0
## 50 51 52
## pH_0 Bicarbonato_0 AcidoLactico_0
## 53 54 55
## AcidoPiruvico_0 week_12 CargaViral_12
## 56 60 62
## CD4A_12 CD4P_12 CD8A_12
## 63 64 65
## CD8P_12 Hematocrito_12 Hemoglobina_12
## 66 67 68
## Plaquetas_12 Leucocitos_12 LinfosTotales_12
## 69 70 71
## Glucosa_mg_12 Urea_mg_12 Creatinina_mumol_12
## 72 73 74
## Sodio_12 Potasio_12 Cloro_12
## 75 76 77
## Calcio_12 Bilirrubina_mumol_12 GPT_12
## 78 79 80
## GOT_12 GGT_12 ProteinasTotales_12
## 81 82 83
## Albumina_12 Colesterol_mg_12 LDL_mg_12
## 84 85 86
## HDL_mg_12 Trigliceridos_mg_12 Amilasa_12
## 87 88 89
## pH_12 Bicarbonato_12 AcidoLactico_12
## 90 91 92
## AcidoPiruvico_12 week_24 CargaViral_24
## 93 97 99
## CD4A_24 CD4P_24 CD8A_24
## 100 101 102
## CD8P_24 Hematocrito_24 Hemoglobina_24
## 103 104 105
## Plaquetas_24 Leucocitos_24 LinfosTotales_24
## 106 107 108
## Glucosa_mg_24 Urea_mg_24 Creatinina_mumol_24
## 109 110 111
## Sodio_24 Potasio_24 Cloro_24
## 112 113 114
## Calcio_24 Bilirrubina_mumol_24 GPT_24
## 115 116 117
## GOT_24 GGT_24 ProteinasTotales_24
## 118 119 120
## Albumina_24 Colesterol_mg_24 LDL_mg_24
## 121 122 123
## HDL_mg_24 Trigliceridos_mg_24 Amilasa_24
## 124 125 126
## pH_24 Bicarbonato_24 AcidoLactico_24
## 127 128 129
## AcidoPiruvico_24 week_36 CargaViral_36
## 130 137 139
## CD4A_36 CD4P_36 CD8A_36
## 140 141 142
## CD8P_36 Hematocrito_36 Hemoglobina_36
## 143 144 145
## Plaquetas_36 Leucocitos_36 LinfosTotales_36
## 146 147 148
## Glucosa_mg_36 Urea_mg_36 Creatinina_mumol_36
## 149 150 151
## Sodio_36 Potasio_36 Cloro_36
## 152 153 154
## Calcio_36 Bilirrubina_mumol_36 GPT_36
## 155 156 157
## GOT_36 GGT_36 ProteinasTotales_36
## 158 159 160
## Albumina_36 Colesterol_mg_36 LDL_mg_36
## 161 162 163
## HDL_mg_36 Trigliceridos_mg_36 Amilasa_36
## 164 165 166
## pH_36 Bicarbonato_36 AcidoLactico_36
## 167 168 169
## AcidoPiruvico_36 week_48 CargaViral_48
## 170 174 176
## CD4A_48 CD4P_48 CD8A_48
## 177 178 179
## CD8P_48 Hematocrito_48 Hemoglobina_48
## 180 181 182
## Plaquetas_48 Leucocitos_48 LinfosTotales_48
## 183 184 185
## Glucosa_mg_48 Urea_mg_48 Creatinina_mumol_48
## 186 187 188
## Sodio_48 Potasio_48 Cloro_48
## 189 190 191
## Calcio_48 Bilirrubina_mumol_48 GPT_48
## 192 193 194
## GOT_48 GGT_48 ProteinasTotales_48
## 195 196 197
## Albumina_48 Colesterol_mg_48 LDL_mg_48
## 198 199 200
## HDL_mg_48 Trigliceridos_mg_48 Amilasa_48
## 201 202 203
## pH_48 Bicarbonato_48 AcidoLactico_48
## 204 205 206
## AcidoPiruvico_48 tpo_vih_meses diff_cd4_48_0
## 207 213 215
## diff_cd4p_48_0 diff_col_48_0 diff_HDL_48_0
## 216 217 218
## diff_LDL_48_0
## 219
Es mostra que l’estructura de la matriu és amb els paràmetres col·locats per cada individu. S’ha de modificar aquest format perquè encaixi amb el format necessari per aplicar amb el codi R on s’afegeix la variable temps i cada variable individual com a columna.
L’objectiu d’aquesta pràctica no és l’anàlisi multivariant i per tant interessa una matriu amb poques variables d’inici, amb variables identificadores i la variable temps. Per fer això es genera un nou dataset a partir de l’original fent les eliminacions comentades, càlculs de temps i en el cas de les variables numèriques s’extreuen semialeatòriament 10 variables:
LAKE0a<- LAKE0[,which(sapply(LAKE0,class)!="haven_labelled" & sapply(LAKE0,class)!="Date")]
LAKE0a$Grupo <- LAKE0$Grupo
LAKE0a$sexo <- LAKE0$sexo
#De les variables numèriques es deixa només amb edad i grup com a idenfiticadors i tractament, i amb les que depenen del temps
#Var no depenents del temps
LAKE1a0 <- LAKE0a[,c(which(sapply(LAKE0a,class)=="character"),which(colnames(LAKE0a) == c("edad","Grupo","sexo")))]
LAKE1a0 <- LAKE1a0[,-3] #S'elimina la var classificació que sembla no tenir sentit aquí
#Es converteix a factors els dos termes identificatius per facilitar la interpretació de les dades i es genera una nova columna identificadora que combina el nusuario i ncap
LAKE1a0[,1] <- as.factor(LAKE1a0[,1][[1]])
LAKE1a0[,2] <- as.factor(LAKE1a0[,2][[1]])
LAKE1a0$Ident <- as.factor(paste(LAKE1a0$nusuario,LAKE1a0$npac))
LAKE1a0$Grupo <- as.factor(LAKE1a0$Grupo)# Es converteix el tractament en una variable del tipus factor
levels(LAKE1a0$Grupo) <- list("Treatm.1"=-1, "Treatm.2"=0)
LAKE1a0$sexo <- as.factor(LAKE1a0$sexo) # Es converteix la variable sexo en una variable del tipus factor
levels(LAKE1a0$sexo) <- list("Hombre"=1, "Mujer"=2)
#Var depenents del temps
LAKE1b0 <- LAKE0a[,which(sapply(LAKE0a,class)=="numeric")]
LAKE1b0 <- LAKE1b0[,-grep("diff_",colnames(LAKE1b0))]#S'elimina les variables diff_ que poden fer confondre en aquest punt de reorganització de dades
nvart_string <- length(vart_string <- c("_0","_12","_24","_36","_48"))
vart_pos <- as.vector(sapply(vart_string,function(s){
grep(paste(s),colnames(LAKE1b0))
}))
LAKE1b0 <- LAKE1b0[,vart_pos]#Es conserva només les variables que acaben en _0, _12, _24, _36 i _48 corresponents a les depenents del temps
nvar <- ncol(LAKE1b0)/length(vart_string)
var_names <- str_remove(colnames(LAKE1b0)[1:nvar],"_0")
var_names <- var_names[-c(which(var_names=="week"),which(var_names=="CargaViral"))]#Elimine la variable week i la CargaViral (aquesta última es recupera a posteriori)
#Es tria de manera semialeatoria algunes variables de les ja seleccionades
set.seed(7982)
nsel <- 3
(var_names_sel <- c("CargaViral",var_names[sample(1:length(var_names),nsel,replace=F)]))
## [1] "CargaViral" "Albumina" "CD4P" "Cloro"
nvar_sel <- length(var_names_sel)
LAKE1b0_sel <- LAKE1b0[,c(paste(var_names_sel,rep(vart_string,each=length(var_names_sel)),sep=""))]
LAKE1ab0_prov <-cbind(LAKE1a0,LAKE1b0_sel)
LAKE1ab0_prov_noNA <- na.omit(LAKE1ab0_prov, cols=((ncol(LAKE1a0)+1):ncol(LAKE1ab0_prov_noNA))) #només a les variables numèriques es busca els NA
LAKE1a0_noNA <- LAKE1ab0_prov_noNA[,1:ncol(LAKE1a0)]
LAKE1b0_noNA <- LAKE1ab0_prov_noNA[,(ncol(LAKE1a0)+1):ncol(LAKE1ab0_prov_noNA)]
#Es reorrganitzen les variables (es crea una funció per poder utilitzar en altres casos)
reorg <- function(df,nvar,nTA){
reorg1 <- lapply(1:nvar,function(v){
var_v <- df[,v]
for(r in 1:(nTA-1)){
var_v <- rbind(var_v,df[,(v+(nvar*r))])
}
return(var_v)
})
reorg2 <- sapply(1:nvar,function(v){
as.vector(reorg1[[v]])
})
return(reorg2)
}
reorgLAKE1b0_noNA <- reorg(LAKE1b0_noNA,nvar_sel,nvart_string)
LAKE1b_noNA <- reorgLAKE1b0_noNA
colnames(LAKE1b_noNA) <- var_names_sel
LAKE1a_noNA <- LAKE1a0_noNA[rep(1:nrow(LAKE1a0_noNA),each=nvart_string),]
LAKE1a_noNA <- droplevels(LAKE1a_noNA) #Es netegen les variables factor de nivells no utilitzats
LAKE1a_noNA$Tiempo <- rep(seq(from=0,to=48,by=12),nlevels(LAKE1a_noNA$Ident))#Es genera la variable temps corresponent a les setmanes de la mesura
head(LAKE1_noNA <- cbind(LAKE1a_noNA,LAKE1b_noNA))
str(LAKE1_noNA)
## 'data.frame': 40 obs. of 11 variables:
## $ nusuario : Factor w/ 2 levels "aocampo","gcarosi": 1 1 1 1 1 1 1 1 1 1 ...
## $ npac : Factor w/ 7 levels "001","002","003",..: 2 2 2 2 2 3 3 3 3 3 ...
## $ edad : num 29 29 29 29 29 33 33 33 33 33 ...
## $ Grupo : Factor w/ 2 levels "Treatm.1","Treatm.2": 1 1 1 1 1 1 1 1 1 1 ...
## $ sexo : Factor w/ 2 levels "Hombre","Mujer": 1 1 1 1 1 1 1 1 1 1 ...
## $ Ident : Factor w/ 8 levels "aocampo 001",..: 2 2 2 2 2 3 3 3 3 3 ...
## $ Tiempo : num 0 12 24 36 48 0 12 24 36 48 ...
## $ CargaViral: num 48400 40 40 40 40 190000 87 40 47 40 ...
## $ Albumina : num 4.39 4.54 4.51 4.67 4.59 4.09 4.06 4.08 3.94 4.27 ...
## $ CD4P : num 12.6 15 26 25 26 19 22 24 29 24 ...
## $ Cloro : num 107 103 104 104 108 105 104 103 106 103 ...
LAKE1 <- LAKE1_noNA #Es deixa un nom més fàcil d'utilitzar
#S'elimina de moment npac i nusuari ja que ja s'ha creat la Ident
LAKE1 <- LAKE1[,3:ncol(LAKE1)]
Al depurar el dataset per acomodar-lo als propòsits d’aquesta part pràctica s’ha quedat en un conjunt de dades amb 8 subjectes i 9 variables.
Es realitzen els anàlisis descriptius bàsics al dataset resultant per donar una visió més clara de les variables que intervenen:
summary(LAKE1)
## edad Grupo sexo Ident Tiempo
## Min. :29.00 Treatm.1:20 Hombre:35 aocampo 001: 5 Min. : 0
## 1st Qu.:32.00 Treatm.2:20 Mujer : 5 aocampo 002: 5 1st Qu.:12
## Median :34.00 aocampo 003: 5 Median :24
## Mean :35.12 aocampo 004: 5 Mean :24
## 3rd Qu.:37.50 gcarosi 003: 5 3rd Qu.:36
## Max. :42.00 gcarosi 006: 5 Max. :48
## (Other) :10
## CargaViral Albumina CD4P Cloro
## Min. : 40.0 Min. :3.100 Min. : 3.90 Min. : 96.0
## 1st Qu.: 40.0 1st Qu.:4.075 1st Qu.:17.00 1st Qu.:101.0
## Median : 50.0 Median :4.400 Median :23.30 Median :104.0
## Mean : 38735.7 Mean :4.302 Mean :22.11 Mean :103.1
## 3rd Qu.: 194.2 3rd Qu.:4.598 3rd Qu.:26.00 3rd Qu.:105.0
## Max. :500000.0 Max. :5.060 Max. :35.00 Max. :108.0
##
by(LAKE1,LAKE1$Grupo,summary)
## LAKE1$Grupo: Treatm.1
## edad Grupo sexo Ident Tiempo
## Min. :29.0 Treatm.1:20 Hombre:20 aocampo 002:5 Min. : 0
## 1st Qu.:32.0 Treatm.2: 0 Mujer : 0 aocampo 003:5 1st Qu.:12
## Median :34.5 aocampo 004:5 Median :24
## Mean :35.0 gcarosi 007:5 Mean :24
## 3rd Qu.:37.5 aocampo 001:0 3rd Qu.:36
## Max. :42.0 gcarosi 003:0 Max. :48
## (Other) :0
## CargaViral Albumina CD4P Cloro
## Min. : 40.0 Min. :3.600 Min. : 3.90 Min. : 96
## 1st Qu.: 40.0 1st Qu.:4.075 1st Qu.:17.00 1st Qu.:103
## Median : 43.5 Median :4.390 Median :22.15 Median :104
## Mean : 37051.2 Mean :4.327 Mean :20.77 Mean :104
## 3rd Qu.: 136.5 3rd Qu.:4.560 3rd Qu.:24.32 3rd Qu.:106
## Max. :500000.0 Max. :4.890 Max. :29.00 Max. :108
##
## --------------------------------------------------------
## LAKE1$Grupo: Treatm.2
## edad Grupo sexo Ident Tiempo
## Min. :32.00 Treatm.1: 0 Hombre:15 aocampo 001:5 Min. : 0
## 1st Qu.:32.00 Treatm.2:20 Mujer : 5 gcarosi 003:5 1st Qu.:12
## Median :33.50 gcarosi 006:5 Median :24
## Mean :35.25 gcarosi 010:5 Mean :24
## 3rd Qu.:36.75 aocampo 002:0 3rd Qu.:36
## Max. :42.00 aocampo 003:0 Max. :48
## (Other) :0
## CargaViral Albumina CD4P Cloro
## Min. : 40.0 Min. :3.100 Min. : 8.40 Min. : 98.0
## 1st Qu.: 50.0 1st Qu.:4.067 1st Qu.:18.70 1st Qu.:100.0
## Median : 50.0 Median :4.470 Median :23.90 Median :102.0
## Mean : 40420.2 Mean :4.277 Mean :23.45 Mean :102.2
## 3rd Qu.: 194.2 3rd Qu.:4.707 3rd Qu.:28.75 3rd Qu.:104.0
## Max. :303409.0 Max. :5.060 Max. :35.00 Max. :107.0
##
by(LAKE1,LAKE1$sexo,summary)
## LAKE1$sexo: Hombre
## edad Grupo sexo Ident Tiempo
## Min. :29.00 Treatm.1:20 Hombre:35 aocampo 002:5 Min. : 0
## 1st Qu.:32.00 Treatm.2:15 Mujer : 0 aocampo 003:5 1st Qu.:12
## Median :35.00 aocampo 004:5 Median :24
## Mean :35.57 gcarosi 003:5 Mean :24
## 3rd Qu.:42.00 gcarosi 006:5 3rd Qu.:36
## Max. :42.00 gcarosi 007:5 Max. :48
## (Other) :5
## CargaViral Albumina CD4P Cloro
## Min. : 40 Min. :3.100 Min. : 3.90 Min. : 96.0
## 1st Qu.: 40 1st Qu.:4.035 1st Qu.:16.50 1st Qu.:100.5
## Median : 50 Median :4.400 Median :22.30 Median :104.0
## Mean : 41060 Mean :4.303 Mean :21.01 Mean :103.1
## 3rd Qu.: 195 3rd Qu.:4.645 3rd Qu.:25.20 3rd Qu.:105.0
## Max. :500000 Max. :5.060 Max. :31.60 Max. :108.0
##
## --------------------------------------------------------
## LAKE1$sexo: Mujer
## edad Grupo sexo Ident Tiempo
## Min. :32 Treatm.1:0 Hombre:0 aocampo 001:5 Min. : 0
## 1st Qu.:32 Treatm.2:5 Mujer :5 aocampo 002:0 1st Qu.:12
## Median :32 aocampo 003:0 Median :24
## Mean :32 aocampo 004:0 Mean :24
## 3rd Qu.:32 gcarosi 003:0 3rd Qu.:36
## Max. :32 gcarosi 006:0 Max. :48
## (Other) :0
## CargaViral Albumina CD4P Cloro
## Min. : 40 Min. :4.16 Min. :23.0 Min. :102.0
## 1st Qu.: 40 1st Qu.:4.19 1st Qu.:26.0 1st Qu.:102.0
## Median : 84 Median :4.20 Median :31.0 Median :104.0
## Mean : 22467 Mean :4.30 Mean :29.8 Mean :103.4
## 3rd Qu.: 170 3rd Qu.:4.46 3rd Qu.:34.0 3rd Qu.:104.0
## Max. :112000 Max. :4.49 Max. :35.0 Max. :105.0
##
table(subset(LAKE1, LAKE1$Grupo=="Treatm.1")$sexo)
##
## Hombre Mujer
## 20 0
table(subset(LAKE1, LAKE1$Grupo=="Treatm.2")$sexo)
##
## Hombre Mujer
## 15 5
Excepte la variable de càrrega viral que sembla moure’s d’una manera diferent per cada individu, la resta de variables numèriques no mostra en aquest punt una diferència que es percebi amb claredat en el resum de dades. Sí que es veu en aquest cas que la proporció de sexes és 100 %, 0 % (Homes/Dones) pel tractament 1 i 75 %, 25 % (Homes/Dones) pel cas del tractament 2 (en cas que no hi hagi repartiment equitatiu pot ser un problema per analitzar el significat real del factor sexe).
Per analitzar possibles relacions entre les diferents variables es pot fer una primera aproximació de càlcul de correlacions i la visualització de l’anomenat scatterplot matrix:
pairs.panels(LAKE1[,sapply(LAKE1,class)=="numeric"],
method = "pearson", # correlation method
hist.col = "#00AFBB",
density = TRUE, # show density plots
ellipses = TRUE # show correlation ellipses
)
LAKE1cor <- as.data.frame(cor(LAKE1[,sapply(LAKE1,class)=="numeric"]))
LAKE1cor2 <- LAKE1cor[which(LAKE1cor$CargaViral!=1),]
LAKE1maxpred <- rownames(LAKE1cor2[which.max(abs(LAKE1cor2$CargaViral)),])
D’aquest gràfic i càlculs interessa saber quines són les variables que tenen una correlació més alta amb la resposta de la càrrega viral que en aquest cas en resulta una correlació amb les variables edad, Tiempo, Albumina, CD4P, Cloro de 0.22, -0.54, -0.31, -0.58, -0.33 respectivament, prenent com a predictor principal la variable de correlació més alta en valor absolut: CD4P.
Com a proba addicional i veient l’histograma de dades de la variable càrrega viral, es repeteix la matriu de correlació eliminant el temps 0 de tots els subjectes per veure com es modifiquen els resultats:
LAKE1r <- subset(LAKE1,LAKE1$Tiempo!=0)
pairs.panels(LAKE1r[,sapply(LAKE1r,class)=="numeric"],
method = "pearson", # correlation method
hist.col = "#00AFBB",
density = TRUE, # show density plots
ellipses = TRUE # show correlation ellipses
)
La proporció és similar pel que de moment no es modificarà res. Tot i així és possible que s’hagi d’aplicar alguna transformació per poder explicar la variabilitat de la càrrega viral amb models més simples.
p01<-ggplot(LAKE1, aes(x=CargaViral)) +
geom_vline(aes(xintercept=mean(LAKE1$CargaViral)),
color=randomColor(), linetype="dashed", size=1) +
geom_histogram()
p02<-ggplot(LAKE1, aes(x="", y=CargaViral, fill=randomColor())) +
geom_boxplot() + geom_dotplot(binaxis='y', stackdir='center', dotsize=0.5) +
theme(legend.position="none")
p03<-ggplot(LAKE1, aes(sample = CargaViral)) + stat_qq() + stat_qq_line()
grid.arrange(p01,p02,p03)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
La distribució d’aquests valor s’allunya clarament de la normalitat. Es mostra si amb la transformació logarítmica canvia una mica aquests primers anàlisis descriptius:
LAKE1$LogCargaViral <- log(LAKE1$CargaViral)
pairs.panels(LAKE1[,sapply(LAKE1,class)=="numeric"],
method = "pearson", # correlation method
hist.col = "#00AFBB",
density = TRUE, # show density plots
ellipses = TRUE # show correlation ellipses
)
p11<-ggplot(LAKE1, aes(x=LogCargaViral)) +
geom_vline(aes(xintercept=mean(LAKE1$LogCargaViral)),
color=randomColor(), linetype="dashed", size=1) +
geom_histogram()
p12<-ggplot(LAKE1, aes(x="", y=LogCargaViral, fill=randomColor())) +
geom_boxplot() + geom_dotplot(binaxis='y', stackdir='center', dotsize=0.5) +
theme(legend.position="none")
p13<-ggplot(LAKE1, aes(sample = LogCargaViral)) + stat_qq() + stat_qq_line()
grid.arrange(p11,p12,p13)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
Sembla millorar lleugerament. Es prova addicionalment també la transformació de Box-Cox que utilitza 2 funcions una logarítmica i una potencial per determinar la millor transformació en un hipotètic model simple depenent del temps:
bct<-boxcox(CargaViral~Tiempo,lambda = seq(-2, 2, length = 10),data=LAKE1,plotit=T)
lambda<-bct$x[which.max(bct$y)]
lambda
## [1] -0.5050505
y <- LAKE1$CargaViral
gm<-geometric.mean(y)
LAKE1$CargaViral_t<-(y^lambda-1)/(lambda*gm^(lambda-1))
LAKE1$LogCargaViral <- log(LAKE1$CargaViral)
pairs.panels(LAKE1[,sapply(LAKE1,class)=="numeric"],
method = "pearson", # correlation method
hist.col = "#00AFBB",
density = TRUE, # show density plots
ellipses = TRUE # show correlation ellipses
)
p21<-ggplot(LAKE1, aes(x=CargaViral_t)) +
geom_vline(aes(xintercept=mean(LAKE1$CargaViral_t)),
color=randomColor(), linetype="dashed", size=1) +
geom_histogram()
p22<-ggplot(LAKE1, aes(x="", y=CargaViral_t, fill=randomColor())) +
geom_boxplot() + geom_dotplot(binaxis='y', stackdir='center', dotsize=0.5) +
theme(legend.position="none")
p23<-ggplot(LAKE1, aes(sample = CargaViral_t)) + stat_qq() + stat_qq_line()
grid.arrange(p21,p22,p23)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
Millora lleugerament respecte a la transformació logarítmica, tot i que els inconvenients potser no justifiquen mantenir aquesta transformació. En principi es mantindrà únicament la logarítmica i en algun cas es pot fer alguna prova per si aquesta transformació aporta una millora significativa.
Lake1SCPlot.a <- ggplot(LAKE1, aes(x=Tiempo, y=CargaViral, color=Grupo, shape=Grupo)) +
geom_point()+
geom_smooth(method=lm)+
labs(title="Càrrega viral s/ tract.",
x="Tiempo", y = "CargaViral")+
scale_color_brewer(palette="Accent")
Lake1SCPlot.b <- ggplot(LAKE1, aes(x=Tiempo, y=LogCargaViral, color=Grupo, shape=Grupo)) +
geom_point()+
geom_smooth(method=lm)+
labs(title="Càrrega viral s/ tract.",
x="Tiempo", y = "Log CargaViral")+
scale_color_brewer(palette="Accent")
Lake1SCPlot.c <- ggplot(LAKE1, aes(x=Tiempo, y=CargaViral_t, color=Grupo, shape=Grupo)) +
geom_point()+
geom_smooth(method=lm)+
labs(title="Càrrega viral s/ tract.",
x="Tiempo", y = "CargaViral_t")+
scale_color_brewer(palette="Accent")
grid.arrange(Lake1SCPlot.a,Lake1SCPlot.b,Lake1SCPlot.c)
Lake1BPlot.a <- ggplot(LAKE1, aes(x=Grupo, y=CargaViral, fill=Grupo)) +
geom_boxplot()+
geom_dotplot(binaxis='y', stackdir='center', dotsize=1)+
labs(title="Càrrega viral s/ tract.",
x="Tractament", y = "CargaViral")+
scale_fill_brewer(palette="Accent")
Lake1BPlot.b <- ggplot(LAKE1, aes(x=Grupo, y=LogCargaViral, fill=Grupo)) +
geom_boxplot()+
geom_dotplot(binaxis='y', stackdir='center', dotsize=1)+
labs(title="Càrrega viral s/ tract.",
x="Tractament", y = "LogCargaViral")+
scale_fill_brewer(palette="Accent")
Lake1BPlot.c <- ggplot(LAKE1, aes(x=Grupo, y=CargaViral_t, fill=Grupo)) +
geom_boxplot()+
geom_dotplot(binaxis='y', stackdir='center', dotsize=1)+
labs(title="Càrrega viral s/ tract.",
x="Tractament", y = "CargaViral_t")+
scale_fill_brewer(palette="Accent")
grid.arrange(Lake1BPlot.a,Lake1BPlot.b,Lake1BPlot.c)
## `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`.
S’ha deixat un model amb pocs individus i poques variables per poder il·lustrar el context habitual dels estudis d’estabilitat on es solen tenir pocs lots de referència pels primers models de predicció i poques variables predictores a part del temps.
Abans de començar a entrar amb models més complexos es comença ajustant el model més simple per mínims quadrats, però en aquest cas s’aplica per ser més pràctics una aproximació mitjançant algoritmes que segueixen el model stepwise el qual és una combinació de les comparacions ANOVA entre afegir predictors al model base o eliminar-ne al model complet per arribar al model més òptim (des del punt de vista d’aquest sistema). Habitualment no s’arriba al model òptim, però servirà com a punt de partida per elaborar els models més complexos:
# Fit the full model
Lake1.fm0 <- lm(CargaViral ~ . - LogCargaViral - CargaViral_t - Ident, data = LAKE1)
# Stepwise regression model
Lake1.step.model0 <- stepAIC(Lake1.fm0, direction = "both",
trace = F)
(Lake1.step.model0s <- summary(Lake1.step.model0))
##
## Call:
## lm(formula = CargaViral ~ Tiempo + CD4P + Cloro, data = LAKE1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -113727 -55558 -7185 42815 276156
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 910289 478303 1.903 0.0650 .
## Tiempo -1868 958 -1.950 0.0591 .
## CD4P -4985 2505 -1.990 0.0542 .
## Cloro -6948 4748 -1.463 0.1521
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## s: 81390 on 36 degrees of freedom
## Multiple R-squared: 0.4226,
## Adjusted R-squared: 0.3745
## F-statistic: 8.783 on 3 and 36 DF, p-value: 0.0001673
anova(Lake1.step.model0)
# Fit the full model
Lake1.fm1 <- lm(LogCargaViral ~ . - CargaViral - CargaViral_t - Ident, data = LAKE1)
# Stepwise regression model
Lake1.step.model1 <- stepAIC(Lake1.fm1, direction = "both",
trace = F)
(Lake1.step.model1s <- summary(Lake1.step.model1))
##
## Call:
## lm(formula = LogCargaViral ~ Tiempo, data = LAKE1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.4563 -1.6775 -0.0412 1.6012 4.3657
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.7566 0.5966 14.678 < 2e-16 ***
## Tiempo -0.1343 0.0203 -6.616 8.16e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## s: 2.178 on 38 degrees of freedom
## Multiple R-squared: 0.5353,
## Adjusted R-squared: 0.5231
## F-statistic: 43.78 on 1 and 38 DF, p-value: 8.162e-08
# Fit the full model
Lake1.fm2 <- lm(CargaViral_t ~ . - CargaViral - LogCargaViral - Ident, data = LAKE1)
# Stepwise regression model
Lake1.step.model2 <- stepAIC(Lake1.fm2, direction = "both",
trace = F)
(Lake1.step.model2s <- summary(Lake1.step.model2))
##
## Call:
## lm(formula = CargaViral_t ~ Tiempo, data = LAKE1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -647.6 -238.6 11.4 272.2 352.4
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7834.432 81.191 96.494 < 2e-16 ***
## Tiempo -21.784 2.762 -7.887 1.61e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## s: 296.5 on 38 degrees of freedom
## Multiple R-squared: 0.6208,
## Adjusted R-squared: 0.6108
## F-statistic: 62.2 on 1 and 38 DF, p-value: 1.606e-09
Lake1.st1 <- (cbind(round(Lake1.step.model0s$coefficients[-1,-3],4),
rbind(round(Lake1.step.model1s$coefficients[-1,-3],4),rep("--",3),rep("--",3)),
rbind(round(Lake1.step.model2s$coefficients[-1,-3],4),rep("--",3),rep("--",3))))
kable(Lake1.st1) %>%
kable_styling("striped") %>%
add_header_above(c(" " = 1, "M.Simple sense transf. 1" = 3, "M.Simple log transf" = 3, "M.Simple Box-Cox t" = 3))
Estimate | Std. Error | Pr(>|t|) | Estimate | Std. Error | Pr(>|t|) | Estimate | Std. Error | Pr(>|t|) | |
---|---|---|---|---|---|---|---|---|---|
Tiempo | -1867.5819 | 957.9577 | 0.0591 | -0.1343 | 0.0203 | 0 | -21.7842 | 2.7622 | 0 |
CD4P | -4985.2783 | 2505.0233 | 0.0542 | – | – | – | – | – | – |
Cloro | -6947.9474 | 4748.509 | 0.1521 | – | – | – | – | – | – |
Lake1.st2 <- data.frame(Model=c("M.Simple sense transf","M.Simple log transf","M.Simple Box-Cox t"))
Lake1.st2$"Mediana residual" <- c(median(Lake1.step.model0s$residuals),
median(Lake1.step.model1s$residuals),
median(Lake1.step.model2s$residuals))
kable(Lake1.st2) %>%
kable_styling("striped")
Model | Mediana residual |
---|---|
M.Simple sense transf | -7184.6798137 |
M.Simple log transf | -0.0411551 |
M.Simple Box-Cox t | 11.4023266 |
Com es pot observar en els 3 modelatges de model simple:
El que també s’observa és que de moment no hi ha cap prova que diferenciï el factor de tractament amb els paràmetres que s’estan utilitzant pel model, pel que no es pren de moment la consideració que hi hagi diferències entre tractaments (en l’estudi original hauria estat l’objectiu, però en aquest es pren l’objectiu d’aconseguir un model de predicció el màxim de exacte i precís per poder explicar correctament i predir les respostes futures en base al conjunt de dades seleccionat).
Lake1.SM0.sdres <- rstandard(Lake1.step.model0)
Lake1.SM0.studres <- studres(Lake1.step.model0)
par(mfrow=c(1,1))
qqnorm(Lake1.step.model0$residuals,
ylab="Raw Residuals",
xlab="Normal Scores",
main="Model Original - NPP of Raw Residuals")
qqline(Lake1.step.model0$residuals)
qqnorm(Lake1.SM0.sdres,
ylab="Standardized Residuals",
xlab="Normal Scores",
main="Model Original - NPP of Standardized Residuals")
qqline(Lake1.SM0.sdres)
qq01<-qqnorm(Lake1.SM0.studres,
ylab="Studentized Residuals",
xlab="Normal Scores",
main="Model Original - NPP of Studentized Residuals")
qqline(Lake1.SM0.studres)
S0 <- shapiro.test(Lake1.step.model0$residuals)
S0
##
## Shapiro-Wilk normality test
##
## data: Lake1.step.model0$residuals
## W = 0.92419, p-value = 0.01046
Lake1.SM1.sdres <- rstandard(Lake1.step.model1)
Lake1.SM1.studres <- studres(Lake1.step.model1)
par(mfrow=c(1,1))
qqnorm(Lake1.step.model1$residuals,
ylab="Raw Residuals",
xlab="Normal Scores",
main="Model Log - NPP of Raw Residuals")
qqline(Lake1.step.model1$residuals)
qqnorm(Lake1.SM1.sdres,
ylab="Standardized Residuals",
xlab="Normal Scores",
main="Model Log - NPP of Standardized Residuals")
qqline(Lake1.SM1.sdres)
qqnorm(Lake1.SM1.studres,
ylab="Studentized Residuals",
xlab="Normal Scores",
main="Model Log - NPP of Studentized Residuals")
qqline(Lake1.SM1.studres)
S1 <- shapiro.test(Lake1.step.model1$residuals)
S1
##
## Shapiro-Wilk normality test
##
## data: Lake1.step.model1$residuals
## W = 0.94994, p-value = 0.07547
Lake1.SM2.sdres <- rstandard(Lake1.step.model2)
Lake1.SM2.studres <- studres(Lake1.step.model2)
par(mfrow=c(1,1))
qqnorm(Lake1.step.model2$residuals,
ylab="Raw Residuals",
xlab="Normal Scores",
main="Model Box-Cox - NPP of Raw Residuals")
qqline(Lake1.step.model2$residuals)
qqnorm(Lake1.SM2.sdres,
ylab="Standardized Residuals",
xlab="Normal Scores",
main="Model Box-Cox - NPP of Standardized Residuals")
qqline(Lake1.SM2.sdres)
qqnorm(Lake1.SM2.studres,
ylab="Studentized Residuals",
xlab="Normal Scores",
main="Model Box-Cox - NPP of Studentized Residuals")
qqline(Lake1.SM2.studres)
S2 <- shapiro.test(Lake1.step.model2$residuals)
S2
##
## Shapiro-Wilk normality test
##
## data: Lake1.step.model2$residuals
## W = 0.91901, p-value = 0.00717
par(mfrow=c(1,3))
qqnorm(Lake1.SM0.sdres,
ylab="Standardized Residuals",
xlab="Normal Scores",
main="Model Original - Standardized Res.",cex.main=1)
qqline(Lake1.SM0.sdres)
qqnorm(Lake1.SM1.sdres,
ylab="Standardized Residuals",
xlab="Normal Scores",
main="Model Log - Standardized Res.",cex.main=1)
qqline(Lake1.SM1.sdres)
qqnorm(Lake1.SM2.sdres,
ylab="Standardized Residuals",
xlab="Normal Scores",
main="Model Box-Cox - Standardized Res.",cex.main=1)
qqline(Lake1.SM2.sdres)
Sorprenentment els 3 models entren en el que es podria classificar com una distribució normal dels residus del model tot i que no queda suportat per la prova de shapiro amb p-valors del model original, logtransformat i box-cox de 0.0105, 0.0755 i 0.0072 respectivament.
pH01 <- ggplot(LAKE1, aes(x=fitted(Lake1.step.model0), y=residuals(Lake1.step.model0),
shape=Grupo, color=Grupo)) +
geom_jitter(position=position_jitter(1), cex=1.5) +
geom_hline(yintercept=0, linetype="dashed", color = "black") +
labs(title="Fitted vs Residuals (Homocedasticity test)",
x="Fitted", y = "Residuals")
pH02 <-ggplot(LAKE1, aes(x=fitted(Lake1.step.model0), y=sqrt(abs(residuals(Lake1.step.model0))),
shape=Grupo, color=Grupo)) +
geom_jitter(position=position_jitter(1), cex=1.5) +
geom_hline(yintercept=0, linetype="dashed", color = "black") +
labs(title="Fitted vs Residuals (Homocedasticity test)",
x="Fitted", y = expression(sqrt("|residuals|")))
grid.arrange(pH01, pH02,ncol=2)
plot(Lake1.step.model0,which=1)
BP01<-bptest(Lake1.step.model0)
BP01
##
## studentized Breusch-Pagan test
##
## data: Lake1.step.model0
## BP = 17.345, df = 3, p-value = 6e-04
pH11 <- ggplot(LAKE1, aes(x=fitted(Lake1.step.model1), y=residuals(Lake1.step.model1),
shape=Grupo, color=Grupo)) +
geom_jitter(position=position_jitter(1), cex=1.5) +
geom_hline(yintercept=0, linetype="dashed", color = "black") +
labs(title="Fitted vs Residuals (Homocedasticity test)",
x="Fitted", y = "Residuals")
pH12 <-ggplot(LAKE1, aes(x=fitted(Lake1.step.model1), y=sqrt(abs(residuals(Lake1.step.model1))),
shape=Grupo, color=Grupo)) +
geom_jitter(position=position_jitter(1), cex=1.5) +
geom_hline(yintercept=0, linetype="dashed", color = "black") +
labs(title="Fitted vs Residuals (Homocedasticity test)",
x="Fitted", y = expression(sqrt("|residuals|")))
grid.arrange(pH11, pH12,ncol=2)
plot(Lake1.step.model1,which=1)
BP11<-bptest(Lake1.step.model1)
BP11
##
## studentized Breusch-Pagan test
##
## data: Lake1.step.model1
## BP = 19.474, df = 1, p-value = 1.02e-05
pH21 <- ggplot(LAKE1, aes(x=fitted(Lake1.step.model2), y=residuals(Lake1.step.model2),
shape=Grupo, color=Grupo)) +
geom_jitter(position=position_jitter(1), cex=1.5) +
geom_hline(yintercept=0, linetype="dashed", color = "black") +
labs(title="Fitted vs Residuals (Homocedasticity test)",
x="Fitted", y = "Residuals")
pH22 <-ggplot(LAKE1, aes(x=fitted(Lake1.step.model2), y=sqrt(abs(residuals(Lake1.step.model2))),
shape=Grupo, color=Grupo)) +
geom_jitter(position=position_jitter(1), cex=1.5) +
geom_hline(yintercept=0, linetype="dashed", color = "black") +
labs(title="Fitted vs Residuals (Homocedasticity test)",
x="Fitted", y = expression(sqrt("|residuals|")))
grid.arrange(pH21, pH22,ncol=2)
plot(Lake1.step.model2,which=1)
BP21<-bptest(Lake1.step.model2)
BP21
##
## studentized Breusch-Pagan test
##
## data: Lake1.step.model2
## BP = 5.718, df = 1, p-value = 0.01679
par(mfrow=c(1,3))
plot(Lake1.step.model0,which=1)
plot(Lake1.step.model1,which=1)
plot(Lake1.step.model2,which=1)
Els gràfics es veu que:
En aquest punt es podria plantejar algun altre tipus de transformació, plantejar un model no lineal o intentar aplicar funcions dins el model que modulin la variància del model. D’acord amb el que es planteja en el treball interessa en aquest punt posar en pràctica les funcions de modulació per intentar millorar el model, però això es veurà en la part d’aplicació de models mixtes.
S’observa addicionalment que la variància no es diferencia de manera significativa aparentment entre els dos grups de tractament en els 3 models ajustats simples:
by(residuals(Lake1.step.model0),LAKE1$Grupo,var)
## LAKE1$Grupo: Treatm.1
## [1] 7328369519
## --------------------------------------------------------
## LAKE1$Grupo: Treatm.2
## [1] 5213866561
by(residuals(Lake1.step.model1),LAKE1$Grupo,var)
## LAKE1$Grupo: Treatm.1
## [1] 4.704479
## --------------------------------------------------------
## LAKE1$Grupo: Treatm.2
## [1] 4.686232
by(residuals(Lake1.step.model2),LAKE1$Grupo,var)
## LAKE1$Grupo: Treatm.1
## [1] 98189.66
## --------------------------------------------------------
## LAKE1$Grupo: Treatm.2
## [1] 71892.55
En el cas concret del model original hi ha un valor que a nivell de residus estandarditzats o estudentitzats apareix amb una distància fora del comú tant per els gràfics de normalitat com els de homogeneïtat de residus:
LAKE1[which.max(qq01$y),]
LAKE1[which.max(sqrt(abs(residuals(Lake1.step.model0)))),]
Habitualment s’investigaria si és un possible outlier i els motius com es va obtenir el valor per veure si hi ha algun error que s’hagi de corregir. En aquest cas s’opta per fer una eliminació del valor d’aquesta observació i veure les diferències (es grafica directe els QQ plot de residus estandarditzats i estudentitzats per veure més clarament els valors aberrants):
Lake1.step.model0s
##
## Call:
## lm(formula = CargaViral ~ Tiempo + CD4P + Cloro, data = LAKE1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -113727 -55558 -7185 42815 276156
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 910289 478303 1.903 0.0650 .
## Tiempo -1868 958 -1.950 0.0591 .
## CD4P -4985 2505 -1.990 0.0542 .
## Cloro -6948 4748 -1.463 0.1521
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## s: 81390 on 36 degrees of freedom
## Multiple R-squared: 0.4226,
## Adjusted R-squared: 0.3745
## F-statistic: 8.783 on 3 and 36 DF, p-value: 0.0001673
LAKE1c <- LAKE1[-which.max(qq01$y),]
# Fit the full model
Lake1c.fm0 <- lm(CargaViral ~ . - LogCargaViral - CargaViral_t - Ident, data = LAKE1c)
# Stepwise regression model
Lake1c.step.model0 <- stepAIC(Lake1c.fm0, direction = "both",
trace = F)
(Lake1c.step.model0s <- summary(Lake1c.step.model0))
##
## Call:
## lm(formula = CargaViral ~ Tiempo + Albumina, data = LAKE1c)
##
## Residuals:
## Min 1Q Median 3Q Max
## -82200 -40919 -7309 29408 183024
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 210820.1 92558.5 2.278 0.028785 *
## Tiempo -2351.3 568.6 -4.135 0.000203 ***
## Albumina -29172.7 21093.4 -1.383 0.175175
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## s: 59440 on 36 degrees of freedom
## Multiple R-squared: 0.3469,
## Adjusted R-squared: 0.3106
## F-statistic: 9.561 on 2 and 36 DF, p-value: 0.0004673
Lake1c.SM0.sdres <- rstandard(Lake1c.step.model0)
Lake1c.SM0.studres <- studres(Lake1c.step.model0)
par(mfrow=c(1,2))
qqnorm(Lake1c.SM0.sdres,
ylab="Standardized Residuals",
xlab="Normal Scores",
main="Model Original - NPP of Standardized Residuals")
qqline(Lake1c.SM0.sdres)
qq02<-qqnorm(Lake1c.SM0.studres, #Es converteix en objecte per investigar el possible outlier a continuació
ylab="Studentized Residuals",
xlab="Normal Scores",
main="Model Original - NPP of Studentized Residuals")
qqline(Lake1c.SM0.studres)
S0 <- shapiro.test(Lake1c.step.model0$residuals)
S0
##
## Shapiro-Wilk normality test
##
## data: Lake1c.step.model0$residuals
## W = 0.89417, p-value = 0.001508
LAKE1c[which.max(qq02$y),]
Es mostra un gran canvi ja que el sistema deixa de tenir en compte els 2 predictors que tenia en compte abans Tiempo, CD4P, Cloro i té en compte ara Tiempo, Albumina, però es continua veient que al gràfic QQ estudentitzat i estandarditzat punts que estan exageradament per fora de la mostra. Aquí cal plantejar-se si aquests punts són tots outliers, si formen part d’una població diferent o si el model li falta un ajust addicional per contemplar aquesta diferenciació.
Si s’eliminés el segon i tercer punts del model s’obtindrien els següents resultats:
Lake1.step.model0s
##
## Call:
## lm(formula = CargaViral ~ Tiempo + CD4P + Cloro, data = LAKE1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -113727 -55558 -7185 42815 276156
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 910289 478303 1.903 0.0650 .
## Tiempo -1868 958 -1.950 0.0591 .
## CD4P -4985 2505 -1.990 0.0542 .
## Cloro -6948 4748 -1.463 0.1521
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## s: 81390 on 36 degrees of freedom
## Multiple R-squared: 0.4226,
## Adjusted R-squared: 0.3745
## F-statistic: 8.783 on 3 and 36 DF, p-value: 0.0001673
Lake1c.step.model0s
##
## Call:
## lm(formula = CargaViral ~ Tiempo + Albumina, data = LAKE1c)
##
## Residuals:
## Min 1Q Median 3Q Max
## -82200 -40919 -7309 29408 183024
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 210820.1 92558.5 2.278 0.028785 *
## Tiempo -2351.3 568.6 -4.135 0.000203 ***
## Albumina -29172.7 21093.4 -1.383 0.175175
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## s: 59440 on 36 degrees of freedom
## Multiple R-squared: 0.3469,
## Adjusted R-squared: 0.3106
## F-statistic: 9.561 on 2 and 36 DF, p-value: 0.0004673
LAKE1c2 <- LAKE1c[-which.max(qq02$y),]
# Fit the full model
Lake1c2.fm0 <- lm(CargaViral ~ . - LogCargaViral - CargaViral_t - Ident, data = LAKE1c2)
# Stepwise regression model
Lake1c2.step.model0 <- stepAIC(Lake1c2.fm0, direction = "both",
trace = F)
(Lake1c2.step.model0s <- summary(Lake1c2.step.model0))
##
## Call:
## lm(formula = CargaViral ~ Grupo + Tiempo + Cloro, data = LAKE1c2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -73430 -27076 -6571 16815 155899
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -424385.0 337613.6 -1.257 0.217316
## Grupo(Treatm.1) -11930.9 8379.0 -1.424 0.163586
## Tiempo -1808.5 464.3 -3.895 0.000437 ***
## Cloro 4735.1 3259.2 1.453 0.155435
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## s: 47120 on 34 degrees of freedom
## Multiple R-squared: 0.3507,
## Adjusted R-squared: 0.2934
## F-statistic: 6.122 on 3 and 34 DF, p-value: 0.001913
Lake1c2.SM0.sdres <- rstandard(Lake1c2.step.model0)
Lake1c2.SM0.studres <- studres(Lake1c2.step.model0)
par(mfrow=c(1,2))
qqnorm(Lake1c2.SM0.sdres,
ylab="Standardized Residuals",
xlab="Normal Scores",
main="Model Original - NPP of Standardized Residuals")
qqline(Lake1c2.SM0.sdres)
qq03<-qqnorm(Lake1c2.SM0.studres, #Es converteix en objecte per investigar el possible outlier a continuació
ylab="Studentized Residuals",
xlab="Normal Scores",
main="Model Original - NPP of Studentized Residuals")
qqline(Lake1c2.SM0.studres)
S0 <- shapiro.test(Lake1c2.step.model0$residuals)
S0
##
## Shapiro-Wilk normality test
##
## data: Lake1c2.step.model0$residuals
## W = 0.86346, p-value = 0.000276
LAKE1c2[which.max(qq03$y),]
LAKE1c3 <- LAKE1c2[-which.max(qq02$y),]
# Fit the full model
Lake1c3.fm0 <- lm(CargaViral ~ . - LogCargaViral - CargaViral_t - Ident, data = LAKE1c3)
# Stepwise regression model
Lake1c3.step.model0 <- stepAIC(Lake1c3.fm0, direction = "both",
trace = F)
(Lake1c3.step.model0s <- summary(Lake1c3.step.model0))
##
## Call:
## lm(formula = CargaViral ~ Grupo + Tiempo + Cloro, data = LAKE1c3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -75864 -30134 -4811 19688 150741
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -437829.3 335813.3 -1.304 0.201328
## Grupo(Treatm.1) -13618.4 8450.4 -1.612 0.116579
## Tiempo -1882.1 465.7 -4.042 0.000299 ***
## Cloro 4897.7 3242.8 1.510 0.140481
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## s: 46850 on 33 degrees of freedom
## Multiple R-squared: 0.3751,
## Adjusted R-squared: 0.3183
## F-statistic: 6.604 on 3 and 33 DF, p-value: 0.001282
Lake1c3.SM0.sdres <- rstandard(Lake1c3.step.model0)
Lake1c3.SM0.studres <- studres(Lake1c3.step.model0)
par(mfrow=c(1,2))
qqnorm(Lake1c3.SM0.sdres,
ylab="Standardized Residuals",
xlab="Normal Scores",
main="Model Original - NPP of Standardized Residuals")
qqline(Lake1c3.SM0.sdres)
qq04<-qqnorm(Lake1c3.SM0.studres, #Es converteix en objecte per investigar el possible outlier a continuació
ylab="Studentized Residuals",
xlab="Normal Scores",
main="Model Original - NPP of Studentized Residuals")
qqline(Lake1c3.SM0.studres)
S0 <- shapiro.test(Lake1c3.step.model0$residuals)
S0
##
## Shapiro-Wilk normality test
##
## data: Lake1c3.step.model0$residuals
## W = 0.87288, p-value = 0.0005645
LAKE1c3[which.max(qq04$y),]
En el tercer dels modelatges curiosament el sistema afegeix la variable tractament de Grup que fins ara s’havien descartat com a significativa i també és interessant veure que tots els possibles outliers provenen de la mateixa categoria de la variable nusuario
del dataset original. Per comprovar possibles diferències es fan alguns gràfics descriptius de les 3 transformacions:
LAKE1_test <- LAKE1
LAKE1_test$nusuario <- LAKE1_noNA$nusuario
table(LAKE1_test$nusuario)
##
## aocampo gcarosi
## 20 20
Lake1_testSCPlot.a <- ggplot(LAKE1_test, aes(x=Tiempo, y=CargaViral, color=nusuario, shape=nusuario)) +
geom_point()+
geom_smooth(method=lm,se=F)+
labs(title="Càrrega viral s/ nusuario",
x="Tiempo", y = "CargaViral")+
scale_color_brewer(palette="Accent")
Lake1_testSCPlot.b <- ggplot(LAKE1_test, aes(x=Tiempo, y=CargaViral, color=nusuario, shape=Grupo)) +
geom_point()+
geom_smooth(method=lm,se=F)+
labs(title="Càrrega viral s/ tract. y nusuario",
x="Tiempo", y = "CargaViral")+
scale_color_brewer(palette="Accent")
grid.arrange(Lake1_testSCPlot.a,Lake1_testSCPlot.b)
Lake1_testBPlot.a <- ggplot(LAKE1_test, aes(x=nusuario, y=CargaViral, fill=nusuario)) +
geom_boxplot()+
geom_dotplot(binaxis='y', stackdir='center', dotsize=1)+
labs(title="Càrrega viral s/ nusuario",
x="Tractament", y = "CargaViral")+
scale_fill_brewer(palette="Accent")
Lake1_testBPlot.b <- ggplot(LAKE1_test, aes(x=nusuario, y=CargaViral, fill=Grupo)) +
geom_boxplot()+
geom_dotplot(binaxis='y', stackdir='center', dotsize=1)+
labs(title="Càrrega viral s/ tract. y nusuario",
x="Tractament", y = "CargaViral")+
scale_fill_brewer(palette="Accent")
grid.arrange(Lake1_testBPlot.a,Lake1_testBPlot.b)
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
Log:
Lake1_testSCPlot.a <- ggplot(LAKE1_test, aes(x=Tiempo, y=LogCargaViral, color=nusuario, shape=nusuario)) +
geom_point()+
geom_smooth(method=lm,se=F)+
labs(title="Càrrega viral s/ nusuario",
x="Tiempo", y = "LogCargaViral")+
scale_color_brewer(palette="Accent")
Lake1_testSCPlot.b <- ggplot(LAKE1_test, aes(x=Tiempo, y=LogCargaViral, color=nusuario, shape=Grupo)) +
geom_point()+
geom_smooth(method=lm,se=F)+
labs(title="Càrrega viral s/ tract. y nusuario",
x="Tiempo", y = "LogCargaViral")+
scale_color_brewer(palette="Accent")
grid.arrange(Lake1_testSCPlot.a,Lake1_testSCPlot.b)
Lake1_testBPlot.a <- ggplot(LAKE1_test, aes(x=nusuario, y=LogCargaViral, fill=nusuario)) +
geom_boxplot()+
geom_dotplot(binaxis='y', stackdir='center', dotsize=1)+
labs(title="Càrrega viral s/ nusuario",
x="Tractament", y = "LogCargaViral")+
scale_fill_brewer(palette="Accent")
Lake1_testBPlot.b <- ggplot(LAKE1_test, aes(x=nusuario, y=LogCargaViral, fill=Grupo)) +
geom_boxplot()+
geom_dotplot(binaxis='y', stackdir='center', dotsize=1)+
labs(title="Càrrega viral s/ tract. y nusuario",
x="Tractament", y = "LogCargaViral")+
scale_fill_brewer(palette="Accent")
grid.arrange(Lake1_testBPlot.a,Lake1_testBPlot.b)
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
Box-Cox
Lake1_testSCPlot.a <- ggplot(LAKE1_test, aes(x=Tiempo, y=CargaViral_t, color=nusuario, shape=nusuario)) +
geom_point()+
geom_smooth(method=lm,se=F)+
labs(title="Càrrega viral s/ nusuario",
x="Tiempo", y = "CargaViral_t")+
scale_color_brewer(palette="Accent")
Lake1_testSCPlot.b <- ggplot(LAKE1_test, aes(x=Tiempo, y=CargaViral_t, color=nusuario, shape=Grupo)) +
geom_point()+
geom_smooth(method=lm,se=F)+
labs(title="Càrrega viral s/ tract. y nusuario",
x="Tiempo", y = "CargaViral_t")+
scale_color_brewer(palette="Accent")
grid.arrange(Lake1_testSCPlot.a,Lake1_testSCPlot.b)
Lake1_testBPlot.a <- ggplot(LAKE1_test, aes(x=nusuario, y=CargaViral_t, fill=nusuario)) +
geom_boxplot()+
geom_dotplot(binaxis='y', stackdir='center', dotsize=1)+
labs(title="Càrrega viral s/ nusuario",
x="Tractament", y = "CargaViral_t")+
scale_fill_brewer(palette="Accent")
Lake1_testBPlot.b <- ggplot(LAKE1_test, aes(x=nusuario, y=CargaViral_t, fill=Grupo)) +
geom_boxplot()+
geom_dotplot(binaxis='y', stackdir='center', dotsize=1)+
labs(title="Càrrega viral s/ tract. y nusuario",
x="Tractament", y = "CargaViral_t")+
scale_fill_brewer(palette="Accent")
grid.arrange(Lake1_testBPlot.a,Lake1_testBPlot.b)
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
LAKE1c2_test <- LAKE1_test[-c(which.max(qq01$y),which.max(qq02$y)),]
Lake1c2_testBPlot.a <- ggplot(LAKE1c2_test, aes(x=nusuario, y=CargaViral_t, fill=nusuario)) +
geom_boxplot()+
geom_dotplot(binaxis='y', stackdir='center', dotsize=1)+
labs(title="Càrrega viral s/ nusuario (- 2 outl.)",
x="Tractament", y = "CargaViral_t")+
scale_fill_brewer(palette="Accent")
Lake1c2_testBPlot.b <- ggplot(LAKE1c2_test, aes(x=nusuario, y=CargaViral_t, fill=Grupo)) +
geom_boxplot()+
geom_dotplot(binaxis='y', stackdir='center', dotsize=1)+
labs(title="Càrrega viral s/ tract. y nusuario (- 2 outl.)",
x="Tractament", y = "CargaViral_t")+
scale_fill_brewer(palette="Accent")
grid.arrange(Lake1c2_testBPlot.a,Lake1c2_testBPlot.b)
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
No queda clara una diferenciació significativa entre aquest 2 grups de nusuario
, però si que s’observa algunes punts interessants que podrien explicar el que s’ha vist en aquest apartat:
nusuario
que té els valors més extremats i això explicaria possiblement perquè els 3 primers valors aberrants que s’eliminarien son d’aquest grup. Possiblement si es seguís eliminant valors d’aquesta categoria en algun moment es trobarien els primers aberrants de l’altre grup de nusuario
.aocampo
pren més pes i també canvia la proporció entre tractaments pel primer grup. Tot i així es queda en un model no balancejat de temps i nusuario
pel que també es perdria certesa.Sembla que en els dos grups hi ha un conjunt de valors relativament gran que es mostren com aberrants i influents ja que al eliminar-los fan variar significativament el model en cada eliminació. Eliminar-los no seria un procediment correcte sense una justificació com ara una variable oculta que el convertís en una nova població. Com que no es té aquesta evidència, s’opta per no eliminar cap d’ells i amb avançar en paral·lel amb els conjunts transformats que mitigarien aquest efecte.
En els models on apareix més d’un predictor es pot donar el cas que algun predictor sigui una combinació lineal dels altres donant lloc a una matriu \(X^T X\) singular o encara que no hi hagi una combinació lineal perfecte una aproximació a la matriu singular. En aquests casos anomenats de multicolinealitat hi pot haver un problema al tenir més d’un càlcul de mínims quadrats possibles pel coeficient del predictor i això sovint provoca imprecisions en l’estimació de \(\beta\).
A continuació es realitzen algunes proves en el model on apareix més d’un predictor per detectar si hi ha algun problema en aquest sentit:
LAKE1col <- LAKE1[,attr(Lake1.step.model0$terms,"term.labels")]
pairs.panels(LAKE1col,
method = "pearson", # correlation method
hist.col = "#00AFBB",
density = TRUE, # show density plots
ellipses = TRUE # show correlation ellipses
)
LAKE1colcor <- as.data.frame(cor(LAKE1col))
LAKE1colcor2 <- LAKE1colcor[which(LAKE1colcor$Tiempo!=1),] #Es mostra el més gran relacionat amb Temps
LAKE1colmaxpred <- rownames(LAKE1colcor2[which.max(abs(LAKE1colcor2$Tiempo)),])
Semblaria en aquest primer diagnòstic que la variable amb correlació més gran amb la principal covariant Tiempo
seria CD4P.
(LAKE1col.lms<-summary(LAKE1col.lm <- lm(Tiempo ~ .,data=LAKE1col)))
##
## Call:
## lm(formula = Tiempo ~ ., data = LAKE1col)
##
## Residuals:
## Min 1Q Median 3Q Max
## -25.132 -10.896 -1.341 9.187 26.683
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.9199 82.0115 0.255 0.800
## CD4P 1.5745 0.3432 4.587 5.01e-05 ***
## Cloro -0.3077 0.8133 -0.378 0.707
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## s: 13.97 on 37 degrees of freedom
## Multiple R-squared: 0.3735,
## Adjusted R-squared: 0.3396
## F-statistic: 11.03 on 2 and 37 DF, p-value: 0.0001753
També ho s’observa en aquest resultat al comprovar els valors de significació dels diferents predictors secundaris CD4P, Cloro amb p-valors de 1e-04, 0.7073 respectivament.
x <- model.matrix(Lake1.step.model0)[,-1]
e <- eigen(t(x) %*% x)
e$val
## [1] 469390.468 11237.683 1009.149
(k<-sqrt(e$val[1]/e$val)[-1])
## [1] 6.462919 21.566989
Ratios bastant petits i essent per sota de 30 els ratios amb CD4P, Cloro, sobretot el ratio de CD4P amb un ratio de 6.46, indicant possible multicolinealitat.
vif(x)
## Tiempo CD4P Cloro
## 1.596046 1.714961 1.097428
En aquest cas no es mostra una gran inflació de variància en cap dels predictors, pel que no s’aconsegueix una prova molt definitiva de presència de colinealitat.
Es mostra les diferències en l’ajust del model simple:
Lake1.step.model0s
##
## Call:
## lm(formula = CargaViral ~ Tiempo + CD4P + Cloro, data = LAKE1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -113727 -55558 -7185 42815 276156
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 910289 478303 1.903 0.0650 .
## Tiempo -1868 958 -1.950 0.0591 .
## CD4P -4985 2505 -1.990 0.0542 .
## Cloro -6948 4748 -1.463 0.1521
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## s: 81390 on 36 degrees of freedom
## Multiple R-squared: 0.4226,
## Adjusted R-squared: 0.3745
## F-statistic: 8.783 on 3 and 36 DF, p-value: 0.0001673
(Lake1.step.model0col.s<-summary(Lake1.step.model0col<-update(Lake1.step.model0,paste(". ~ . -",colnames(x)[which.min(k)+1],sep=" "))))
##
## Call:
## lm(formula = CargaViral ~ Tiempo + Cloro, data = LAKE1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -113986 -59944 -17035 36231 321205
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1090786.9 488052.3 2.235 0.031540 *
## Tiempo -3015.5 794.8 -3.794 0.000532 ***
## Cloro -9499.9 4751.5 -1.999 0.052954 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## s: 84580 on 37 degrees of freedom
## Multiple R-squared: 0.3591,
## Adjusted R-squared: 0.3244
## F-statistic: 10.36 on 2 and 37 DF, p-value: 0.0002667
(Lake1.step.model0col.an<-anova(Lake1.step.model0col,Lake1.step.model0))
Per un marge molt estret no es consideren significativament diferents els dos models contrastats amb anova amb un p-valor de 0.0542, però sí que es pot veure que en el segon model la variable temps pren més significació passant de un p-valor de 0.0591 a 5e-04, encara que també es redueix el coeficient de determinació \(R^2\) absolut i ajustat de 0.4226/0.3745 a 0.3591/0.3244, perdent % d’explicabilitat de la variància de la resposta.
S’està doncs en una situació de compromís on no hi ha una prova definitiva entre les fetes que doni una direcció molt clara. Com que l’objectiu és crear un model de predicció es sacrificarà tenir incertesa en els coeficients dels predictors en front de tenir un % d’explicabilitat més gran de la resposta agafant el model amb més variables, encara que hi hagi certa colinealitat present entre Tiempo, CD4P, Cloro. En tot cas en les properes passes s’ha de tenir en compte com una opció la d’eliminar aquest predictor amb colinealitat en cas que es compensi amb altres estructures més eficients.
És d’interès comprovar l’efecte que té el factor individu amb la covariant temps per plantejar-nos si en propers passos es considera l’efecte aleatori d’aquest. Per fer això un dels sistemes és afegir el factor individu dins el model i construir el terme creuat amb la covariant temps per veure si aquest factor és significatiu, és a dir si la variació de les observacions amb la covariant temps depenen del factor individu o no (seria l’equivalent a la prova de paral·lelisme dels models de regressió):
#Model original
Lake1.step.model0Is <- summary(Lake1.step.model0I<-update(Lake1.step.model0,. ~ . + Ident*Tiempo)) #S'afegeix el factor creuat al model original
(Lake1.step.model0I.an<-Anova(Lake1.step.model0I))
summary(aov(formula(paste(Lake1.step.model0I$call[2])),data=LAKE1))
## Df Sum Sq Mean Sq F value Pr(>F)
## Tiempo 1 1.197e+11 1.197e+11 23.227 8.16e-05 ***
## CD4P 1 4.065e+10 4.065e+10 7.889 0.0102 *
## Cloro 1 1.418e+10 1.418e+10 2.752 0.1113
## Ident 7 7.980e+10 1.140e+10 2.212 0.0732 .
## Tiempo:Ident 7 4.529e+10 6.470e+09 1.256 0.3165
## Residuals 22 1.134e+11 5.153e+09
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Model amb correcció per colinealitat
Lake1.step.model0colIs <- summary(Lake1.step.model0colI<-update(Lake1.step.model0col,. ~ . + Ident*Tiempo)) #S'afegeix el factor creuat al model original
(Lake1.step.model0colI.an<-Anova(Lake1.step.model0colI))
summary(aov(formula(paste(Lake1.step.model0colI$call[2])),data=LAKE1))
## Df Sum Sq Mean Sq F value Pr(>F)
## Tiempo 1 1.197e+11 1.197e+11 14.573 0.000884 ***
## Cloro 1 2.860e+10 2.860e+10 3.482 0.074851 .
## Ident 7 1.620e+10 2.314e+09 0.282 0.954660
## Tiempo:Ident 7 5.959e+10 8.512e+09 1.036 0.433518
## Residuals 23 1.889e+11 8.213e+09
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Model Log
Lake1.step.model1Is <- summary(Lake1.step.model1I<-update(Lake1.step.model1,. ~ . + Ident*Tiempo)) #S'afegeix el factor creuat al model original
(Lake1.step.model1I.an<-Anova(Lake1.step.model1I))
summary(aov(formula(paste(Lake1.step.model1I$call[2])),data=LAKE1))
## Df Sum Sq Mean Sq F value Pr(>F)
## Tiempo 1 207.75 207.75 31.240 9.43e-06 ***
## Ident 7 10.09 1.44 0.217 0.978
## Tiempo:Ident 7 10.64 1.52 0.229 0.974
## Residuals 24 159.60 6.65
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Model Box-Cox
Lake1.step.model2Is <- summary(Lake1.step.model2I<-update(Lake1.step.model2,. ~ . + Ident*Tiempo)) #S'afegeix el factor creuat al model original
(Lake1.step.model2I.an<-Anova(Lake1.step.model2I))
summary(aov(formula(paste(Lake1.step.model2I$call[2])),data=LAKE1))
## Df Sum Sq Mean Sq F value Pr(>F)
## Tiempo 1 5466812 5466812 46.169 5e-07 ***
## Ident 7 438041 62577 0.528 0.804
## Tiempo:Ident 7 60066 8581 0.072 0.999
## Residuals 24 2841826 118409
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Es considera els resultats de p-valor del factor interacció que inferiors al 5% de significància com a significativa la variació de pendent per individu. En aquest cas els resultats dels p-valors dels diferents models Original, Original corregit per multicolinealitat, Log transformat i amb transformació Box-Cox són de 0.3165, 0.4335, 0.9744 i 0.9992 respectivament.
Segons els resultats en principi no s’aplicarien mesures per ajustar el model a les variacions de pendents dependent del factor individu.
En el model original s’ha considerat més variables predictores a part de la variable temps. Seria interessant veure si la interacció entre aquestes variables té alguna significació que millori l’explicabilitat del model i per això es fan els ajustos corresponents per veure les proves ANOVA de significància. Es fa també el procés de stepwise amb el model complet d’interaccions:
#Original
Inter.0 <- attr(Lake1.step.model0s$terms,"term.labels")
Lake1.step.model0int.s <- summary(Lake1.step.model0int <- update(Lake1.step.model0,paste(". ~ . +",paste(Inter.0,collapse="*"),sep=" ")))
(Lake1.step.model0int2.s2 <- summary(Lake1.step.model0int2 <- stepAIC(Lake1.step.model0int, direction = "both", trace = F)))
##
## Call:
## lm(formula = CargaViral ~ Tiempo + CD4P + Cloro + Tiempo:CD4P +
## Tiempo:Cloro + CD4P:Cloro + Tiempo:CD4P:Cloro, data = LAKE1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -73263 -18903 -1618 24695 147784
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6447992.94 1022547.10 6.306 4.49e-07 ***
## Tiempo -176830.68 58543.48 -3.021 0.00493 **
## CD4P -343452.68 74069.49 -4.637 5.69e-05 ***
## Cloro -60755.43 10234.76 -5.936 1.31e-06 ***
## Tiempo:CD4P 9148.16 2830.75 3.232 0.00285 **
## Tiempo:Cloro 1649.19 574.98 2.868 0.00725 **
## CD4P:Cloro 3282.16 732.55 4.480 8.93e-05 ***
## Tiempo:CD4P:Cloro -86.72 27.74 -3.127 0.00375 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## s: 47430 on 32 degrees of freedom
## Multiple R-squared: 0.8257,
## Adjusted R-squared: 0.7875
## F-statistic: 21.65 on 7 and 32 DF, p-value: 1.856e-10
(Lake1.step.model0int2.an2<-Anova(Lake1.step.model0int2))
(Lake1.step.model0int2.an1_2<-anova(Lake1.step.model0int2,Lake1.step.model0))
if (Lake1.step.model0int2.an1_2$`Pr(>F)`[2]>0.05){
intM0 <- Lake1.step.model0
intM0s <- Lake1.step.model0s
} else {
intM0 <- Lake1.step.model0int2
intM0s <- Lake1.step.model0int2.s2
}
intM0col <- Lake1.step.model0col
intM0cols <- Lake1.step.model0col.s
En el cas dels models transformats no seria necessari en principi ja que només ha quedat un predictor com a model òptim i per tant no té altres predictors per creuar-se. Tot i així es fa la prova on apareixen totes les variables amb interaccions per comprobar la diferència entre models:
#Log
#Step1
Lake1.step.model1int.s <- summary(Lake1.step.model1int <- update(Lake1.step.model1,paste(". ~ . +",paste(Inter.0,collapse="*"),sep=" ")))
(Lake1.step.model1int2.s2 <- summary(Lake1.step.model1int2 <- stepAIC(Lake1.step.model1int, direction = "both", trace = F)))
##
## Call:
## lm(formula = LogCargaViral ~ Tiempo + CD4P + Tiempo:CD4P, data = LAKE1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.7900 -1.0384 -0.2722 0.7064 4.2465
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.642765 1.468711 9.289 4.29e-11 ***
## Tiempo -0.350318 0.069405 -5.047 1.30e-05 ***
## CD4P -0.272305 0.078710 -3.460 0.00141 **
## Tiempo:CD4P 0.010529 0.002956 3.562 0.00106 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## s: 1.89 on 36 degrees of freedom
## Multiple R-squared: 0.6687,
## Adjusted R-squared: 0.6411
## F-statistic: 24.22 on 3 and 36 DF, p-value: 9.352e-09
(Lake1.step.model1int2.an2<-Anova(Lake1.step.model1int2))
(Lake1.step.model1int2.an1_2 <- anova(Lake1.step.model1int2,Lake1.step.model1))
if (Lake1.step.model1int2.an1_2$`Pr(>F)`[2]>0.05){
intM1 <- Lake1.step.model1
intM1s <- Lake1.step.model1s
} else {
intM1 <- Lake1.step.model1int2
intM1s <- Lake1.step.model1int2.s2
}
#Box-Cox
#Step1
Lake1.step.model2int.s <- summary(Lake1.step.model2int <- update(Lake1.step.model2,paste(". ~ . +",paste(Inter.0,collapse="*"),sep=" ")))
(Lake1.step.model2int2.s2 <- summary(Lake1.step.model2int2 <- stepAIC(Lake1.step.model2int, direction = "both", trace = F)))
##
## Call:
## lm(formula = CargaViral_t ~ Tiempo + CD4P + Tiempo:CD4P, data = LAKE1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -631.26 -157.19 31.58 187.95 448.34
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8267.5540 215.9970 38.276 < 2e-16 ***
## Tiempo -47.3143 10.2070 -4.635 4.55e-05 ***
## CD4P -23.5364 11.5756 -2.033 0.0494 *
## Tiempo:CD4P 1.1663 0.4347 2.683 0.0110 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## s: 277.9 on 36 degrees of freedom
## Multiple R-squared: 0.6843,
## Adjusted R-squared: 0.6579
## F-statistic: 26.01 on 3 and 36 DF, p-value: 3.979e-09
(Lake1.step.model2int2.an2<-Anova(Lake1.step.model2int2))
(Lake1.step.model2int2.an1_2 <- anova(Lake1.step.model2int2,Lake1.step.model2))
if (Lake1.step.model2int2.an1_2$`Pr(>F)`[2]>0.05){
intM2 <- Lake1.step.model2
intM2s <- Lake1.step.model2s
} else {
intM2 <- Lake1.step.model2int2
intM2s <- Lake1.step.model2int2.s2
}
intM1s_mod <- as.data.frame(round(intM1s$coefficients[-1,-3],4))
intM1s_mod[4,]<-intM1s_mod[3,]
intM1s_mod[c(3,5:7),]<-NA
intM2s_mod <- as.data.frame(round(intM2s$coefficients[-1,-3],4))
intM2s_mod[4,]<-intM2s_mod[3,]
intM2s_mod[c(3,5:7),]<-NA
Lake1.it1 <- (cbind(round(intM0s$coefficients[-1,-3],4),
intM1s_mod,
intM2s_mod))
kable(Lake1.it1) %>%
kable_styling(bootstrap_options = "striped", full_width = F) %>%
add_header_above(c(" " = 1, "M.Int sense transf. 1" = 3, "M.Int log transf" = 3, "M.Int Box-Cox t" = 3))
Estimate | Std. Error | Pr(>|t|) | Estimate | Std. Error | Pr(>|t|) | Estimate | Std. Error | Pr(>|t|) | |
---|---|---|---|---|---|---|---|---|---|
Tiempo | -176830.6843 | 58543.4804 | 0.0049 | -0.3503 | 0.0694 | 0.0000 | -47.3143 | 10.2070 | 0.0000 |
CD4P | -343452.6770 | 74069.4864 | 0.0001 | -0.2723 | 0.0787 | 0.0014 | -23.5364 | 11.5756 | 0.0494 |
Cloro | -60755.4316 | 10234.7569 | 0.0000 | NA | NA | NA | NA | NA | NA |
Tiempo:CD4P | 9148.1570 | 2830.7465 | 0.0028 | 0.0105 | 0.0030 | 0.0011 | 1.1663 | 0.4347 | 0.0110 |
Tiempo:Cloro | 1649.1938 | 574.9848 | 0.0072 | NA | NA | NA | NA | NA | NA |
CD4P:Cloro | 3282.1583 | 732.5469 | 0.0001 | NA | NA | NA | NA | NA | NA |
Tiempo:CD4P:Cloro | -86.7203 | 27.7370 | 0.0038 | NA | NA | NA | NA | NA | NA |
Lake1.it2 <- data.frame(Model=c("M.Int sense transf","M.Int log transf","M.Int Box-Cox t"))
Lake1.it2$"Mediana residual" <- c(median(intM0s$residuals),
median(intM1s$residuals),
median(intM2s$residuals))
kable(Lake1.it2) %>%
kable_styling("striped")
Model | Mediana residual |
---|---|
M.Int sense transf | -1617.8894864 |
M.Int log transf | -0.2722391 |
M.Int Box-Cox t | 31.5771159 |
S’ha tingut en compte per un costat les interaccions entre els termes existents dels models simples optimitzats i, excepte en el cas de reducció per colinealitat, s’ha comprobat l’efecte de l’aparició dels factors d’interacció per paràmetres que no estaven en el model òptim simple en els de variable de resposta transformada. A continuació es resumeix el model trobat com a òptim al fer els càlculs i la p-valor del model òptim d’interacció amb el model simple òptim del qual parteixen:
S’accepta per tant de moment els models òptims amb la presència d’interaccions i es torna a generar alguns dels diagnòstics per observar les diferències:
Lake1Int.SM0.sdres <- rstandard(intM0)
Lake1Int.SM0.studres <- studres(intM0)
par(mfrow=c(2,2))
qqnorm(Lake1.SM0.sdres,
ylab="Standardized Residuals",
xlab="Normal Scores",
main="M. Original - Standardized Residuals")
qqline(Lake1.SM0.sdres)
qqnorm(Lake1.SM0.studres,
ylab="Studentized Residuals",
xlab="Normal Scores",
main="M. Original - Studentized Residuals")
qqline(Lake1.SM0.studres)
qqnorm(Lake1Int.SM0.sdres,
ylab="Standardized Residuals",
xlab="Normal Scores",
main="M. Original Int. - Standardized Residuals")
qqline(Lake1Int.SM0.sdres)
qqnorm(Lake1Int.SM0.studres,
ylab="Studentized Residuals",
xlab="Normal Scores",
main="M. Original Int. - Studentized Residuals")
qqline(Lake1Int.SM0.studres)
S0
##
## Shapiro-Wilk normality test
##
## data: Lake1c3.step.model0$residuals
## W = 0.87288, p-value = 0.0005645
(S0Int <- shapiro.test(intM0$residuals))
##
## Shapiro-Wilk normality test
##
## data: intM0$residuals
## W = 0.92773, p-value = 0.01359
Lake1.SM0col.sdres <- rstandard(Lake1.step.model0col)
Lake1.SM0col.studres <- studres(Lake1.step.model0col)
Lake1Int.SM0col.sdres <- rstandard(intM0col)
Lake1Int.SM0col.studres <- studres(intM0col)
par(mfrow=c(2,2))
qqnorm(Lake1.SM0col.sdres,
ylab="Standardized Residuals",
xlab="Normal Scores",
main="M. Original Correcció Col. - Standardized Residuals")
qqline(Lake1.SM0.sdres)
qqnorm(Lake1.SM0col.studres,
ylab="Studentized Residuals",
xlab="Normal Scores",
main="M. Original Correcció Col. - Studentized Residuals")
qqline(Lake1.SM0.studres)
qqnorm(Lake1Int.SM0col.sdres,
ylab="Standardized Residuals",
xlab="Normal Scores",
main="M. Original Int. Correcció Col. - Standardized Residuals")
qqline(Lake1Int.SM0col.sdres)
qqnorm(Lake1Int.SM0col.studres,
ylab="Studentized Residuals",
xlab="Normal Scores",
main="M. Original Int. Correcció Col. - Studentized Residuals")
qqline(Lake1Int.SM0col.studres)
(S0col <- shapiro.test(Lake1.step.model0col$residuals))
##
## Shapiro-Wilk normality test
##
## data: Lake1.step.model0col$residuals
## W = 0.88007, p-value = 0.0005286
(S0colInt <- shapiro.test(intM0col$residuals))
##
## Shapiro-Wilk normality test
##
## data: intM0col$residuals
## W = 0.88007, p-value = 0.0005286
Lake1Int.SM1.sdres <- rstandard(intM1)
Lake1Int.SM1.studres <- studres(intM1)
par(mfrow=c(2,2))
qqnorm(Lake1.SM1.sdres,
ylab="Standardized Residuals",
xlab="Normal Scores",
main="M. Log - Standardized Residuals")
qqline(Lake1.SM1.sdres)
qqnorm(Lake1.SM1.studres,
ylab="Studentized Residuals",
xlab="Normal Scores",
main="M. Log - Studentized Residuals")
qqline(Lake1.SM1.studres)
qqnorm(Lake1Int.SM1.sdres,
ylab="Standardized Residuals",
xlab="Normal Scores",
main="M. Log Int.- Standardized Residuals")
qqline(Lake1Int.SM1.sdres)
qqnorm(Lake1Int.SM1.studres,
ylab="Studentized Residuals",
xlab="Normal Scores",
main="M. Log Int.- Studentized Residuals")
qqline(Lake1Int.SM1.studres)
S1
##
## Shapiro-Wilk normality test
##
## data: Lake1.step.model1$residuals
## W = 0.94994, p-value = 0.07547
(S1Int <- shapiro.test(intM1$residuals))
##
## Shapiro-Wilk normality test
##
## data: intM1$residuals
## W = 0.9557, p-value = 0.1192
Lake1Int.SM2.sdres <- rstandard(intM2)
Lake1Int.SM2.studres <- studres(intM2)
par(mfrow=c(2,2))
qqnorm(Lake1.SM2.sdres,
ylab="Standardized Residuals",
xlab="Normal Scores",
main="M. Box-Cox - Standardized Residuals")
qqline(Lake1.SM2.sdres)
qqnorm(Lake1.SM2.studres,
ylab="Studentized Residuals",
xlab="Normal Scores",
main="M. Box-Cox - Studentized Residuals")
qqline(Lake1.SM2.studres)
qqnorm(Lake1Int.SM2.sdres,
ylab="Standardized Residuals",
xlab="Normal Scores",
main="M. Box-Cox Int.- Standardized Residuals")
qqline(Lake1Int.SM2.sdres)
qqnorm(Lake1Int.SM2.studres,
ylab="Studentized Residuals",
xlab="Normal Scores",
main="M. Box-Cox Int.- Studentized Residuals")
qqline(Lake1Int.SM2.studres)
S2
##
## Shapiro-Wilk normality test
##
## data: Lake1.step.model2$residuals
## W = 0.91901, p-value = 0.00717
(S2Int <- shapiro.test(intM2$residuals))
##
## Shapiro-Wilk normality test
##
## data: intM2$residuals
## W = 0.96332, p-value = 0.2172
Tot i que en general sembla millorar la normalitat pel test de Shapiro en els gràfics queda clar pel model original (tant complet com corregit per la colinealitat), però no tant aparent en els models transformats.
pHInt01 <- ggplot(LAKE1, aes(x=fitted(intM0), y=residuals(intM0),
shape=Grupo, color=Grupo)) +
geom_jitter(position=position_jitter(1), cex=1.5) +
geom_hline(yintercept=0, linetype="dashed", color = "black") +
labs(title="Fitted vs Residuals (Int.Homocedasticity test)",
x="Fitted", y = "Residuals")
pHInt02 <-ggplot(LAKE1, aes(x=fitted(intM0), y=sqrt(abs(residuals(intM0))),
shape=Grupo, color=Grupo)) +
geom_jitter(position=position_jitter(1), cex=1.5) +
geom_hline(yintercept=0, linetype="dashed", color = "black") +
labs(title="Fitted vs Residuals (Int.Homocedasticity test)",
x="Fitted", y = expression(sqrt("|residuals|")))
grid.arrange(pH01, pH02, pHInt01, pHInt02,ncol=2)
par(mfrow=c(1,2))
plot(Lake1.step.model0,which=1)
plot(intM0,which=1)
BP01
##
## studentized Breusch-Pagan test
##
## data: Lake1.step.model0
## BP = 17.345, df = 3, p-value = 6e-04
(BPInt01<-bptest(intM0))
##
## studentized Breusch-Pagan test
##
## data: intM0
## BP = 8.9796, df = 7, p-value = 0.2541
pHcol01 <- ggplot(LAKE1, aes(x=fitted(Lake1.step.model0col), y=residuals(Lake1.step.model0col),
shape=Grupo, color=Grupo)) +
geom_jitter(position=position_jitter(1), cex=1.5) +
geom_hline(yintercept=0, linetype="dashed", color = "black") +
labs(title="Fitted vs Residuals (Col.Homocedasticity test)",
x="Fitted", y = "Residuals")
pHcol02 <-ggplot(LAKE1, aes(x=fitted(Lake1.step.model0col), y=sqrt(abs(residuals(Lake1.step.model0col))),
shape=Grupo, color=Grupo)) +
geom_jitter(position=position_jitter(1), cex=1.5) +
geom_hline(yintercept=0, linetype="dashed", color = "black") +
labs(title="Fitted vs Residuals (Col.Homocedasticity test)",
x="Fitted", y = expression(sqrt("|residuals|")))
pHcolInt01 <- ggplot(LAKE1, aes(x=fitted(intM0col), y=residuals(intM0col),
shape=Grupo, color=Grupo)) +
geom_jitter(position=position_jitter(1), cex=1.5) +
geom_hline(yintercept=0, linetype="dashed", color = "black") +
labs(title="Fitted vs Residuals (colInt.Homocedasticity test)",
x="Fitted", y = "Residuals")
pHcolInt02 <-ggplot(LAKE1, aes(x=fitted(intM0col), y=sqrt(abs(residuals(intM0col))),
shape=Grupo, color=Grupo)) +
geom_jitter(position=position_jitter(1), cex=1.5) +
geom_hline(yintercept=0, linetype="dashed", color = "black") +
labs(title="Fitted vs Residuals (colInt.Homocedasticity test)",
x="Fitted", y = expression(sqrt("|residuals|")))
grid.arrange(pHcol01, pHcol02, pHcolInt01, pHcolInt02,ncol=2)
par(mfrow=c(1,2))
plot(Lake1.step.model0col,which=1)
plot(intM0col,which=1)
(BPcol01<-bptest(Lake1.step.model0col))
##
## studentized Breusch-Pagan test
##
## data: Lake1.step.model0col
## BP = 12.423, df = 2, p-value = 0.002006
(BPcolInt01<-bptest(intM0col))
##
## studentized Breusch-Pagan test
##
## data: intM0col
## BP = 12.423, df = 2, p-value = 0.002006
pHInt11 <- ggplot(LAKE1, aes(x=fitted(intM1), y=residuals(intM1),
shape=Grupo, color=Grupo)) +
geom_jitter(position=position_jitter(1), cex=1.5) +
geom_hline(yintercept=0, linetype="dashed", color = "black") +
labs(title="Fitted vs Residuals (Log.Int.Hom. test)",
x="Fitted", y = "Residuals")
pHInt12 <-ggplot(LAKE1, aes(x=fitted(intM1), y=sqrt(abs(residuals(intM1))),
shape=Grupo, color=Grupo)) +
geom_jitter(position=position_jitter(1), cex=1.5) +
geom_hline(yintercept=0, linetype="dashed", color = "black") +
labs(title="Fitted vs Residuals (Log.Int.Hom. test)",
x="Fitted", y = expression(sqrt("|residuals|")))
grid.arrange(pH11, pH12,pHInt11,pHInt12,ncol=2)
plot(Lake1.step.model1,which=1)
plot(intM1,which=1)
BP11
##
## studentized Breusch-Pagan test
##
## data: Lake1.step.model1
## BP = 19.474, df = 1, p-value = 1.02e-05
(BPInt11<-bptest(intM1))
##
## studentized Breusch-Pagan test
##
## data: intM1
## BP = 10.752, df = 3, p-value = 0.01314
pHInt21 <- ggplot(LAKE1, aes(x=fitted(intM2), y=residuals(intM2),
shape=Grupo, color=Grupo)) +
geom_jitter(position=position_jitter(1), cex=1.5) +
geom_hline(yintercept=0, linetype="dashed", color = "black") +
labs(title="Fitted vs Residuals (Int.BC.Hom test)",
x="Fitted", y = "Residuals")
pHInt22 <-ggplot(LAKE1, aes(x=fitted(intM2), y=sqrt(abs(residuals(intM2))),
shape=Grupo, color=Grupo)) +
geom_jitter(position=position_jitter(1), cex=1.5) +
geom_hline(yintercept=0, linetype="dashed", color = "black") +
labs(title="Fitted vs Residuals (Int.BC.Hom test)",
x="Fitted", y = expression(sqrt("|residuals|")))
grid.arrange(pH21, pH22,pHInt21,pHInt22,ncol=2)
plot(Lake1.step.model2,which=1)
plot(intM2,which=1)
BP21
##
## studentized Breusch-Pagan test
##
## data: Lake1.step.model2
## BP = 5.718, df = 1, p-value = 0.01679
(BPint21<-bptest(intM2))
##
## studentized Breusch-Pagan test
##
## data: intM2
## BP = 5.8317, df = 3, p-value = 0.1201
En els dos models sense transformació, independentment del test d’homogeneïtat els gràfics segueixen mostrant un patró que no correspondria a la homocedasticitat de residus i que podria ser preocupant per assumir aquesta condició. En els altres dos models que s’havien transformat, el cas logarítmic millora lleugerament tot o que encara no es pot assegurar veient els gràfics que hi hagi la condició, i el cas del transformat de box cox potser és l’únic cas on no es veuria tan clarament l’heteroscedasticitat i inclús milloraria en els models amb els termes d’interacció.
En resum els models generats per les proves de diagnòstic encara poden patir les conseqüències de no cumplir les 2 assumpcions principals en els resultats d’estimació i altres contrastos d’hipòtesis, pel que s’ha d’anar en compte al fer-los servir.
#Es simplifica els noms dels models inicials per poder utilitzar-los més comodament en el codi
M0<-Lake1.step.model0
M0col<-Lake1.step.model0col
M1<-Lake1.step.model1
M2<-Lake1.step.model2
M0s<-Lake1.step.model0s
M0cols<-Lake1.step.model0col.s
M1s<-Lake1.step.model1s
M2s<-Lake1.step.model2s
Pels resultats que s’han anat obtinguent semblaria que es tracta en una possible variant de la simulació ja vista identificada com a He2 amb paral·lelisme i diferències entre individus, tot i que en aquest cas és possible que sigui un model més complex al afegir la heteroscedasticitat que s’ha anat comprovant i els efectes addicionals en el cas del model original.
Si s’analitza amb més detall el context en el que es troba el model s’ha de considerar que el conjunt utilitzat és un subgrup dels individus del model original, que a la vegada ja es una subpoblació de la població en general que pugui tenir el VIH. Cal tenir e compte que l’objectiu d’aquest model de predicció és fer inferència en la població general d’afectats de VIH i no només les de la mostra estudiada, i per tal de fer això s’ha de tenir en compte que si es realitzessin més experiments amb altres conjunts d’individus la resposta podria variar. Per poder representar això és necessari afegir el component aleatori que aporta l’efecte individu amb la seva pròpia variància dins el model.
Per fer aquesta actualització del model s’optarà com a les simulacions per tenir en compte 3 models amb diferents plantejaments del càlcul òptim: El model clàssic de mínims quadrats amb l’error aleatori, el model de màxima versemblança ML i el model de la màxima versemblança restringida REML. Així en els següents models es veurà si canvia molt la bondat d’ajust respecte als models amb interaccions afegint els termes de efectes aleatoris i tinguent en compte els 2 sistemes addicionals d’estimació de resultats òptims.
Seguint un model semblant a com s’ha procedit en el codi de programació de la simulació es genera una funció per facilitar la obtenció de resultats de diferents conjunts en cas necessari:
Update_LM_RI <- function(LM,RI_var,df) { #Objecte del tipus LM
fmOLS <- update(LM$call$formula,paste(" ~ . + Error(",RI_var,")"))
ML_I_RI_sOLS <- summary(ML_I_RI_OLS <- aov(fmOLS, data = df))
fmMLf <- LM$call$formula
fmMLr <- formula(paste("~1|",RI_var))
ML_I_RI_sMLE <- summary(ML_I_RI_MLE <- lme(fmMLf, random = fmMLr, data = df, method="ML"))
ML_I_RI_sREML <- summary(ML_I_RI_REML <- lme(fmMLf, random = fmMLr, data = df, 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 (Within)","MLE I RI","REML I RI",
"OLS I RI base","MLE I RI base","REML I RI base")
return(l)
}
M0_RI<-Update_LM_RI(M0,"Ident",LAKE1)
eMLE_M0_RI<-exactRLRT(M0_RI_MLE<-lme(M0$call$formula, random = ~1|Ident, data = LAKE1, 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
eREML_M0_RI<-exactRLRT(M0_RI_REML<-lme(M0$call$formula, random = ~1|Ident, data = LAKE1, method="REML"))
## Warning in model.matrix.default(~m$groups[[n.levels - i + 1]] - 1,
## contrasts.arg = c("contr.treatment", : non-list contrasts argument ignored
intM0_RI<-Update_LM_RI(intM0,"Ident",LAKE1)
eMLE_intM0_RI<-exactRLRT(intM0_RI_MLE<-lme(intM0$call$formula, random = ~1|Ident, data = LAKE1, 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
eREML_intM0_RI<-exactRLRT(intM0_RI_REML<-lme(intM0$call$formula, random = ~1|Ident, data = LAKE1, method="REML"))
## Warning in model.matrix.default(~m$groups[[n.levels - i + 1]] - 1,
## contrasts.arg = c("contr.treatment", : non-list contrasts argument ignored
M0col_RI<-Update_LM_RI(M0col,"Ident",LAKE1)
eMLE_M0col_RI<-exactRLRT(M0col_RI_MLE<-lme(M0col$call$formula, random = ~1|Ident, data = LAKE1, 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
eREML_M0col_RI<-exactRLRT(M0col_RI_REML<-lme(M0col$call$formula, random = ~1|Ident, data = LAKE1, method="REML"))
## Warning in model.matrix.default(~m$groups[[n.levels - i + 1]] - 1,
## contrasts.arg = c("contr.treatment", : non-list contrasts argument ignored
intM0col_RI<-Update_LM_RI(intM0col,"Ident",LAKE1)
eMLE_intM0col_RI<-exactRLRT(intM0col_RI_MLE<-lme(intM0col$call$formula, random = ~1|Ident, data = LAKE1, 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
eREML_intM0col_RI<-exactRLRT(intM0col_RI_REML<-lme(intM0col$call$formula, random = ~1|Ident, data = LAKE1, method="REML"))
## Warning in model.matrix.default(~m$groups[[n.levels - i + 1]] - 1,
## contrasts.arg = c("contr.treatment", : non-list contrasts argument ignored
M1_RI<-Update_LM_RI(M1,"Ident",LAKE1)
eMLE_M1_RI<-exactRLRT(M1_RI_MLE<-lme(M1$call$formula, random = ~1|Ident, data = LAKE1, 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
eREML_M1_RI<-exactRLRT(M1_RI_REML<-lme(M1$call$formula, random = ~1|Ident, data = LAKE1, method="REML"))
## Warning in model.matrix.default(~m$groups[[n.levels - i + 1]] - 1,
## contrasts.arg = c("contr.treatment", : non-list contrasts argument ignored
intM1_RI<-Update_LM_RI(intM1,"Ident",LAKE1)
eMLE_intM1_RI<-exactRLRT(intM1_RI_MLE<-lme(intM1$call$formula, random = ~1|Ident, data = LAKE1, 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
eREML_intM1_RI<-exactRLRT(intM1_RI_REML<-lme(intM1$call$formula, random = ~1|Ident, data = LAKE1, method="REML"))
## Warning in model.matrix.default(~m$groups[[n.levels - i + 1]] - 1,
## contrasts.arg = c("contr.treatment", : non-list contrasts argument ignored
M2_RI<-Update_LM_RI(M2,"Ident",LAKE1)
eMLE_M2_RI<-exactRLRT(M2_RI_MLE<-lme(M2$call$formula, random = ~1|Ident, data = LAKE1, 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
eREML_M2_RI<-exactRLRT(M2_RI_REML<-lme(M2$call$formula, random = ~1|Ident, data = LAKE1, method="REML"))
## Warning in model.matrix.default(~m$groups[[n.levels - i + 1]] - 1,
## contrasts.arg = c("contr.treatment", : non-list contrasts argument ignored
intM2_RI<-Update_LM_RI(intM2,"Ident",LAKE1)
eMLE_intM2_RI<-exactRLRT(intM2_RI_MLE<-lme(intM2$call$formula, random = ~1|Ident, data = LAKE1, 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
eREML_intM2_RI<-exactRLRT(intM2_RI_REML<-lme(intM2$call$formula, random = ~1|Ident, data = LAKE1, method="REML"))
## Warning in model.matrix.default(~m$groups[[n.levels - i + 1]] - 1,
## contrasts.arg = c("contr.treatment", : non-list contrasts argument ignored
NomsModels <- c("M.Base","M.Base Int.",
"M.Corr.Col","M.Corr.Col. Int.",
"M.LogTrans.","M.LogTrans.Int.",
"M.Box-Cox","M.Box-Cox Int.")
RLRT_resum_RI <- data.frame(Model_Lineal = paste(rep(NomsModels,each=2),c("MLE","REML")))
RLRT_list_RI <- list(eMLE_M0_RI,eREML_M0_RI,eMLE_intM0_RI,eREML_intM0_RI,
eMLE_M0col_RI,eREML_M0col_RI,eMLE_intM0col_RI,eREML_intM0col_RI,
eMLE_M1_RI,eREML_M1_RI,eMLE_intM1_RI,eREML_intM1_RI,
eMLE_M2_RI,eREML_M2_RI,eMLE_intM2_RI,eREML_intM2_RI)
RLRT_resum_RI$"p-valor RLRT" <- sapply(1:length(RLRT_list_RI),function(l1){
if (class(RLRT_list_RI[[l1]])!="htest"){
return("-")
} else {
round(RLRT_list_RI[[l1]]$p.value,4)
}
})
kable(RLRT_resum_RI, booktabs = T) %>% kable_styling(latex_options = c("striped","condensed"))
Model_Lineal | p-valor RLRT |
---|---|
M.Base MLE | 1 |
M.Base REML | 1 |
M.Base Int. MLE | 1 |
M.Base Int. REML | 1 |
M.Corr.Col MLE | 1 |
M.Corr.Col REML | 1 |
M.Corr.Col. Int. MLE | 1 |
M.Corr.Col. Int. REML | 1 |
M.LogTrans. MLE | 1 |
M.LogTrans. REML | 1 |
M.LogTrans.Int. MLE | 1 |
M.LogTrans.Int. REML | 1 |
M.Box-Cox MLE | 1 |
M.Box-Cox REML | 1 |
M.Box-Cox Int. MLE | 1 |
M.Box-Cox Int. REML | 1 |
Tot i que la significació de l’efecte aleatori RI surt no significatiu a nivell més concret pot ser que signifiquin algun canvi.
Tot i tenir un resultat de no significació en tots els models amb l’efecte aleatori de l’intercepció depenent de l’individu, es recupera la funció ja utilitzada per evaluar alguns indicadors 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 = "Error.std. residual (Sigma)")) +
theme(axis.text.x = element_text(angle=80,size = 10,hjust=0))
return(list(R1,R2))
}
M_RI_listA0 <- list(M0,intM0,M0col,intM0col,
M0_RI[1:3],intM0_RI[1:3],M0col_RI[1:3],intM0col_RI[1:3])
M_RI_listA <- ConvList(M_RI_listA0)
## Warning in if (class(List[[l1]]) == "list") {: the condition has length > 1
## and only the first element will be used
## Warning in if (class(List[[l1]]) == "list") {: the condition has length > 1
## and only the first element will be used
## Warning in if (class(List[[l1]]) == "list") {: the condition has length > 1
## and only the first element will be used
## Warning in if (class(List[[l1]]) == "list") {: the condition has length > 1
## and only the first element will be used
M_RI_listB0 <- list(M1,intM1,M2,intM2,
M1_RI[1:3],intM1_RI[1:3],M2_RI[1:3],intM2_RI[1:3])
M_RI_listB <- ConvList(M_RI_listB0)
## Warning in if (class(List[[l1]]) == "list") {: the condition has length > 1
## and only the first element will be used
## Warning in if (class(List[[l1]]) == "list") {: the condition has length > 1
## and only the first element will be used
## Warning in if (class(List[[l1]]) == "list") {: the condition has length > 1
## and only the first element will be used
## Warning in if (class(List[[l1]]) == "list") {: the condition has length > 1
## and only the first element will be used
names(M_RI_listA)[1:4] <- NomsModels[1:4]
names(M_RI_listA)[(4+1):length(M_RI_listA)]<-paste(rep(NomsModels[1:4],each=3),rep(paste("RI",c("OLS","MLE","REML")),4))
names(M_RI_listB)[1:4] <- NomsModels[5:8]
names(M_RI_listB)[(4+1):length(M_RI_listB)]<-paste(rep(NomsModels[5:8],each=3),rep(paste("RI",c("OLS","MLE","REML")),4))
(DiagA_RI <- ML_Resum(M_RI_listA,"Lake1"))
## [[1]]
## Sigma AIC BIC
## Lake1--M.Base 81386.03 1023.8574 1032.3018
## Lake1--M.Base Int. 47433.69 983.9564 999.1563
## Lake1--M.Corr.Col 84579.43 1026.0323 1032.7878
## Lake1--M.Corr.Col. Int. 84579.43 1026.0323 1032.7878
## Lake1--M.Base RI OLS 73964.95 813.1882 819.0511
## Lake1--M.Base RI MLE 77209.57 1025.8574 1035.9907
## Lake1--M.Base RI REML 81386.03 954.0256 963.5267
## Lake1--M.Base Int. RI OLS 50420.37 791.9142 803.6401
## Lake1--M.Base Int. RI MLE 42425.99 985.9564 1002.8452
## Lake1--M.Base Int. RI REML 47433.69 873.8382 888.4956
## Lake1--M.Corr.Col RI OLS 91010.46 825.5455 829.9427
## Lake1--M.Corr.Col RI MLE 81345.89 1028.0323 1036.4767
## Lake1--M.Corr.Col RI REML 84579.43 973.3637 981.4183
## Lake1--M.Corr.Col. Int. RI OLS 91010.46 825.5455 829.9427
## Lake1--M.Corr.Col. Int. RI MLE 81345.89 1028.0323 1036.4767
## Lake1--M.Corr.Col. Int. RI REML 84579.43 973.3637 981.4183
##
## [[2]]
BestA <- unique(c(order(DiagA_RI[[1]]$AIC,decreasing=F)[1:4],
order(DiagA_RI[[1]]$Sigma,decreasing=F)[1:4]))
ML_Resum(M_RI_listA[sort(BestA)],"Lake1") #Comparativa més selectiva
## [[1]]
## Sigma AIC BIC
## Lake1--M.Base Int. 47433.69 983.9564 999.1563
## Lake1--M.Base RI OLS 73964.95 813.1882 819.0511
## Lake1--M.Base Int. RI OLS 50420.37 791.9142 803.6401
## Lake1--M.Base Int. RI MLE 42425.99 985.9564 1002.8452
## Lake1--M.Base Int. RI REML 47433.69 873.8382 888.4956
## Lake1--M.Corr.Col RI OLS 91010.46 825.5455 829.9427
## Lake1--M.Corr.Col. Int. RI OLS 91010.46 825.5455 829.9427
##
## [[2]]
(DiagB_RI<-ML_Resum(M_RI_listB,"Lake1"))
## [[1]]
## Sigma AIC BIC
## Lake1--M.LogTrans. 2.1784 179.7516 184.8182
## Lake1--M.LogTrans.Int. 1.8898 170.2175 178.6619
## Lake1--M.Box-Cox 296.4675 572.8184 577.8850
## Lake1--M.Box-Cox Int. 277.9218 569.4879 577.9322
## Lake1--M.LogTrans. RI OLS 2.3434 148.2996 151.2311
## Lake1--M.LogTrans. RI MLE 2.1233 181.7516 188.5071
## Lake1--M.LogTrans. RI REML 2.1784 188.0539 194.6042
## Lake1--M.LogTrans.Int. RI OLS 1.9420 138.1412 144.0041
## Lake1--M.LogTrans.Int. RI MLE 1.7928 172.2175 182.3508
## Lake1--M.LogTrans.Int. RI REML 1.8898 192.9877 202.4888
## Lake1--M.Box-Cox RI OLS 305.9566 460.0965 463.0279
## Lake1--M.Box-Cox RI MLE 288.9608 574.8184 581.5739
## Lake1--M.Box-Cox RI REML 296.4675 561.4673 568.0177
## Lake1--M.Box-Cox Int. RI OLS 283.5571 457.0965 462.9594
## Lake1--M.Box-Cox Int. RI MLE 263.6598 571.4879 581.6211
## Lake1--M.Box-Cox Int. RI REML 277.9218 552.3310 561.8321
##
## [[2]]
BestB <- unique(c(order(DiagB_RI[[1]]$AIC,decreasing=F)[1:4],
order(DiagB_RI[[1]]$Sigma,decreasing=F)[1:4]))
ML_Resum(M_RI_listB[sort(BestB)],"Lake1") #Comparativa més selectiva
## [[1]]
## Sigma AIC BIC
## Lake1--M.LogTrans.Int. 1.8898 170.2175 178.6619
## Lake1--M.LogTrans. RI OLS 2.3434 148.2996 151.2311
## Lake1--M.LogTrans.Int. RI OLS 1.9420 138.1412 144.0041
## Lake1--M.LogTrans.Int. RI MLE 1.7928 172.2175 182.3508
## Lake1--M.LogTrans.Int. RI REML 1.8898 192.9877 202.4888
##
## [[2]]
#Deixem només de base els d'interacció i fora el de correcció de colinealitat també:
(DiagAi_RI <- ML_Resum(M_RI_listA[-c(1,3:7,11:16)],"Lake1"))
## [[1]]
## Sigma AIC BIC
## Lake1--M.Base Int. 47433.69 983.9564 999.1563
## Lake1--M.Base Int. RI OLS 50420.37 791.9142 803.6401
## Lake1--M.Base Int. RI MLE 42425.99 985.9564 1002.8452
## Lake1--M.Base Int. RI REML 47433.69 873.8382 888.4956
##
## [[2]]
BestAi <- unique(c(order(DiagAi_RI[[1]]$AIC,decreasing=F)[1:4],
order(DiagAi_RI[[1]]$Sigma,decreasing=F)[1:4]))
ML_Resum((M_RI_listA[-c(1,3:7,11:16)])[sort(BestAi)],"Lake1") #Comparativa més selectiva
## [[1]]
## Sigma AIC BIC
## Lake1--M.Base Int. 47433.69 983.9564 999.1563
## Lake1--M.Base Int. RI OLS 50420.37 791.9142 803.6401
## Lake1--M.Base Int. RI MLE 42425.99 985.9564 1002.8452
## Lake1--M.Base Int. RI REML 47433.69 873.8382 888.4956
##
## [[2]]
(DiagBi_RI<-ML_Resum(M_RI_listB[-c(1,3,5:7,11:13)],"Lake1"))
## [[1]]
## Sigma AIC BIC
## Lake1--M.LogTrans.Int. 1.8898 170.2175 178.6619
## Lake1--M.Box-Cox Int. 277.9218 569.4879 577.9322
## Lake1--M.LogTrans.Int. RI OLS 1.9420 138.1412 144.0041
## Lake1--M.LogTrans.Int. RI MLE 1.7928 172.2175 182.3508
## Lake1--M.LogTrans.Int. RI REML 1.8898 192.9877 202.4888
## Lake1--M.Box-Cox Int. RI OLS 283.5571 457.0965 462.9594
## Lake1--M.Box-Cox Int. RI MLE 263.6598 571.4879 581.6211
## Lake1--M.Box-Cox Int. RI REML 277.9218 552.3310 561.8321
##
## [[2]]
BestBi <- unique(c(order(DiagBi_RI[[1]]$AIC,decreasing=F)[1:4],
order(DiagBi_RI[[1]]$Sigma,decreasing=F)[1:4]))
ML_Resum((M_RI_listB[-c(1,3,5:7,11:13)])[sort(BestBi)],"Lake1") #Comparativa més selectiva
## [[1]]
## Sigma AIC BIC
## Lake1--M.LogTrans.Int. 1.8898 170.2175 178.6619
## Lake1--M.LogTrans.Int. RI OLS 1.9420 138.1412 144.0041
## Lake1--M.LogTrans.Int. RI MLE 1.7928 172.2175 182.3508
## Lake1--M.LogTrans.Int. RI REML 1.8898 192.9877 202.4888
##
## [[2]]
En els primers gràfics s’ha comparat tots els models, per una banda sense transformació i per l’altra amb les transformacions. Posteriorment s’ha fet una selecció dels models amb els indicadors AIC/BIC i la variància residual més baixos per veure quins eren els millors models a tenir en compte:
Es calcula el coeficient de determinació \(R^2\). En el cas dels models amb efectes aleatoris es calcula a partir del coeficient de determinació condicionat que suma el que aporta els efectes fixos i els efectes aleatoris (en el cas d’aquest models només s’ha trobat fàcilment com calcular-ho en R amb els models calculats per ML i dins d’aquí en els casos amb termes d’interacció es fa una aproximació mitjançant la correlació lineal entre la resposta i els valors ajustats):
M_RI_listA_R2 <- list(M0,intM0,M0col,intM0col,M0_RI_MLE,M0_RI_REML,
intM0_RI_MLE,intM0_RI_REML,M0col_RI_MLE,M0col_RI_REML,
intM0col_RI_MLE,intM0col_RI_REML)
names(M_RI_listA_R2)[1:4] <- NomsModels[1:4]
names(M_RI_listA_R2)[(4+1):length(M_RI_listA_R2)]<-paste(rep(NomsModels[1:4],each=2),rep(paste("RI",c("MLE","REML")),4))
M_RI_dfA_R2 <- data.frame(Model = names(M_RI_listA_R2))
M_RI_dfA_R2$R2 <- sapply(1:length(M_RI_listA_R2),function(M){
if (class(M_RI_listA_R2[[M]])[1]=="lmm"){
summary(M_RI_listA_R2[[M]])$r.squared
} else if (class(M_RI_listA_R2[[M]])=="lme"){
if(all(paste(attr(M_RI_listA_R2[[M]]$terms,"predvars")[-c(1:2)]) ==
paste(attr(M_RI_listA_R2[[M]]$terms,"term.labels")))){
rsquared(M_RI_listA_R2[[M]])$Conditional
} else {
summary(lm(formula(paste(paste("M_RI_listA_R2[[",M,"]]$data$",
attr(M_RI_listA_R2[[M]]$terms,"variables")[[2]],sep=""),
"~","M_RI_listA_R2[[",M,"]]$fitted"),sep="")))$r.squared
}
}
})
## Warning in paste(attr(M_RI_listA_R2[[M]]$terms, "predvars")[-c(1:2)]) == :
## longer object length is not a multiple of shorter object length
## Warning in paste(attr(M_RI_listA_R2[[M]]$terms, "predvars")[-c(1:2)]) == :
## longer object length is not a multiple of shorter object length
R2_A <- ggplot(data=M_RI_dfA_R2, aes(x=Model, y=R2, fill=Model)) +
geom_bar(stat="identity", position=position_dodge())+
scale_fill_brewer(palette="Set3") + theme_minimal() +
labs(y = "R2 Coeff.") +
theme(axis.text.x = element_text(angle=80,size = 10,hjust=0))
M_RI_listB_R2 <- list(M1,intM1,M2,intM2,M1_RI_MLE,
M1_RI_REML,intM1_RI_MLE,intM1_RI_REML,
M2_RI_MLE,M2_RI_REML,intM2_RI_MLE,intM2_RI_REML)
names(M_RI_listB_R2)[1:4] <- NomsModels[5:8]
names(M_RI_listB_R2)[(4+1):length(M_RI_listB_R2)]<-paste(rep(NomsModels[5:8],each=2),rep(paste("RI",c("MLE","REML")),4))
M_RI_dfB_R2 <- data.frame(Model = names(M_RI_listB_R2))
M_RI_dfB_R2 <- data.frame(Model = names(M_RI_listB_R2))
M_RI_dfB_R2$R2 <- sapply(1:length(M_RI_listB_R2),function(M){
if (class(M_RI_listB_R2[[M]])[1]=="lmm"){
summary(M_RI_listB_R2[[M]])$r.squared
} else if (class(M_RI_listB_R2[[M]])=="lme"){
if(all(paste(attr(M_RI_listB_R2[[M]]$terms,"predvars")[-c(1:2)]) ==
paste(attr(M_RI_listB_R2[[M]]$terms,"term.labels")))){
rsquared(M_RI_listB_R2[[M]])$Conditional
} else {
summary(lm(formula(paste(paste("M_RI_listB_R2[[",M,"]]$data$",
attr(M_RI_listB_R2[[M]]$terms,"variables")[[2]],sep=""),
"~","M_RI_listB_R2[[",M,"]]$fitted"),sep="")))$r.squared
}
}
})
## Warning in paste(attr(M_RI_listB_R2[[M]]$terms, "predvars")[-c(1:2)]) == :
## longer object length is not a multiple of shorter object length
## Warning in paste(attr(M_RI_listB_R2[[M]]$terms, "predvars")[-c(1:2)]) == :
## longer object length is not a multiple of shorter object length
## Warning in paste(attr(M_RI_listB_R2[[M]]$terms, "predvars")[-c(1:2)]) == :
## longer object length is not a multiple of shorter object length
## Warning in paste(attr(M_RI_listB_R2[[M]]$terms, "predvars")[-c(1:2)]) == :
## longer object length is not a multiple of shorter object length
R2_B <- ggplot(data=M_RI_dfB_R2, aes(x=Model, y=R2, fill=Model)) +
geom_bar(stat="identity", position=position_dodge())+
scale_fill_brewer(palette="Set3") + theme_minimal() +
labs(y = "R2 Coeff.") +
theme(axis.text.x = element_text(angle=80,size = 10,hjust=0))
R2_A
R2_B
BestA_R2 <- order(M_RI_dfA_R2$R2,decreasing=T)[1:4]
BestB_R2 <- order(M_RI_dfB_R2$R2,decreasing=T)[1:4]
M_RI_df_sel_R2 <- rbind(M_RI_dfA_R2[BestA_R2,],M_RI_dfB_R2[BestB_R2,])
R2_sel <- ggplot(data=M_RI_df_sel_R2, aes(x=Model, y=R2, fill=Model)) +
geom_bar(stat="identity", position=position_dodge())+
scale_fill_brewer(palette="Set3") + theme_minimal() +
labs(y = "R2 Coeff.") +
theme(axis.text.x = element_text(angle=80,size = 10,hjust=0))
R2_sel
Com en l’anterior cas s’ha fet dos gràfics amb tots els models i un tercer gràfic on s’ha seleccionat les \(R^2\) més elevades. En aquest cas es pot veure que s’han seleccionat el model base amb interaccions i el Box-Cox transformat amb interaccions que també dona informació de quins models tenen un % d’explicabilitat més elevat de la resposta, però tampoc sembla haver diferències amb els models on s’afegeix l’efecte aleatori RI.
Tot i que en el test realitzat d’interacció entre factor i covariant ha donat el resultat del paral·lelisme de pendent entre individus, seria interessant veure com es comporta si s’afegeix en el model una component aleatòria d’aquest paràmetre ja que en el supòsit que el test de interacció fet no fos correcte la interpretació, igual que la resposta dels individus, la degradació del paràmetre observat en front al temps podria no ser igual per tots els individus. Al fer la inferència a partir de la mostra podria ser que es trobés que hi ha una degradació equivalent per tots els individus o al contrari que cada un es comporta diferent. Això donaria peu també a suposar un possible efecte aleatori en la pendent respecte a la covariant temps de cada individu i per tant també amb un altre efecte de variància dins el model, encara que pels resultats obtinguts no s’espera que sigui així.
Update_LM_RIS <- function(LM,RI_var,RIS_var,df) { #Objecte del tipus LM
fmOLS <- update(LM$call$formula,paste(" ~ . + Error(",RI_var,"/",RIS_var,")"))
ML_I_RIS_sOLS <- summary(ML_I_RIS_OLS <- aov(fmOLS, data = df))
fmMLf <- LM$call$formula
fmMLr <- formula(paste("~1+",RIS_var,"|",RI_var))
ML_I_RIS_sMLE <- summary(ML_I_RIS_MLE <- lme(fmMLf, random = fmMLr, control=ctrl, data = df, method="ML"))
ML_I_RIS_sREML <- summary(ML_I_RIS_REML <- lme(fmMLf, random = fmMLr, control=ctrl, data = df, method="REML"))
l <- list(ML_I_RIS_OLS$Within,ML_I_RIS_sMLE,ML_I_RIS_sREML,ML_I_RIS_OLS,ML_I_RIS_MLE,ML_I_RIS_REML)
names(l) <- c("OLS I RIS (Within)","MLE I RIS","REML I RIS",
"OLS I RIS base","MLE I RIS base","REML I RIS base")
return(l)
}
A diferència del cas anterior no és tan fàcil aplicar la prova RLRT al haver dos efectes aleatoris. Es fa la prova ANOVA dels models corresponents ML amb 1 efecte aleatori (intercepció) i 2 (pendent i intercepció):
#NO ES PODEN FER ANOVAS DE RESULTATS DE FUNCIONS JA QUE NECESSITA LES FORMULES ORIGINALS, ES FA INDIVIDUALMENT AMB ELS MODELS AJUSTATS PER REML
M0_RIS<-Update_LM_RIS(M0,"Ident","Tiempo",LAKE1)
M0_RIS_REML<-lme(M0$call$formula, random = ~1+Tiempo|Ident, control=ctrl, data = LAKE1, method="REML")
(M0_RI_RIS_anREML <- anova(M0_RI_REML,M0_RIS_REML))
intM0_RIS<-Update_LM_RIS(intM0,"Ident","Tiempo",LAKE1)
intM0_RIS_REML<-lme(intM0$call$formula, random = ~1+Tiempo|Ident, control=ctrl, data = LAKE1, method="REML")
(intM0_RI_RIS_anREML <- anova(intM0_RI_REML,intM0_RIS_REML))
M0col_RIS<-Update_LM_RIS(M0col,"Ident","Tiempo",LAKE1)
M0col_RIS_REML<-lme(M0col$call$formula, random = ~1+Tiempo|Ident, control=ctrl, data = LAKE1, method="REML")
(M0col_RI_RIS_anREML <- anova(M0col_RI_REML,M0col_RIS_REML))
intM0col_RIS<-Update_LM_RIS(intM0col,"Ident","Tiempo",LAKE1)
intM0col_RIS_REML<-lme(intM0col$call$formula, random = ~1+Tiempo|Ident, control=ctrl, data = LAKE1, method="REML")
(intM0col_RI_RIS_anREML <- anova(intM0col_RI_REML,intM0col_RIS_REML))
M1_RIS<-Update_LM_RIS(M1,"Ident","Tiempo",LAKE1)
M1_RIS_REML<-lme(M1$call$formula, random = ~1+Tiempo|Ident, control=ctrl, data = LAKE1, method="REML")
(M1_RI_RIS_anREML <- anova(M1_RI_REML,M1_RIS_REML))
intM1_RIS<-Update_LM_RIS(intM1,"Ident","Tiempo",LAKE1)
intM1_RIS_REML<-lme(intM1$call$formula, random = ~1+Tiempo|Ident, control=ctrl, data = LAKE1, method="REML")
(intM1_RI_RIS_anREML <- anova(intM1_RI_REML,intM1_RIS_REML))
M2_RIS<-Update_LM_RIS(M2,"Ident","Tiempo",LAKE1)
M2_RIS_REML<-lme(M2$call$formula, random = ~1+Tiempo|Ident, control=ctrl, data = LAKE1, method="REML")
(M2_RI_RIS_anREML <- anova(M2_RI_REML,M2_RIS_REML))
intM2_RIS<-Update_LM_RIS(intM2,"Ident","Tiempo",LAKE1)
intM2_RIS_REML<-lme(intM2$call$formula, random = ~1+Tiempo|Ident, control=ctrl, data = LAKE1, method="REML")
(intM2_RI_RIS_anREML <- anova(intM2_RI_REML,intM2_RIS_REML))
#Generació dels objectes MLE i comparativa
M0_RIS_MLE<-lme(M0$call$formula, random = ~1+Tiempo|Ident, control=ctrl, data = LAKE1, method="ML")
(M0_RI_RIS_anMLE <- anova(M0_RI_MLE,M0_RIS_MLE))
intM0_RIS<-Update_LM_RIS(intM0,"Ident","Tiempo",LAKE1)
intM0_RIS_MLE<-lme(intM0$call$formula, random = ~1+Tiempo|Ident, control=ctrl, data = LAKE1, method="ML")
(intM0_RI_RIS_anMLE <- anova(intM0_RI_MLE,intM0_RIS_MLE))
M0col_RIS<-Update_LM_RIS(M0col,"Ident","Tiempo",LAKE1)
M0col_RIS_MLE<-lme(M0col$call$formula, random = ~1+Tiempo|Ident, control=ctrl, data = LAKE1, method="ML")
(M0col_RI_RIS_anMLE <- anova(M0col_RI_MLE,M0col_RIS_MLE))
intM0col_RIS<-Update_LM_RIS(intM0col,"Ident","Tiempo",LAKE1)
intM0col_RIS_MLE<-lme(intM0col$call$formula, random = ~1+Tiempo|Ident, control=ctrl, data = LAKE1, method="ML")
(intM0col_RI_RIS_anMLE <- anova(intM0col_RI_MLE,intM0col_RIS_MLE))
M1_RIS<-Update_LM_RIS(M1,"Ident","Tiempo",LAKE1)
M1_RIS_MLE<-lme(M1$call$formula, random = ~1+Tiempo|Ident, control=ctrl, data = LAKE1, method="ML")
(M1_RI_RIS_anMLE <- anova(M1_RI_MLE,M1_RIS_MLE))
intM1_RIS<-Update_LM_RIS(intM1,"Ident","Tiempo",LAKE1)
intM1_RIS_MLE<-lme(intM1$call$formula, random = ~1+Tiempo|Ident, control=ctrl, data = LAKE1, method="ML")
(intM1_RI_RIS_anMLE <- anova(intM1_RI_MLE,intM1_RIS_MLE))
M2_RIS<-Update_LM_RIS(M2,"Ident","Tiempo",LAKE1)
M2_RIS_MLE<-lme(M2$call$formula, random = ~1+Tiempo|Ident, control=ctrl, data = LAKE1, method="ML")
(M2_RI_RIS_anMLE <- anova(M2_RI_MLE,M2_RIS_MLE))
intM2_RIS<-Update_LM_RIS(intM2,"Ident","Tiempo",LAKE1)
intM2_RIS_MLE<-lme(intM2$call$formula, random = ~1+Tiempo|Ident, control=ctrl, data = LAKE1, method="ML")
(intM2_RI_RIS_anMLE <- anova(intM2_RI_MLE,intM2_RIS_MLE))
RI_RIS_an <- data.frame(Model = paste(rep(NomsModels,2),rep(c("RI_RIS_MLE","RI_RIS_REML"),each=length(NomsModels))))
RI_RIS_an$Model.base <- rep(NomsModels,2)
RI_RIS_an$ML.method <- rep(c("RI_RIS_MLE","RI_RIS_REML"),each=length(NomsModels))
RI_RIS_an$p.valor <- c(M0_RI_RIS_anMLE$`p-value`[2], intM0_RI_RIS_anMLE$`p-value`[2], M0col_RI_RIS_anMLE$`p-value`[2],
intM0col_RI_RIS_anMLE$`p-value`[2], M1_RI_RIS_anMLE$`p-value`[2], intM1_RI_RIS_anMLE$`p-value`[2],
M2_RI_RIS_anMLE$`p-value`[2], intM2_RI_RIS_anMLE$`p-value`[2],
M0_RI_RIS_anREML$`p-value`[2], intM0_RI_RIS_anREML$`p-value`[2], M0col_RI_RIS_anREML$`p-value`[2],
intM0col_RI_RIS_anREML$`p-value`[2], M1_RI_RIS_anREML$`p-value`[2], intM1_RI_RIS_anREML$`p-value`[2],
M2_RI_RIS_anREML$`p-value`[2], intM2_RI_RIS_anREML$`p-value`[2])
RI_RIS_an_plot <- ggplot(data=RI_RIS_an, aes(x=Model.base, y=p.valor, fill=ML.method)) +
geom_bar(stat="identity", position=position_dodge())+
scale_fill_brewer(palette="Blues") + theme_minimal() +
labs(y = "RI/RIS anova p-valor") +
theme(axis.text.x = element_text(angle=80,size = 10,hjust=0))
RI_RIS_an_plot
En la majoria de models segueix essent la significació baixa (p-valor alt), però es veu curiosament que hi ha alguns models on es torna o queda proper a la significància, concretament si es pren els que tenen p-valors per sota de 0.25 serien M.Base Int. RI_RIS_MLE, M.LogTrans.Int. RI_RIS_MLE, M.Base Int. RI_RIS_REML, M.LogTrans.Int. RI_RIS_REML amb p-valors de 0.0541, 0.0236, 0.0702, 0.0215 respectivament.
Es tornen a generar els diagnòstics per aquests models amb efecte RIS:
M_RIS_listA0 <- list(M0,intM0,M0col,intM0col,
M0_RIS[1:3],intM0_RIS[1:3],M0col_RIS[1:3],intM0col_RIS[1:3])
M_RIS_listA <- ConvList(M_RIS_listA0)
## Warning in if (class(List[[l1]]) == "list") {: the condition has length > 1
## and only the first element will be used
## Warning in if (class(List[[l1]]) == "list") {: the condition has length > 1
## and only the first element will be used
## Warning in if (class(List[[l1]]) == "list") {: the condition has length > 1
## and only the first element will be used
## Warning in if (class(List[[l1]]) == "list") {: the condition has length > 1
## and only the first element will be used
M_RIS_listB0 <- list(M1,intM1,M2,intM2,
M1_RIS[1:3],intM1_RIS[1:3],M2_RIS[1:3],intM2_RIS[1:3])
M_RIS_listB <- ConvList(M_RIS_listB0)
## Warning in if (class(List[[l1]]) == "list") {: the condition has length > 1
## and only the first element will be used
## Warning in if (class(List[[l1]]) == "list") {: the condition has length > 1
## and only the first element will be used
## Warning in if (class(List[[l1]]) == "list") {: the condition has length > 1
## and only the first element will be used
## Warning in if (class(List[[l1]]) == "list") {: the condition has length > 1
## and only the first element will be used
names(M_RIS_listA)[1:4] <- NomsModels[1:4]
names(M_RIS_listA)[(4+1):length(M_RIS_listA)]<-paste(rep(NomsModels[1:4],each=3),rep(paste("RIS",c("OLS","MLE","REML")),4))
names(M_RIS_listB)[1:4] <- NomsModels[5:8]
names(M_RIS_listB)[(4+1):length(M_RIS_listB)]<-paste(rep(NomsModels[5:8],each=3),rep(paste("RIS",c("OLS","MLE","REML")),4))
(DiagA_RIS <- ML_Resum(M_RIS_listA,"Lake1"))
## [[1]]
## Sigma AIC BIC
## Lake1--M.Base 81386.03 1023.8574 1032.3018
## Lake1--M.Base Int. 47433.69 983.9564 999.1563
## Lake1--M.Corr.Col 84579.43 1026.0323 1032.7878
## Lake1--M.Corr.Col. Int. 84579.43 1026.0323 1032.7878
## Lake1--M.Base RIS OLS 71784.03 608.7288 612.2630
## Lake1--M.Base RIS MLE 73628.98 1029.9517 1043.4628
## Lake1--M.Base RIS REML 75424.20 957.7893 970.4574
## Lake1--M.Base Int. RIS OLS 38378.48 581.8568 590.1032
## Lake1--M.Base Int. RIS MLE 32069.02 984.1226 1004.3891
## Lake1--M.Base Int. RIS REML 35948.69 872.5253 890.1141
## Lake1--M.Corr.Col RIS OLS 90625.78 618.9833 621.3394
## Lake1--M.Corr.Col RIS MLE 81343.91 1032.0368 1043.8590
## Lake1--M.Corr.Col RIS REML 81534.59 977.4262 988.7026
## Lake1--M.Corr.Col. Int. RIS OLS 90625.78 618.9833 621.3394
## Lake1--M.Corr.Col. Int. RIS MLE 81343.91 1032.0368 1043.8590
## Lake1--M.Corr.Col. Int. RIS REML 81534.59 977.4262 988.7026
##
## [[2]]
BestA2 <- unique(c(order(DiagA_RIS[[1]]$AIC,decreasing=F)[1:4],
order(DiagA_RIS[[1]]$Sigma,decreasing=F)[1:4]))
ML_Resum(M_RIS_listA[sort(BestA2)],"Lake1") #Comparativa més selectiva
## [[1]]
## Sigma AIC BIC
## Lake1--M.Base Int. 47433.69 983.9564 999.1563
## Lake1--M.Base RIS OLS 71784.03 608.7288 612.2630
## Lake1--M.Base Int. RIS OLS 38378.48 581.8568 590.1032
## Lake1--M.Base Int. RIS MLE 32069.02 984.1226 1004.3891
## Lake1--M.Base Int. RIS REML 35948.69 872.5253 890.1141
## Lake1--M.Corr.Col RIS OLS 90625.78 618.9833 621.3394
## Lake1--M.Corr.Col. Int. RIS OLS 90625.78 618.9833 621.3394
##
## [[2]]
(DiagB_RIS<-ML_Resum(M_RIS_listB,"Lake1"))
## [[1]]
## Sigma AIC BIC
## Lake1--M.LogTrans. 2.1784 179.7516 184.8182
## Lake1--M.LogTrans.Int. 1.8898 170.2175 178.6619
## Lake1--M.Box-Cox 296.4675 572.8184 577.8850
## Lake1--M.Box-Cox Int. 277.9218 569.4879 577.9322
## Lake1--M.LogTrans. RIS OLS 2.5788 115.5800 116.7580
## Lake1--M.LogTrans. RIS MLE 2.1233 185.7520 195.8853
## Lake1--M.LogTrans. RIS REML 2.1784 192.0543 201.8799
## Lake1--M.LogTrans.Int. RIS OLS 1.5054 91.6563 95.1905
## Lake1--M.LogTrans.Int. RIS MLE 1.3307 168.7252 182.2362
## Lake1--M.LogTrans.Int. RIS REML 1.3916 189.3124 201.9806
## Lake1--M.Box-Cox RIS OLS 344.1067 350.4747 351.6528
## Lake1--M.Box-Cox RIS MLE 288.9329 578.8238 588.9571
## Lake1--M.Box-Cox RIS REML 296.4283 565.4728 575.2983
## Lake1--M.Box-Cox Int. RIS OLS 265.9007 340.0107 343.5448
## Lake1--M.Box-Cox Int. RIS MLE 249.5092 575.5950 589.1060
## Lake1--M.Box-Cox Int. RIS REML 255.5007 556.0538 568.7220
##
## [[2]]
BestB2 <- unique(c(order(DiagB_RIS[[1]]$AIC,decreasing=F)[1:4],
order(DiagB_RIS[[1]]$Sigma,decreasing=F)[1:4]))
ML_Resum(M_RIS_listB[sort(BestB2)],"Lake1") #Comparativa més selectiva
## [[1]]
## Sigma AIC BIC
## Lake1--M.LogTrans.Int. 1.8898 170.2175 178.6619
## Lake1--M.LogTrans. RIS OLS 2.5788 115.5800 116.7580
## Lake1--M.LogTrans.Int. RIS OLS 1.5054 91.6563 95.1905
## Lake1--M.LogTrans.Int. RIS MLE 1.3307 168.7252 182.2362
## Lake1--M.LogTrans.Int. RIS REML 1.3916 189.3124 201.9806
##
## [[2]]
#Deixem només de base els d'interacció i fora el de correcció de colinealitat també:
(DiagAi_RIS <- ML_Resum(M_RIS_listA[-c(1,3:7,11:16)],"Lake1"))
## [[1]]
## Sigma AIC BIC
## Lake1--M.Base Int. 47433.69 983.9564 999.1563
## Lake1--M.Base Int. RIS OLS 38378.48 581.8568 590.1032
## Lake1--M.Base Int. RIS MLE 32069.02 984.1226 1004.3891
## Lake1--M.Base Int. RIS REML 35948.69 872.5253 890.1141
##
## [[2]]
BestAi <- unique(c(order(DiagAi_RIS[[1]]$AIC,decreasing=F)[1:4],
order(DiagAi_RIS[[1]]$Sigma,decreasing=F)[1:4]))
ML_Resum((M_RIS_listA[-c(1,3:7,11:16)])[sort(BestAi)],"Lake1") #Comparativa més selectiva
## [[1]]
## Sigma AIC BIC
## Lake1--M.Base Int. 47433.69 983.9564 999.1563
## Lake1--M.Base Int. RIS OLS 38378.48 581.8568 590.1032
## Lake1--M.Base Int. RIS MLE 32069.02 984.1226 1004.3891
## Lake1--M.Base Int. RIS REML 35948.69 872.5253 890.1141
##
## [[2]]
(DiagBi_RIS<-ML_Resum(M_RIS_listB[-c(1,3,5:7,11:13)],"Lake1"))
## [[1]]
## Sigma AIC BIC
## Lake1--M.LogTrans.Int. 1.8898 170.2175 178.6619
## Lake1--M.Box-Cox Int. 277.9218 569.4879 577.9322
## Lake1--M.LogTrans.Int. RIS OLS 1.5054 91.6563 95.1905
## Lake1--M.LogTrans.Int. RIS MLE 1.3307 168.7252 182.2362
## Lake1--M.LogTrans.Int. RIS REML 1.3916 189.3124 201.9806
## Lake1--M.Box-Cox Int. RIS OLS 265.9007 340.0107 343.5448
## Lake1--M.Box-Cox Int. RIS MLE 249.5092 575.5950 589.1060
## Lake1--M.Box-Cox Int. RIS REML 255.5007 556.0538 568.7220
##
## [[2]]
BestBi <- unique(c(order(DiagBi_RIS[[1]]$AIC,decreasing=F)[1:4],
order(DiagBi_RIS[[1]]$Sigma,decreasing=F)[1:4]))
ML_Resum((M_RIS_listB[-c(1,3,5:7,11:13)])[sort(BestBi)],"Lake1") #Comparativa més selectiva
## [[1]]
## Sigma AIC BIC
## Lake1--M.LogTrans.Int. 1.8898 170.2175 178.6619
## Lake1--M.LogTrans.Int. RIS OLS 1.5054 91.6563 95.1905
## Lake1--M.LogTrans.Int. RIS MLE 1.3307 168.7252 182.2362
## Lake1--M.LogTrans.Int. RIS REML 1.3916 189.3124 201.9806
##
## [[2]]
Pels models seleccionats amb una variància residual més baixa així com els indicadors AIC/BIC més baixos curiosament es pot veure que pels 2 models que havien donat una significació alta el fet d’afegir l’efecte aleatori RIS, aparentment es nota d’una manera més efectiva la millora en aquests paràmetres al comparar-los amb la referència base. No és una millora espectacular, però podria ser suficient per fer prediccions més precises.
Es calcula com abans el coeficient de determinació \(R^2\). En el cas dels models amb efectes aleatoris es calcula a partir del coeficient de determinació condicionat que suma el que aporta els efectes fixos i els efectes aleatoris (en el cas d’aquest models només s’ha trobat fàcilment com calcular-ho en R amb els models calculats per ML i dins d’aquí en els casos amb termes d’interacció es fa una aproximació mitjançant la correlació lineal entre la resposta i els valors ajustats):
M_RIS_listA_R2 <- list(M0,intM0,M0col,intM0col,M0_RIS_MLE,M0_RIS_REML,
intM0_RIS_MLE,intM0_RIS_REML,M0col_RIS_MLE,M0col_RIS_REML,
intM0col_RIS_MLE,intM0col_RIS_REML)
names(M_RIS_listA_R2)[1:4] <- NomsModels[1:4]
names(M_RIS_listA_R2)[(4+1):length(M_RIS_listA_R2)]<-paste(rep(NomsModels[1:4],each=2),rep(paste("RIS",c("MLE","REML")),4))
M_RIS_dfA_R2 <- data.frame(Model = names(M_RIS_listA_R2))
M_RIS_dfA_R2$R2 <- sapply(1:length(M_RIS_listA_R2),function(M){
if (class(M_RIS_listA_R2[[M]])[1]=="lmm"){
summary(M_RIS_listA_R2[[M]])$r.squared
} else if (class(M_RIS_listA_R2[[M]])=="lme"){
if(all(paste(attr(M_RIS_listA_R2[[M]]$terms,"predvars")[-c(1:2)]) ==
paste(attr(M_RIS_listA_R2[[M]]$terms,"term.labels")))){
rsquared(M_RIS_listA_R2[[M]])$Conditional
} else {
summary(lm(formula(paste(paste("M_RIS_listA_R2[[",M,"]]$data$",
attr(M_RIS_listA_R2[[M]]$terms,"variables")[[2]],sep=""),
"~","M_RIS_listA_R2[[",M,"]]$fitted"),sep="")))$r.squared
}
}
})
## Warning in paste(attr(M_RIS_listA_R2[[M]]$terms, "predvars")[-c(1:2)]) == :
## longer object length is not a multiple of shorter object length
## Warning in paste(attr(M_RIS_listA_R2[[M]]$terms, "predvars")[-c(1:2)]) == :
## longer object length is not a multiple of shorter object length
R2_A2 <- ggplot(data=M_RIS_dfA_R2, aes(x=Model, y=R2, fill=Model)) +
geom_bar(stat="identity", position=position_dodge())+
scale_fill_brewer(palette="Set3") + theme_minimal() +
labs(y = "R2 Coeff.") +
theme(axis.text.x = element_text(angle=80,size = 10,hjust=0))
M_RIS_listB_R2 <- list(M1,intM1,M2,intM2,M1_RIS_MLE,
M1_RIS_REML,intM1_RIS_MLE,intM1_RIS_REML,
M2_RIS_MLE,M2_RIS_REML,intM2_RIS_MLE,intM2_RIS_REML)
names(M_RIS_listB_R2)[1:4] <- NomsModels[5:8]
names(M_RIS_listB_R2)[(4+1):length(M_RIS_listB_R2)]<-paste(rep(NomsModels[5:8],each=2),rep(paste("RIS",c("MLE","REML")),4))
M_RIS_dfB_R2 <- data.frame(Model = names(M_RIS_listB_R2))
M_RIS_dfB_R2 <- data.frame(Model = names(M_RIS_listB_R2))
M_RIS_dfB_R2$R2 <- sapply(1:length(M_RIS_listB_R2),function(M){
if (class(M_RIS_listB_R2[[M]])[1]=="lmm"){
summary(M_RIS_listB_R2[[M]])$r.squared
} else if (class(M_RIS_listB_R2[[M]])=="lme"){
if(all(paste(attr(M_RIS_listB_R2[[M]]$terms,"predvars")[-c(1:2)]) ==
paste(attr(M_RIS_listB_R2[[M]]$terms,"term.labels")))){
rsquared(M_RIS_listB_R2[[M]])$Conditional
} else {
summary(lm(formula(paste(paste("M_RIS_listB_R2[[",M,"]]$data$",
attr(M_RIS_listB_R2[[M]]$terms,"variables")[[2]],sep=""),
"~","M_RIS_listB_R2[[",M,"]]$fitted"),sep="")))$r.squared
}
}
})
## Warning in paste(attr(M_RIS_listB_R2[[M]]$terms, "predvars")[-c(1:2)]) == :
## longer object length is not a multiple of shorter object length
## Warning in paste(attr(M_RIS_listB_R2[[M]]$terms, "predvars")[-c(1:2)]) == :
## longer object length is not a multiple of shorter object length
## Warning in paste(attr(M_RIS_listB_R2[[M]]$terms, "predvars")[-c(1:2)]) == :
## longer object length is not a multiple of shorter object length
## Warning in paste(attr(M_RIS_listB_R2[[M]]$terms, "predvars")[-c(1:2)]) == :
## longer object length is not a multiple of shorter object length
R2_B2 <- ggplot(data=M_RIS_dfB_R2, aes(x=Model, y=R2, fill=Model)) +
geom_bar(stat="identity", position=position_dodge())+
scale_fill_brewer(palette="Set3") + theme_minimal() +
labs(y = "R2 Coeff.") +
theme(axis.text.x = element_text(angle=80,size = 10,hjust=0))
R2_A2
R2_B2
BestA2_R2 <- order(M_RIS_dfA_R2$R2,decreasing=T)[1:4]
BestB2_R2 <- order(M_RIS_dfB_R2$R2,decreasing=T)[1:4]
M_RIS_df_sel_R2 <- rbind(M_RIS_dfA_R2[BestA2_R2,],M_RIS_dfB_R2[BestB2_R2,])
R2_sel2 <- ggplot(data=M_RIS_df_sel_R2, aes(x=Model, y=R2, fill=Model)) +
geom_bar(stat="identity", position=position_dodge())+
scale_fill_brewer(palette="Set3") + theme_minimal() +
labs(y = "R2 Coeff.") +
theme(axis.text.x = element_text(angle=80,size = 10,hjust=0))
R2_sel2
#Deixem només els models amb interacció
M_RIS_dfA_R2i <- M_RIS_dfA_R2[-c(1,3:6,9:12),]
R2_A2i <- ggplot(data=M_RIS_dfA_R2i, aes(x=Model, y=R2, fill=Model)) +
geom_bar(stat="identity", position=position_dodge())+
scale_fill_brewer(palette="Set3") + theme_minimal() +
labs(y = "R2 Coeff.") +
theme(axis.text.x = element_text(angle=80,size = 10,hjust=0))
M_RIS_dfB_R2i <- M_RIS_dfB_R2[-c(1,3,5:6,9:10),]
R2_B2i <- ggplot(data=M_RIS_dfB_R2i, aes(x=Model, y=R2, fill=Model)) +
geom_bar(stat="identity", position=position_dodge())+
scale_fill_brewer(palette="Set3") + theme_minimal() +
labs(y = "R2 Coeff.") +
theme(axis.text.x = element_text(angle=80,size = 10,hjust=0))
R2_A2i
R2_B2i
M_RIS_df_sel_R2i <- rbind(M_RIS_dfA_R2i,M_RIS_dfB_R2i)
R2_sel2i <- ggplot(data=M_RIS_df_sel_R2i, aes(x=Model, y=R2, fill=Model)) +
geom_bar(stat="identity", position=position_dodge())+
scale_fill_brewer(palette="Set3") + theme_minimal() +
labs(y = "R2 Coeff.") +
theme(axis.text.x = element_text(angle=80,size = 10,hjust=0))
R2_sel2i
Sí que s’observa aquesta vegada una millora més aparent en el coeficient \(R^2\), concretament si es pren els 4 models amb un coeficient més alt serien M.Base Int. RIS REML, M.Base Int. RIS MLE, M.LogTrans.Int. RIS REML, M.LogTrans.Int. RIS MLE amb coeficients de 0.9211, 0.9208, 0.8564, 0.8556 respectivament. Curiosament segueixen essent els que denotaven millora amb l’efecte aleatori.
Per els següents apartats de codi es fa una selecció dels models en base a les proves que s’han anat fent en aquest últim pas pel que es tindrà en compte els models següents:
#Els models que s¡utilitzaran per continuar amb l'anàlisi son
Models.sel <- list(intM0_RIS_MLE,intM0_RIS_REML,
intM1_RIS_MLE,intM1_RIS_REML,
intM2_RIS_MLE,intM2_RIS_REML)
names(Models.sel) <- c("M.Int.RIS.MLE","M.Int.RIS.REML",
"M.Logt.Int.RIS.MLE","M.Logt.Int.RIS.REML",
"M.BCt.Int.RIS.MLE","M.BCt.Int.RIS.REML")
En aquest apartat on es denota més efecte per l’efecte aleatori RIS també és interessant veure si ha canviat gaire el diagnòstic clàssic del model lineal. Es fa en els models seleccionats únicament:
Lake1IntRIS_MLE.SM0.sdres <- residuals(intM0_RIS_MLE,type = "pearson")
Lake1IntRIS_REML.SM0.sdres <- residuals(intM0_RIS_REML,type = "pearson")
par(mfrow=c(1,3))
qqnorm(Lake1Int.SM0.sdres,
ylab="Standardized Res.",
xlab="Normal Scores",
main="M. Original Int. - Standardized Res.")
qqline(Lake1Int.SM0.sdres)
qqnorm(Lake1IntRIS_MLE.SM0.sdres,
ylab="Standardized Res.",
xlab="Normal Scores",
main="M. Original Int RIS MLE - Standardized Res.")
qqline(Lake1Int.SM0.sdres)
qqnorm(Lake1IntRIS_REML.SM0.sdres,
ylab="Standardized Res.",
xlab="Normal Scores",
main="M. Original Int RIS REML - Standardized Res.")
qqline(Lake1Int.SM0.sdres)
S0Int
##
## Shapiro-Wilk normality test
##
## data: intM0$residuals
## W = 0.92773, p-value = 0.01359
(S0IntRIS_MLE <- shapiro.test(intM0_RIS_MLE$residuals))
##
## Shapiro-Wilk normality test
##
## data: intM0_RIS_MLE$residuals
## W = 0.93961, p-value = 0.0009178
(S0IntRIS_REML <- shapiro.test(intM0_RIS_REML$residuals))
##
## Shapiro-Wilk normality test
##
## data: intM0_RIS_REML$residuals
## W = 0.9398, p-value = 0.0009395
Lake1IntRIS_MLE.SM1.sdres <- residuals(intM1_RIS_MLE,type = "pearson")
Lake1IntRIS_REML.SM1.sdres <- residuals(intM1_RIS_REML,type = "pearson")
par(mfrow=c(1,3))
qqnorm(Lake1Int.SM1.sdres,
ylab="Standardized Res.",
xlab="Normal Scores",
main="M. Logt Int. - Standardized Res.")
qqline(Lake1Int.SM1.sdres)
qqnorm(Lake1IntRIS_MLE.SM1.sdres,
ylab="Standardized Res.",
xlab="Normal Scores",
main="M. Logt Int RIS MLE - Standardized Res.")
qqline(Lake1Int.SM1.sdres)
qqnorm(Lake1IntRIS_REML.SM1.sdres,
ylab="Standardized Res.",
xlab="Normal Scores",
main="M. Logt Int RIS REML - Standardized Res.")
qqline(Lake1Int.SM1.sdres)
S1Int
##
## Shapiro-Wilk normality test
##
## data: intM1$residuals
## W = 0.9557, p-value = 0.1192
(S1IntRIS_MLE <- shapiro.test(intM1_RIS_MLE$residuals))
##
## Shapiro-Wilk normality test
##
## data: intM1_RIS_MLE$residuals
## W = 0.95018, p-value = 0.003532
(S1IntRIS_REML <- shapiro.test(intM1_RIS_REML$residuals))
##
## Shapiro-Wilk normality test
##
## data: intM1_RIS_REML$residuals
## W = 0.95009, p-value = 0.003488
Lake1IntRIS_MLE.SM2.sdres <- residuals(intM2_RIS_MLE,type = "pearson")
Lake1IntRIS_REML.SM2.sdres <- residuals(intM2_RIS_REML,type = "pearson")
par(mfrow=c(1,3))
qqnorm(Lake1Int.SM2.sdres,
ylab="Standardized Res.",
xlab="Normal Scores",
main="M. BCt Int. - Standardized Res.")
qqline(Lake1Int.SM2.sdres)
qqnorm(Lake1IntRIS_MLE.SM2.sdres,
ylab="Standardized Res.",
xlab="Normal Scores",
main="M. BCt Int RIS MLE - Standardized Res.")
qqline(Lake1Int.SM2.sdres)
qqnorm(Lake1IntRIS_REML.SM2.sdres,
ylab="Standardized Res.",
xlab="Normal Scores",
main="M. BCt Int RIS REML - Standardized Res.")
qqline(Lake1Int.SM2.sdres)
S2Int
##
## Shapiro-Wilk normality test
##
## data: intM2$residuals
## W = 0.96332, p-value = 0.2172
(S2IntRIS_MLE <- shapiro.test(intM2_RIS_MLE$residuals))
##
## Shapiro-Wilk normality test
##
## data: intM2_RIS_MLE$residuals
## W = 0.96907, p-value = 0.04949
(S2IntRIS_REML <- shapiro.test(intM2_RIS_REML$residuals))
##
## Shapiro-Wilk normality test
##
## data: intM2_RIS_REML$residuals
## W = 0.97436, p-value = 0.1078
El que es veu en el test Shapiro no seria bo ja que sembla que els models s’allunyen de la normalitat al afegir els efectes aleatoris RIS, tot i així visualment es pot veure que hi ha poques diferències en el QQplot, pel que potser seria bo no prendre’s aquest test de Shapiro com definitiu (en cas que es volgués aprofundir més es podria intentar fer més proves de normalitat).
pHIntRIS_MLE01 <- ggplot(LAKE1, aes(x=fitted(intM0_RIS_MLE), y=residuals(intM0_RIS_MLE),
shape=Grupo, color=Grupo)) +
geom_jitter(position=position_jitter(1), cex=1.5) +
geom_hline(yintercept=0, linetype="dashed", color = "black") +
labs(title="Fitted vs Residuals (IntRIS_MLE.Homocedasticity test)",
x="Fitted", y = "Residuals")
pHIntRIS_MLE02 <-ggplot(LAKE1, aes(x=fitted(intM0_RIS_MLE), y=sqrt(abs(residuals(intM0_RIS_MLE))),
shape=Grupo, color=Grupo)) +
geom_jitter(position=position_jitter(1), cex=1.5) +
geom_hline(yintercept=0, linetype="dashed", color = "black") +
labs(title="Fitted vs Residuals (IntRIS_MLE.Homocedasticity test)",
x="Fitted", y = expression(sqrt("|residuals|")))
pHIntRIS_REML01 <- ggplot(LAKE1, aes(x=fitted(intM0_RIS_REML), y=residuals(intM0_RIS_REML),
shape=Grupo, color=Grupo)) +
geom_jitter(position=position_jitter(1), cex=1.5) +
geom_hline(yintercept=0, linetype="dashed", color = "black") +
labs(title="Fitted vs Residuals (IntRIS_REML.Homocedasticity test)",
x="Fitted", y = "Residuals")
pHIntRIS_REML02 <-ggplot(LAKE1, aes(x=fitted(intM0_RIS_REML), y=sqrt(abs(residuals(intM0_RIS_REML))),
shape=Grupo, color=Grupo)) +
geom_jitter(position=position_jitter(1), cex=1.5) +
geom_hline(yintercept=0, linetype="dashed", color = "black") +
labs(title="Fitted vs Residuals (IntRIS_REML.Homocedasticity test)",
x="Fitted", y = expression(sqrt("|residuals|")))
grid.arrange(pHInt01, pHIntRIS_MLE01,pHIntRIS_REML01,nrow=3)
grid.arrange(pHInt02, pHIntRIS_MLE02,pHIntRIS_REML02,nrow=3)
plot(intM0,which=1)
plot(intM0_RIS_MLE,which=1)
plot(intM0_RIS_REML,which=1)
pHIntRIS_MLE11 <- ggplot(LAKE1, aes(x=fitted(intM1_RIS_MLE), y=residuals(intM1_RIS_MLE),
shape=Grupo, color=Grupo)) +
geom_jitter(position=position_jitter(1), cex=1.5) +
geom_hline(yintercept=0, linetype="dashed", color = "black") +
labs(title="Fitted vs Residuals (Logt.IntRIS_MLE.Homocedasticity test)",
x="Fitted", y = "Residuals")
pHIntRIS_MLE12 <-ggplot(LAKE1, aes(x=fitted(intM1_RIS_MLE), y=sqrt(abs(residuals(intM1_RIS_MLE))),
shape=Grupo, color=Grupo)) +
geom_jitter(position=position_jitter(1), cex=1.5) +
geom_hline(yintercept=0, linetype="dashed", color = "black") +
labs(title="Fitted vs Residuals (Logt.IntRIS_MLE.Homocedasticity test)",
x="Fitted", y = expression(sqrt("|residuals|")))
pHIntRIS_REML11 <- ggplot(LAKE1, aes(x=fitted(intM1_RIS_REML), y=residuals(intM1_RIS_REML),
shape=Grupo, color=Grupo)) +
geom_jitter(position=position_jitter(1), cex=1.5) +
geom_hline(yintercept=0, linetype="dashed", color = "black") +
labs(title="Fitted vs Residuals (Logt.IntRIS_REML.Homocedasticity test)",
x="Fitted", y = "Residuals")
pHIntRIS_REML12 <-ggplot(LAKE1, aes(x=fitted(intM1_RIS_REML), y=sqrt(abs(residuals(intM1_RIS_REML))),
shape=Grupo, color=Grupo)) +
geom_jitter(position=position_jitter(1), cex=1.5) +
geom_hline(yintercept=0, linetype="dashed", color = "black") +
labs(title="Fitted vs Residuals (Logt.IntRIS_REML.Homocedasticity test)",
x="Fitted", y = expression(sqrt("|residuals|")))
grid.arrange(pHInt11, pHIntRIS_MLE11,pHIntRIS_REML11,nrow=3)
grid.arrange(pHInt12, pHIntRIS_MLE12,pHIntRIS_REML12,nrow=3)
plot(intM1,which=1)
plot(intM1_RIS_MLE,which=1)
plot(intM1_RIS_REML,which=1)
pHIntRIS_MLE21 <- ggplot(LAKE1, aes(x=fitted(intM2_RIS_MLE), y=residuals(intM2_RIS_MLE),
shape=Grupo, color=Grupo)) +
geom_jitter(position=position_jitter(1), cex=1.5) +
geom_hline(yintercept=0, linetype="dashed", color = "black") +
labs(title="Fitted vs Residuals (BCt.IntRIS_MLE.Homocedasticity test)",
x="Fitted", y = "Residuals")
pHIntRIS_MLE22 <-ggplot(LAKE1, aes(x=fitted(intM2_RIS_MLE), y=sqrt(abs(residuals(intM2_RIS_MLE))),
shape=Grupo, color=Grupo)) +
geom_jitter(position=position_jitter(1), cex=1.5) +
geom_hline(yintercept=0, linetype="dashed", color = "black") +
labs(title="Fitted vs Residuals (BCt.IntRIS_MLE.Homocedasticity test)",
x="Fitted", y = expression(sqrt("|residuals|")))
pHIntRIS_REML21 <- ggplot(LAKE1, aes(x=fitted(intM2_RIS_REML), y=residuals(intM2_RIS_REML),
shape=Grupo, color=Grupo)) +
geom_jitter(position=position_jitter(1), cex=1.5) +
geom_hline(yintercept=0, linetype="dashed", color = "black") +
labs(title="Fitted vs Residuals (BCt.IntRIS_REML.Homocedasticity test)",
x="Fitted", y = "Residuals")
pHIntRIS_REML22 <-ggplot(LAKE1, aes(x=fitted(intM2_RIS_REML), y=sqrt(abs(residuals(intM2_RIS_REML))),
shape=Grupo, color=Grupo)) +
geom_jitter(position=position_jitter(1), cex=1.5) +
geom_hline(yintercept=0, linetype="dashed", color = "black") +
labs(title="Fitted vs Residuals (BCt.IntRIS_REML.Homocedasticity test)",
x="Fitted", y = expression(sqrt("|residuals|")))
grid.arrange(pHInt21, pHIntRIS_MLE21,pHIntRIS_REML21,nrow=3)
grid.arrange(pHInt22, pHIntRIS_MLE22,pHIntRIS_REML22,nrow=3)
plot(intM2,which=1)
plot(intM2_RIS_MLE,which=1)
plot(intM2_RIS_REML,which=1)
Tampoc es mostra una diferència aparentment significativa en l’homocedasticitat que ja tenien els models de referència. Sembla que segueix fent falta algun tipus de modulació de la variància excepte potser en algun dels models transformats.
Es recupera una de les funcions generades en altres simulacions per obtenir els models amb funcions de variància amb algunes modificaciones per adaptar-la a l’objectiu d’aquest codi:
Comp_W_LM <- function(fLME,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 <- update(fLME, control=ctrl, method="ML")
baseREML <- update(fLME, control=ctrl, method="REML")
l <- list()
fMLk <- "IS"
fREk <- "RIS"
if(WVI==1 | WVI==3){
ML_I_RI_WVI_sMLE <- summary(ML_I_RI_WVI_MLE <- update(baseML, weights = varIdent(form = ~Tiempo)))
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 = ~Tiempo)))
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 = ~Tiempo)))
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 = ~Tiempo)))
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)
}
Com en les simulacions es va trobar que no sempre convergien en una resposta es fa algunes proves per veure en quins models podem ajustar aquestes funcions:
#Models.sel <- list(intM0_RIS_MLE,intM0_RIS_REML, intM1_RIS_MLE,intM1_RIS_REML, intM2_RIS_MLE,intM2_RIS_REML)
intM0_RIS_W <- Comp_W_LM(intM0_RIS_MLE,3,3,0)
## 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
## 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
## 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
intM1_RIS_W <- Comp_W_LM(intM1_RIS_MLE,3,3,3)
## 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
## 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
## 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
intM2_RIS_W <- Comp_W_LM(intM2_RIS_MLE,3,3,2)
## 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
## Warning in logLik.reStruct(object, conLin): Singular precision matrix in
## level -1, block 1
intM0_RIS_Ws0 <- list(Models.sel[c(1:2)],intM0_RIS_W[c(1,3,5,7)])
intM0_RIS_Ws <- ConvList(intM0_RIS_Ws0)
(intM0_QI<-ML_Resum(intM0_RIS_Ws,"CargaV.")) #Es guarda la taula per us posterior
## [[1]]
## Sigma AIC BIC
## CargaV.--M.Int.RIS.MLE 32069.02 984.1226 1004.3891
## CargaV.--M.Int.RIS.REML 35948.69 872.5253 890.1141
## CargaV.--IS MLE RIS VarIdent 32069.02 984.1226 1004.3891
## CargaV.--IS REML RIS VarIdent 35948.69 872.5253 890.1141
## CargaV.--IS MLE RIS VarExp 234497.34 567.8461 589.8015
## CargaV.--IS REML RIS VarExp 198543.08 594.2267 613.2813
##
## [[2]]
ML_Resum(intM0_RIS_Ws[1:4],"CargaV.")[[2]]
intM1_RIS_Ws0 <- list(Models.sel[c(3:4)],intM1_RIS_W[c(1,3,5,7,9,11)])
intM1_RIS_Ws <- ConvList(intM1_RIS_Ws0)
(intM1_QI <- ML_Resum(intM1_RIS_Ws,"LogCargaV."))
## [[1]]
## Sigma AIC BIC
## LogCargaV.--M.Logt.Int.RIS.MLE 1.3307 168.7252 182.2362
## LogCargaV.--M.Logt.Int.RIS.REML 1.3916 189.3124 201.9806
## LogCargaV.--IS MLE RIS VarIdent 1.3307 168.7252 182.2362
## LogCargaV.--IS REML RIS VarIdent 1.3916 189.3124 201.9806
## LogCargaV.--IS MLE RIS VarExp 5.8062 93.8243 109.0242
## LogCargaV.--IS REML RIS VarExp 5.8149 127.6175 141.8692
## LogCargaV.--IS MLE RIS VarPower 0.0000 116.8606 132.0605
## LogCargaV.--IS REML RIS VarPower 0.0000 148.6902 162.9418
##
## [[2]]
ML_Resum(intM1_RIS_Ws[-c(5:6)],"LogCargaV.")[[2]]
intM2_RIS_Ws0 <- list(Models.sel[c(1:2)],intM2_RIS_W[c(1,3,5,7,11)])
intM2_RIS_Ws <- ConvList(intM2_RIS_Ws0)
(intM2_QI <- ML_Resum(intM2_RIS_Ws,"BCtCargaV."))
## [[1]]
## Sigma AIC BIC
## BCtCargaV.--M.Int.RIS.MLE 32069.0244 984.1226 1004.3891
## BCtCargaV.--M.Int.RIS.REML 35948.6941 872.5253 890.1141
## BCtCargaV.--IS MLE RIS VarIdent 249.5092 575.5950 589.1060
## BCtCargaV.--IS REML RIS VarIdent 255.5007 556.0538 568.7220
## BCtCargaV.--IS MLE RIS VarExp 657.1780 557.3659 572.5659
## BCtCargaV.--IS REML RIS VarExp 664.2374 541.3619 555.6136
## BCtCargaV.--IS REML RIS VarPower 0.0000 546.2706 560.5222
##
## [[2]]
ML_Resum(intM2_RIS_Ws[-c(1:2)],"BCtCargaV.")[[2]]
En aquests primers diagnòstics semblaria que la funció varPower és la que té més efectivitat al modular la variància d’una manera molt significativa. S’ha pogut comprobar en els dos models transformats, però no en el model base on no ha convergit en una solució real i on les altres funcions de modulació de variància no semblen aportar cap millora al model.
Es calcula com abans el coeficient de determinació \(R^2\) amb els mateixos criteris aproximats de càlcul:
M0_RIS_W_list_R20 <- list(Models.sel[c(1:2)],intM0_RIS_W[c(2,4,6,8)])
M0_RIS_W_list_R2 <- ConvList(M0_RIS_W_list_R20)
M0_RIS_W_dfA_R2 <- data.frame(Model = paste("M.Int",str_remove(str_remove(str_remove(names(M0_RIS_W_list_R2), "Model"),"Base"),"RE: ")))
M0_RIS_W_dfA_R2$R2 <- sapply(1:length(M0_RIS_W_list_R2),function(M){
if (class(M0_RIS_W_list_R2[[M]])[1]=="lmm"){
summary(M0_RIS_W_list_R2[[M]])$r.squared
} else if (class(M0_RIS_W_list_R2[[M]])=="lme"){
if(all(paste(attr(M0_RIS_W_list_R2[[M]]$terms,"predvars")[-c(1:2)]) ==
paste(attr(M0_RIS_W_list_R2[[M]]$terms,"term.labels")))){
rsquared(M0_RIS_W_list_R2[[M]])$Conditional
} else {
summary(lm(formula(paste(paste("M0_RIS_W_list_R2[[",M,"]]$data$",
attr(M0_RIS_W_list_R2[[M]]$terms,"variables")[[2]],sep=""),
"~","M0_RIS_W_list_R2[[",M,"]]$fitted"),sep="")))$r.squared
}
}
})
## Warning in paste(attr(M0_RIS_W_list_R2[[M]]$terms, "predvars")[-c(1:2)])
## == : longer object length is not a multiple of shorter object length
## Warning in paste(attr(M0_RIS_W_list_R2[[M]]$terms, "predvars")[-c(1:2)])
## == : longer object length is not a multiple of shorter object length
## Warning in paste(attr(M0_RIS_W_list_R2[[M]]$terms, "predvars")[-c(1:2)])
## == : longer object length is not a multiple of shorter object length
## Warning in paste(attr(M0_RIS_W_list_R2[[M]]$terms, "predvars")[-c(1:2)])
## == : longer object length is not a multiple of shorter object length
## Warning in paste(attr(M0_RIS_W_list_R2[[M]]$terms, "predvars")[-c(1:2)])
## == : longer object length is not a multiple of shorter object length
## Warning in paste(attr(M0_RIS_W_list_R2[[M]]$terms, "predvars")[-c(1:2)])
## == : longer object length is not a multiple of shorter object length
R2_M0W <- ggplot(data=M0_RIS_W_dfA_R2, aes(x=Model, y=R2, fill=Model)) +
geom_bar(stat="identity", position=position_dodge())+
scale_fill_brewer(palette="Set3") + theme_minimal() +
labs(y = "R2 Coeff.") +
theme(axis.text.x = element_text(angle=80,size = 10,hjust=0))
R2_M0W
M1_RIS_W_list_R20 <- list(Models.sel[c(1:2)],intM1_RIS_W[c(2,4,6,8,10,12)])
M1_RIS_W_list_R2 <- ConvList(M1_RIS_W_list_R20)
M1_RIS_W_dfA_R2 <- data.frame(Model = paste("Log.M.Int",str_remove(str_remove(str_remove(names(M1_RIS_W_list_R2), "Model"),"Base"),"RE: ")))
M1_RIS_W_dfA_R2$R2 <- sapply(1:length(M1_RIS_W_list_R2),function(M){
if (class(M1_RIS_W_list_R2[[M]])[1]=="lmm"){
summary(M1_RIS_W_list_R2[[M]])$r.squared
} else if (class(M1_RIS_W_list_R2[[M]])=="lme"){
if(all(paste(attr(M1_RIS_W_list_R2[[M]]$terms,"predvars")[-c(1:2)]) ==
paste(attr(M1_RIS_W_list_R2[[M]]$terms,"term.labels")))){
rsquared(M1_RIS_W_list_R2[[M]])$Conditional
} else {
summary(lm(formula(paste(paste("M1_RIS_W_list_R2[[",M,"]]$data$",
attr(M1_RIS_W_list_R2[[M]]$terms,"variables")[[2]],sep=""),
"~","M1_RIS_W_list_R2[[",M,"]]$fitted"),sep="")))$r.squared
}
}
})
## Warning in paste(attr(M1_RIS_W_list_R2[[M]]$terms, "predvars")[-c(1:2)])
## == : longer object length is not a multiple of shorter object length
## Warning in paste(attr(M1_RIS_W_list_R2[[M]]$terms, "predvars")[-c(1:2)])
## == : longer object length is not a multiple of shorter object length
## Warning in paste(attr(M1_RIS_W_list_R2[[M]]$terms, "predvars")[-c(1:2)])
## == : longer object length is not a multiple of shorter object length
## Warning in paste(attr(M1_RIS_W_list_R2[[M]]$terms, "predvars")[-c(1:2)])
## == : longer object length is not a multiple of shorter object length
## Warning in paste(attr(M1_RIS_W_list_R2[[M]]$terms, "predvars")[-c(1:2)])
## == : longer object length is not a multiple of shorter object length
## Warning in paste(attr(M1_RIS_W_list_R2[[M]]$terms, "predvars")[-c(1:2)])
## == : longer object length is not a multiple of shorter object length
## Warning in paste(attr(M1_RIS_W_list_R2[[M]]$terms, "predvars")[-c(1:2)])
## == : longer object length is not a multiple of shorter object length
## Warning in paste(attr(M1_RIS_W_list_R2[[M]]$terms, "predvars")[-c(1:2)])
## == : longer object length is not a multiple of shorter object length
R2_M1W <- ggplot(data=M1_RIS_W_dfA_R2, aes(x=Model, y=R2, fill=Model)) +
geom_bar(stat="identity", position=position_dodge())+
scale_fill_brewer(palette="Set3") + theme_minimal() +
labs(y = "R2 Coeff.") +
theme(axis.text.x = element_text(angle=80,size = 10,hjust=0))
R2_M1W
M2_RIS_W_list_R20 <- list(Models.sel[c(1:2)],intM2_RIS_W[c(2,4,6,8,12)])
M2_RIS_W_list_R2 <- ConvList(M2_RIS_W_list_R20)
M2_RIS_W_dfA_R2 <- data.frame(Model = paste("BCt.M.Int",str_remove(str_remove(str_remove(names(M2_RIS_W_list_R2), "Model"),"Base"),"RE: ")))
M2_RIS_W_dfA_R2$R2 <- sapply(1:length(M2_RIS_W_list_R2),function(M){
if (class(M2_RIS_W_list_R2[[M]])[1]=="lmm"){
summary(M2_RIS_W_list_R2[[M]])$r.squared
} else if (class(M2_RIS_W_list_R2[[M]])=="lme"){
if(all(paste(attr(M2_RIS_W_list_R2[[M]]$terms,"predvars")[-c(1:2)]) ==
paste(attr(M2_RIS_W_list_R2[[M]]$terms,"term.labels")))){
rsquared(M2_RIS_W_list_R2[[M]])$Conditional
} else {
summary(lm(formula(paste(paste("M2_RIS_W_list_R2[[",M,"]]$data$",
attr(M2_RIS_W_list_R2[[M]]$terms,"variables")[[2]],sep=""),
"~","M2_RIS_W_list_R2[[",M,"]]$fitted"),sep="")))$r.squared
}
}
})
## Warning in paste(attr(M2_RIS_W_list_R2[[M]]$terms, "predvars")[-c(1:2)])
## == : longer object length is not a multiple of shorter object length
## Warning in paste(attr(M2_RIS_W_list_R2[[M]]$terms, "predvars")[-c(1:2)])
## == : longer object length is not a multiple of shorter object length
## Warning in paste(attr(M2_RIS_W_list_R2[[M]]$terms, "predvars")[-c(1:2)])
## == : longer object length is not a multiple of shorter object length
## Warning in paste(attr(M2_RIS_W_list_R2[[M]]$terms, "predvars")[-c(1:2)])
## == : longer object length is not a multiple of shorter object length
## Warning in paste(attr(M2_RIS_W_list_R2[[M]]$terms, "predvars")[-c(1:2)])
## == : longer object length is not a multiple of shorter object length
## Warning in paste(attr(M2_RIS_W_list_R2[[M]]$terms, "predvars")[-c(1:2)])
## == : longer object length is not a multiple of shorter object length
## Warning in paste(attr(M2_RIS_W_list_R2[[M]]$terms, "predvars")[-c(1:2)])
## == : longer object length is not a multiple of shorter object length
R2_M2W <- ggplot(data=M2_RIS_W_dfA_R2, aes(x=Model, y=R2, fill=Model)) +
geom_bar(stat="identity", position=position_dodge())+
scale_fill_brewer(palette="Set3") + theme_minimal() +
labs(y = "R2 Coeff.") +
theme(axis.text.x = element_text(angle=80,size = 10,hjust=0))
R2_M2W
En aquesta ocasió no sembla millorar el model des del punt de vista del coeficient de determinació.
M0_RIS_W_list_base <- M0_RIS_W_list_R2
names(M0_RIS_W_list_base) <- str_remove(str_remove(str_remove(names(M0_RIS_W_list_base), "Model"),"Base"),"RE: ")
#M0_RIS_W_list_R2
par(mfrow=c(1,2))
QQ_M0_list <- lapply(1:length(M0_RIS_W_list_base),function(QQ){
IntRIS_MLE_W.SM0.sdres <- residuals(M0_RIS_W_list_base[[QQ]],type = "pearson")
qqnorm(IntRIS_MLE_W.SM0.sdres,
ylab="Standardized Res.",
xlab="Normal Scores",
main=paste("M0",names(M0_RIS_W_list_base[QQ]),"Standard. Res."))
qqline(IntRIS_MLE_W.SM0.sdres)
})
Shap_M0_df <- data.frame(Model = paste("M2",names(M0_RIS_W_list_base)))
Shap_M0_df$Shapiro_p.valor <- Shap_M0_v <- sapply(1:length(M0_RIS_W_list_base),function(Sh){
shapiro.test(M0_RIS_W_list_base[[Sh]]$residuals)$p
})
Shap_M0_df
M1_RIS_W_list_base <- M1_RIS_W_list_R2
names(M1_RIS_W_list_base) <- str_remove(str_remove(str_remove(names(M1_RIS_W_list_base), "Model"),"Base"),"RE: ")
#M1_RIS_W_list_R2
par(mfrow=c(1,2))
QQ_M1_list <- lapply(1:length(M1_RIS_W_list_base),function(QQ){
IntRIS_MLE_W.SM1.sdres <- residuals(M1_RIS_W_list_base[[QQ]],type = "pearson")
qqnorm(IntRIS_MLE_W.SM1.sdres,
ylab="Standardized Res.",
xlab="Normal Scores",
main=paste("M1",names(M1_RIS_W_list_base[QQ]),"Standard. Res."))
qqline(IntRIS_MLE_W.SM1.sdres)
})
Shap_M1_df <- data.frame(Model = paste("M2",names(M1_RIS_W_list_base)))
Shap_M1_df$Shapiro_p.valor <- Shap_M1_v <- sapply(1:length(M1_RIS_W_list_base),function(Sh){
shapiro.test(M1_RIS_W_list_base[[Sh]]$residuals)$p
})
Shap_M1_df
M2_RIS_W_list_base <- M2_RIS_W_list_R2
names(M2_RIS_W_list_base) <- str_remove(str_remove(str_remove(names(M2_RIS_W_list_base), "Model"),"Base"),"RE: ")
#M2_RIS_W_list_R2
par(mfrow=c(1,2))
QQ_M2_list <- lapply(1:length(M2_RIS_W_list_base),function(QQ){
IntRIS_MLE_W.SM2.sdres <- residuals(M2_RIS_W_list_base[[QQ]],type = "pearson")
qqnorm(IntRIS_MLE_W.SM2.sdres,
ylab="Standardized Res.",
xlab="Normal Scores",
main=paste("M2",names(M2_RIS_W_list_base[QQ]),"Standard. Res."))
qqline(IntRIS_MLE_W.SM2.sdres)
})
Shap_M2_df <- data.frame(Model = paste("M2",names(M2_RIS_W_list_base)))
Shap_M2_df$Shapiro_p.valor <- Shap_M2_v <- sapply(1:length(M2_RIS_W_list_base),function(Sh){
shapiro.test(M2_RIS_W_list_base[[Sh]]$residuals)$p
})
Shap_M2_df
El que es veu en el test Shapiro no seria bo ja que sembla que els models s’allunyen de la normalitat al afegir els efectes aleatoris RIS, tot i així visualment es pot veure que hi ha poques diferències en el QQplot, pel que potser seria bo no prendre’s aquest test de Shapiro com definitiu (en cas que es volgués aprofundir més es podria intentar fer més proves de normalitat).
#Variation of “Levene’s Test”
lev.test <- function(lm,id,df){ #Posar la funció, el nom dels subjectes i el dataset
df$Model.F.Res<- residuals(lm) #extracts the residuals and places them in a new column in our original data table
df$Abs.Model.F.Res <-abs(df$Model.F.Res) #creates a new column with the absolute value of the residuals
df$Model.F.Res2 <- df$Abs.Model.F.Res^2 #squares the absolute values of the residuals to provide the more robust estimate
Levene.Model.F <- lm(paste("Model.F.Res2 ~ ",paste(id),sep=""), data=df) #ANOVA of the squared residuals
return(anova(Levene.Model.F))
}
RP_M0_list <- lapply(1:length(M0_RIS_W_list_base),function(RP){
plot(M0_RIS_W_list_base[[RP]],which=1,main=paste("M0",names(M0_RIS_W_list_base)[RP]))
})
marrangeGrob(RP_M0_list, ncol=2, nrow=1)
LT_M0_list <- lapply(1:length(M0_RIS_W_list_base),function(LT){
lev.test(M0_RIS_W_list_base[[LT]],"Ident",LAKE1)
})
names(LT_M0_list)<-names(M0_RIS_W_list_base)
LT_M0_list
## $M.Int.RIS.MLE
## Analysis of Variance Table
##
## Response: Model.F.Res2
## Df Sum Sq Mean Sq F value Pr(>F)
## Ident 7 3.4565e+19 4.9378e+18 2.4131 0.04203 *
## Residuals 32 6.5481e+19 2.0463e+18
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $M.Int.RIS.REML
## Analysis of Variance Table
##
## Response: Model.F.Res2
## Df Sum Sq Mean Sq F value Pr(>F)
## Ident 7 3.4961e+19 4.9944e+18 2.4079 0.0424 *
## Residuals 32 6.6374e+19 2.0742e+18
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $`IS MLE RIS VarIdent `
## Analysis of Variance Table
##
## Response: Model.F.Res2
## Df Sum Sq Mean Sq F value Pr(>F)
## Ident 7 3.4565e+19 4.9378e+18 2.4131 0.04203 *
## Residuals 32 6.5481e+19 2.0463e+18
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $`IS REML RIS VarIdent `
## Analysis of Variance Table
##
## Response: Model.F.Res2
## Df Sum Sq Mean Sq F value Pr(>F)
## Ident 7 3.4961e+19 4.9944e+18 2.4079 0.0424 *
## Residuals 32 6.6374e+19 2.0742e+18
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $`IS MLE RIS VarExp `
## Analysis of Variance Table
##
## Response: Model.F.Res2
## Df Sum Sq Mean Sq F value Pr(>F)
## Ident 7 9.5405e+21 1.3629e+21 0.7219 0.6544
## Residuals 32 6.0417e+22 1.8880e+21
##
## $`IS REML RIS VarExp `
## Analysis of Variance Table
##
## Response: Model.F.Res2
## Df Sum Sq Mean Sq F value Pr(>F)
## Ident 7 9.6273e+21 1.3753e+21 0.723 0.6535
## Residuals 32 6.0873e+22 1.9023e+21
RP_M1_list <- lapply(1:length(M1_RIS_W_list_base),function(RP){
plot(M1_RIS_W_list_base[[RP]],which=1,main=paste("M1",names(M1_RIS_W_list_base)[RP]))
})
marrangeGrob(RP_M1_list, ncol=2, nrow=1)
LT_M1_list <- lapply(1:length(M1_RIS_W_list_base),function(LT){
lev.test(M1_RIS_W_list_base[[LT]],"Ident",LAKE1)
})
names(LT_M1_list)<-paste("Logt",names(M1_RIS_W_list_base))
LT_M1_list
## $`Logt M.Int.RIS.MLE`
## Analysis of Variance Table
##
## Response: Model.F.Res2
## Df Sum Sq Mean Sq F value Pr(>F)
## Ident 7 3.4565e+19 4.9378e+18 2.4131 0.04203 *
## Residuals 32 6.5481e+19 2.0463e+18
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $`Logt M.Int.RIS.REML`
## Analysis of Variance Table
##
## Response: Model.F.Res2
## Df Sum Sq Mean Sq F value Pr(>F)
## Ident 7 3.4961e+19 4.9944e+18 2.4079 0.0424 *
## Residuals 32 6.6374e+19 2.0742e+18
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $`Logt IS MLE RIS VarIdent `
## Analysis of Variance Table
##
## Response: Model.F.Res2
## Df Sum Sq Mean Sq F value Pr(>F)
## Ident 7 34.297 4.8996 1.1832 0.3399
## Residuals 32 132.516 4.1411
##
## $`Logt IS REML RIS VarIdent `
## Analysis of Variance Table
##
## Response: Model.F.Res2
## Df Sum Sq Mean Sq F value Pr(>F)
## Ident 7 33.593 4.7991 1.1992 0.3315
## Residuals 32 128.061 4.0019
##
## $`Logt IS MLE RIS VarExp `
## Analysis of Variance Table
##
## Response: Model.F.Res2
## Df Sum Sq Mean Sq F value Pr(>F)
## Ident 7 800.5 114.36 0.1408 0.9941
## Residuals 32 25990.3 812.20
##
## $`Logt IS REML RIS VarExp `
## Analysis of Variance Table
##
## Response: Model.F.Res2
## Df Sum Sq Mean Sq F value Pr(>F)
## Ident 7 774.4 110.63 0.1393 0.9943
## Residuals 32 25414.5 794.20
##
## $`Logt IS MLE RIS VarPower `
## Analysis of Variance Table
##
## Response: Model.F.Res2
## Df Sum Sq Mean Sq F value Pr(>F)
## Ident 7 361.5 51.64 0.1461 0.9934
## Residuals 32 11312.4 353.51
##
## $`Logt IS REML RIS VarPower `
## Analysis of Variance Table
##
## Response: Model.F.Res2
## Df Sum Sq Mean Sq F value Pr(>F)
## Ident 7 344.9 49.28 0.1477 0.9931
## Residuals 32 10674.8 333.59
RP_M2_list <- lapply(1:length(M2_RIS_W_list_base),function(RP){
plot(M2_RIS_W_list_base[[RP]],which=1,main=paste("M2",names(M2_RIS_W_list_base)[RP]))
})
marrangeGrob(RP_M2_list, ncol=1, nrow=1)
LT_M2_list <- lapply(1:length(M2_RIS_W_list_base),function(LT){
lev.test(M2_RIS_W_list_base[[LT]],"Ident",LAKE1)
})
names(LT_M2_list)<-paste("BCt",names(M2_RIS_W_list_base))
LT_M2_list
## $`BCt M.Int.RIS.MLE`
## Analysis of Variance Table
##
## Response: Model.F.Res2
## Df Sum Sq Mean Sq F value Pr(>F)
## Ident 7 3.4565e+19 4.9378e+18 2.4131 0.04203 *
## Residuals 32 6.5481e+19 2.0463e+18
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $`BCt M.Int.RIS.REML`
## Analysis of Variance Table
##
## Response: Model.F.Res2
## Df Sum Sq Mean Sq F value Pr(>F)
## Ident 7 3.4961e+19 4.9944e+18 2.4079 0.0424 *
## Residuals 32 6.6374e+19 2.0742e+18
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $`BCt IS MLE RIS VarIdent `
## Analysis of Variance Table
##
## Response: Model.F.Res2
## Df Sum Sq Mean Sq F value Pr(>F)
## Ident 7 3.0061e+10 4294400176 0.7051 0.6677
## Residuals 32 1.9489e+11 6090431011
##
## $`BCt IS REML RIS VarIdent `
## Analysis of Variance Table
##
## Response: Model.F.Res2
## Df Sum Sq Mean Sq F value Pr(>F)
## Ident 7 2.4240e+10 3462840051 0.702 0.6701
## Residuals 32 1.5784e+11 4932635905
##
## $`BCt IS MLE RIS VarExp `
## Analysis of Variance Table
##
## Response: Model.F.Res2
## Df Sum Sq Mean Sq F value Pr(>F)
## Ident 7 3.5176e+10 5.0251e+09 0.1016 0.9979
## Residuals 32 1.5822e+12 4.9443e+10
##
## $`BCt IS REML RIS VarExp `
## Analysis of Variance Table
##
## Response: Model.F.Res2
## Df Sum Sq Mean Sq F value Pr(>F)
## Ident 7 3.5565e+10 5.0807e+09 0.1088 0.9973
## Residuals 32 1.4945e+12 4.6704e+10
##
## $`BCt IS REML RIS VarPower `
## Analysis of Variance Table
##
## Response: Model.F.Res2
## Df Sum Sq Mean Sq F value Pr(>F)
## Ident 7 3.3229e+10 4.7471e+09 0.1999 0.9832
## Residuals 32 7.5997e+11 2.3749e+10
En aquest resultat és on es pot veure realment el que s’ha aconseguit amb les funcions de modulació de variància ja que podem recuperar la homocedasticitat, una de les propietats essencials perquè un model no tingui un biaix important en la desviació.
S’hauria de prendre com a òptim el model que reuneixi les propietats més bones, pel que recuperant els resultats últims es pot fer un petit resum. Donat que a priori no hi ha un factor de pes per cada paràmetre de bondat d’ajust estudiada, l’estratègia a seguir serà tenint en compte els paràmetres Sigma, AIC/BIC, \(R^2\), p-valor Shapiro i p-valor Levene-Test(modificat), seleccionar per cada paràmetre 3 models, combinar la selecció i representar els models resultants:
NomsModelsRISW <- paste(rep(paste(rep(paste("M.int RIS",c("MLE","REML")),4),
rep(c("No.W","VarI","VarE","VarP"),each=2)),3),
rep(c("No.t","Logt","BCt"),each=8))
NomsModelsRISW<-NomsModelsRISW[-which(NomsModelsRISW==c("M.int RIS MLE VarP No.t","M.int RIS REML VarP No.t"))]
NomsModelsRISW<-NomsModelsRISW[-which(NomsModelsRISW==c("M.int RIS MLE VarP BCt"))]
ResumRISW <- rbind(intM0_QI[[1]],intM1_QI[[1]],intM2_QI[[1]])
rownames(ResumRISW) <- NomsModelsRISW
ResumRISW$R2 <- c(M0_RIS_W_dfA_R2$R2,M1_RIS_W_dfA_R2$R2,M2_RIS_W_dfA_R2$R2)
ResumRISW$Shapiro_p.valor <- c(Shap_M0_df$Shapiro_p.valor,Shap_M1_df$Shapiro_p.valor,Shap_M2_df$Shapiro_p.valor)
LT_list <- list(LT_M0_list,LT_M1_list,LT_M2_list)
ResumRISW$Levene_p.valor <- unlist(sapply(1:length(LT_list),function(n){
sapply(1:length(LT_list[[n]]),function(i){LT_list[[n]][[i]]$`Pr(>F)`[1]})
}))
#Es genera una funció per convertir i graficar la taula
Convdf_plot <- function(df){
p1<-c(1)
p2<-c(2:3)
p3<-c(4)
p4<-c(5)
p5<-c(6)
p <- list(p1,p2,p3,p4,p5)
pnames<-c("Sigma","AIC i BIC","R2","Shapiro p.valor","Levene p-valor")
ML_table_p <- lapply(1:length(p),function(pvar){
Table <- data.frame(Model = if(length(p[[pvar]])>1){
rep(rownames(df),ncol(df[,p[[pvar]]]))
} else{
rownames(df)
})
Table$Model <- factor(Table$Model, levels = rownames(df))
Table$ParQ <- unlist(df[,p[[pvar]]])
Table$Par <- rep(colnames(df)[p[[pvar]]],each=nrow(df))
return(Table)
})
ML_plot <- lapply(1:length(ML_table_p),function(plot){
ggplot(data=ML_table_p[[plot]], aes(x=Model, y=ParQ, fill=Par)) +
geom_bar(stat="identity", position=position_dodge())+
scale_fill_brewer(palette="Set3") + theme_minimal() +
labs(title="Optimització Model",y = paste(pnames[plot])) +
theme(axis.text.x = element_text(angle=80,size = 10,hjust=0))
})
return(ML_plot)
}
#Selecció #unique
optorder <- c("F","F","F","T","T","T")
BestFi <- sort(unique(as.vector(sapply(1:ncol(ResumRISW),function(v){
order(ResumRISW[,v],decreasing=optorder[v])[1:2]
}))))
Convdf_plot(ResumRISW[BestFi,])
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
#No es considera el shapiro ni els models amb Sigmes massa grans
ResumRISW2 <- ResumRISW[-which(ResumRISW$Sigma>100),]
BestFi2 <- sort(unique(as.vector(sapply(1:ncol(ResumRISW2[,-4]),function(v){
order(ResumRISW2[,v],decreasing=optorder[v])[1:2]
}))))
Convdf_plot(ResumRISW2[BestFi2,])
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
De manera resumida:
En conjunt sembla que seria ideal utilitzar els dos models transformats amb efectes aleatoris i funció de modulació de variància.
Dels models elegits és interessant veure’n les propietats en quant a matriu de variància-covariància i de correlacions:
Opt <- list(M1_RIS_W_list_base[[8]],M2_RIS_W_list_base[[7]])
names(Opt) <- c("M.Logt.RIS.W-VarP.","M.BCt.RIS.W-V")
lapply(1:length(Opt),function(B){
VarCorr(Opt[[B]])
})
## [[1]]
## Ident = pdLogChol(1 + Tiempo)
## Variance StdDev Corr
## (Intercept) 2.008409e-19 4.481527e-10 (Intr)
## Tiempo 5.694890e-15 7.546449e-08 -1
## Residual 2.047354e-14 1.430858e-07
##
## [[2]]
## Ident = pdLogChol(1 + Tiempo)
## Variance StdDev Corr
## (Intercept) 4.485131e-169 6.697112e-85 (Intr)
## Tiempo 3.587467e-172 1.894061e-86 -0.982
## Residual 5.331961e-169 7.302028e-85
lapply(1:length(Opt),function(B){
getVarCov(Opt[[B]],type="conditional")
})
## [[1]]
## Ident aocampo 002
## Conditional variance covariance matrix
## 1 2 3 4 5
## 1 4.2487 0.00000 0.000000 0.000000 0.000000
## 2 0.0000 0.51613 0.000000 0.000000 0.000000
## 3 0.0000 0.00000 0.029884 0.000000 0.000000
## 4 0.0000 0.00000 0.000000 0.056424 0.000000
## 5 0.0000 0.00000 0.000000 0.000000 0.051378
## Standard Deviations: 2.0612 0.71842 0.17287 0.23754 0.22667
##
## [[2]]
## Ident aocampo 002
## Conditional variance covariance matrix
## 1 2 3 4 5
## 1 174080 0 0 0 0
## 2 0 63839 0 0 0
## 3 0 0 22559 0 0
## 4 0 0 0 18890 0
## 5 0 0 0 0 14827
## Standard Deviations: 417.23 252.66 150.2 137.44 121.77
lapply(1:length(Opt),function(B){
a<-getVarCov(Opt[[B]],type="marginal")
b<-cov2cor(a[[1]])
return(list(a,b))
})
## [[1]]
## [[1]][[1]]
## Ident aocampo 002
## Marginal variance covariance matrix
## 1 2 3 4 5
## 1 4.2487e+00 -4.0562e-16 -8.1144e-16 -1.2173e-15 -1.6231e-15
## 2 -4.0562e-16 5.1613e-01 1.6389e-12 2.4586e-12 3.2782e-12
## 3 -8.1144e-16 1.6389e-12 2.9884e-02 4.9184e-12 6.5581e-12
## 4 -1.2173e-15 2.4586e-12 4.9184e-12 5.6424e-02 9.8379e-12
## 5 -1.6231e-15 3.2782e-12 6.5581e-12 9.8379e-12 5.1378e-02
## Standard Deviations: 2.0612 0.71842 0.17287 0.23754 0.22667
##
## [[1]][[2]]
## 1 2 3 4 5
## 1 1.000000e+00 -2.739135e-16 -2.277243e-15 -2.486133e-15 -3.473972e-15
## 2 -2.739135e-16 1.000000e+00 1.319637e-11 1.440687e-11 2.013128e-11
## 3 -2.277243e-15 1.319637e-11 1.000000e+00 1.197748e-10 1.673661e-10
## 4 -2.486133e-15 1.440687e-11 1.197748e-10 1.000000e+00 1.827184e-10
## 5 -3.473972e-15 2.013128e-11 1.673661e-10 1.827184e-10 1.000000e+00
##
##
## [[2]]
## [[2]][[1]]
## Ident aocampo 002
## Marginal variance covariance matrix
## 1 2 3 4 5
## 1 1.7408e+05 2.9903e-169 1.4954e-169 5.1212e-173 -1.4944e-169
## 2 2.9903e-169 6.3839e+04 1.0337e-169 5.5425e-171 -9.2285e-170
## 3 1.4954e-169 1.0337e-169 2.2559e+04 1.1034e-170 -3.5134e-170
## 4 5.1212e-173 5.5425e-171 1.1034e-170 1.8890e+04 2.2016e-170
## 5 -1.4944e-169 -9.2285e-170 -3.5134e-170 2.2016e-170 1.4827e+04
## Standard Deviations: 417.23 252.66 150.2 137.44 121.77
##
## [[2]][[2]]
## 1 2 3 4
## 1 1.000000e+00 2.836525e-174 2.386256e-174 8.930590e-178
## 2 2.836525e-174 1.000000e+00 2.723912e-174 1.596060e-175
## 3 2.386256e-174 2.723912e-174 1.000000e+00 5.345081e-175
## 4 8.930590e-178 1.596060e-175 5.345081e-175 1.000000e+00
## 5 -2.941374e-174 -2.999580e-174 -1.921086e-174 1.315546e-174
## 5
## 1 -2.941374e-174
## 2 -2.999580e-174
## 3 -1.921086e-174
## 4 1.315546e-174
## 5 1.000000e+00
Tot i obtenir resultats diferents per cada conjunt He les conclusions són semblants:
Tot i haver afegit la funció per modular la covariància no sembla tenir un efecte molt gran tal com es veu en la correlació entre temps, però possiblement pot corregir el model encara que sigui lleugerament. En etapes anterirors també s’ha vist com realment l’efecte aleatori RIS pot haver suposat un canvi positiu pel model.
Els models finals són:
lapply(Opt,summary)
## $`M.Logt.RIS.W-VarP.`
## Linear mixed-effects model fit by REML
## Data: LAKE1
## AIC BIC logLik
## 148.6902 162.9418 -65.34508
##
## Random effects:
## Formula: ~1 + Tiempo | Ident
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 4.481527e-10 (Intr)
## Tiempo 7.546449e-08 -1
## Residual 1.430858e-07
##
## Variance function:
## Structure: Power of variance covariate
## Formula: ~fitted(.)
## Parameter estimates:
## power
## 10.53387
## Fixed effects: intM1$call$formula
## Value Std.Error DF t-value p-value
## (Intercept) 6.911785 1.2209466 29 5.661005 0.0000
## Tiempo -0.070609 0.0293028 29 -2.409634 0.0225
## CD4P -0.092617 0.0422670 29 -2.191237 0.0366
## Tiempo:CD4P 0.002147 0.0010477 29 2.049207 0.0496
## Correlation:
## (Intr) Tiempo CD4P
## Tiempo -0.968
## CD4P -0.981 0.964
## Tiempo:CD4P 0.930 -0.981 -0.964
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.0033954 -0.5022150 -0.1524636 0.5019610 3.3206677
##
## Number of Observations: 40
## Number of Groups: 8
##
## $`M.BCt.RIS.W-V`
## Linear mixed-effects model fit by REML
## Data: LAKE1
## AIC BIC logLik
## 546.2706 560.5222 -264.1353
##
## Random effects:
## Formula: ~1 + Tiempo | Ident
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 6.697112e-85 (Intr)
## Tiempo 1.894061e-86 -0.982
## Residual 7.302028e-85
##
## Variance function:
## Structure: Power of variance covariate
## Formula: ~fitted(.)
## Parameter estimates:
## power
## 22.41349
## Fixed effects: intM2$call$formula
## Value Std.Error DF t-value p-value
## (Intercept) 7946.899 392.2934 29 20.257541 0.0000
## Tiempo -24.123 10.5765 29 -2.280817 0.0301
## CD4P -22.654 16.5012 29 -1.372849 0.1803
## Tiempo:CD4P 0.619 0.4382 29 1.412335 0.1685
## Correlation:
## (Intr) Tiempo CD4P
## Tiempo -0.936
## CD4P -0.973 0.912
## Tiempo:CD4P 0.921 -0.977 -0.945
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.33010638 -0.80255701 -0.04421368 0.69241594 1.86325484
##
## Number of Observations: 40
## Number of Groups: 8