if (!require("pacman")) install.packages("pacman")
## Loading required package: pacman
pacman::p_load("tidyverse", "lme4","readxl","beanplot","knitr", "kableExtra","stargazer", "tableone","sjPlot","sjlabelled","sjmisc","gridExtra","emmeans","sjstats","finalfit")

df=read_xlsx(path = "ajedrez por rival.xlsx",sheet = 1) %>% mutate(Player=as.factor(str_c("player ",JUGADOR)))

Estudio1

Queremos estudiar so T y C de un jugador es afectado por

####
#Estudio2
#Estudio1
#Queremos estudiar so T y C de un jugador es afectado por 
#la diferencia de ELO con el rival

df2=df %>% mutate(diffELO=ELO_RIVAL-ELO)

modelo21=lmer(formula = T ~ diffELO +(1|Player),  data=df2)
summary(modelo21)
## Linear mixed model fit by REML ['lmerMod']
## Formula: T ~ diffELO + (1 | Player)
##    Data: df2
## 
## REML criterion at convergence: 380.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8141 -0.7616  0.1860  0.6131  1.8227 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Player   (Intercept) 42582    206.4   
##  Residual             16178    127.2   
## Number of obs: 30, groups:  Player, 6
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 396.2000    87.3854   4.534
## diffELO       0.7239     0.2189   3.307
## 
## Correlation of Fixed Effects:
##         (Intr)
## diffELO 0.000
modelo21 %>% tab_model(show.se = TRUE) %>% return() %$% 
    knitr %>% 
    asis_output()
  T
Predictors Estimates std. Error CI p
(Intercept) 396.20 87.39 224.93 – 567.47 <0.001
diffELO 0.72 0.22 0.29 – 1.15 0.001
Random Effects
σ2 16178.11
τ00 Player 42581.59
ICC 0.72
N Player 6
Observations 30
Marginal R2 / Conditional R2 0.190 / 0.777

E incluyendo la ronda en la partida, es aún mucho más interesante:

modelo21b=lmer(formula = T ~ RONDA+diffELO +(1|Player),  data=df2)
summary(modelo21)
## Linear mixed model fit by REML ['lmerMod']
## Formula: T ~ diffELO + (1 | Player)
##    Data: df2
## 
## REML criterion at convergence: 380.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8141 -0.7616  0.1860  0.6131  1.8227 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Player   (Intercept) 42582    206.4   
##  Residual             16178    127.2   
## Number of obs: 30, groups:  Player, 6
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 396.2000    87.3854   4.534
## diffELO       0.7239     0.2189   3.307
## 
## Correlation of Fixed Effects:
##         (Intr)
## diffELO 0.000

La forma de enviarlo a una revista es más o menos esta:

modelo21b %>% tab_model(show.se = TRUE) %>% return() %$% 
    knitr %>% 
    asis_output()
  T
Predictors Estimates std. Error CI p
(Intercept) 294.30 98.66 100.93 – 487.67 0.003
RONDA 33.97 15.14 4.28 – 63.65 0.025
diffELO 0.74 0.20 0.34 – 1.13 <0.001
Random Effects
σ2 13762.24
τ00 Player 43263.79
ICC 0.76
N Player 6
Observations 30
Marginal R2 / Conditional R2 0.226 / 0.813

Comparemos los dos modelos por si ayuda a tomar una decisión:

tab_model(modelo21,modelo21b,show.se = FALSE,show.aic = TRUE) %>% return() %$% 
    knitr %>% 
    asis_output()
  T T
Predictors Estimates CI p Estimates CI p
(Intercept) 396.20 224.93 – 567.47 <0.001 294.30 100.93 – 487.67 0.003
diffELO 0.72 0.29 – 1.15 0.001 0.74 0.34 – 1.13 <0.001
RONDA 33.97 4.28 – 63.65 0.025
Random Effects
σ2 16178.11 13762.24
τ00 42581.59 Player 43263.79 Player
ICC 0.72 0.76
N 6 Player 6 Player
Observations 30 30
Marginal R2 / Conditional R2 0.190 / 0.777 0.226 / 0.813
AIC 388.203 378.217

Todo lleva a pensar que es más interesante usar el modelo con la RONDA (mejor R2, mejor estimaciones, AIC(Akaike information Criteria) más bajo).

ggplot(df2,aes(x=diffELO,y=T,col=Player,shape=Player),alpha=0.75)+geom_point()+geom_smooth(method="lm",se=TRUE,alpha=0.05)+theme_classic()+
  coord_cartesian(ylim=c(0,1000))

ggplot(df2,aes(x=diffELO,y=T),alpha=0.75)+geom_point(aes(col=Player))+geom_smooth(method="lm")+theme_bw()

Y ahora queremos una tabla con el valor de la tasa de cambio de T en función de diffELO. Queremos mostrar una tasa de cambio diferente para cada jugador. No es que es lo que me parezca mejor para generalizar a los jugadores del mundo (para eso prefiero el modelo anterior). Pero si lo que pretendemos es describir a nuestros jugadores, lo que sigue está bien:

####
#Estudio2
#Estudio1
#Queremos estudiar so T y C de un jugador es afectado por 
#la diferencia de ELO con el rival

df2=df %>% mutate(diffELO=ELO_RIVAL-ELO)

modelo21b=lm(formula = T ~  Player + diffELO:Player,  data=df2)
summary(modelo21b)
## 
## Call:
## lm(formula = T ~ Player + diffELO:Player, data = df2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -244.56  -79.99   21.05   70.32  214.64 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)  
## (Intercept)             330.4956   133.1421   2.482   0.0231 *
## Playerplayer 2          145.6536   152.6036   0.954   0.3525  
## Playerplayer 3          196.4772   149.3890   1.315   0.2049  
## Playerplayer 4          369.3878   147.5947   2.503   0.0222 *
## Playerplayer 5          -16.1616   149.2922  -0.108   0.9150  
## Playerplayer 6         -167.5945   231.6763  -0.723   0.4787  
## Playerplayer 1:diffELO    0.4535     0.6824   0.665   0.5148  
## Playerplayer 2:diffELO    1.5736     0.5534   2.843   0.0108 *
## Playerplayer 3:diffELO    0.7968     0.5418   1.471   0.1586  
## Playerplayer 4:diffELO    0.5557     0.5353   1.038   0.3130  
## Playerplayer 5:diffELO    0.7888     0.5415   1.457   0.1624  
## Playerplayer 6:diffELO    0.2365     0.8394   0.282   0.7814  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 133.4 on 18 degrees of freedom
## Multiple R-squared:  0.7754, Adjusted R-squared:  0.6382 
## F-statistic:  5.65 on 11 and 18 DF,  p-value: 0.0006427

Para la revista quedamos en que le interesaría las estimaciones de las interacciones (no el resto posiblemente), y el \(R^2\)

modelo21b %>% tab_model(show.se = TRUE) %>% return() %$% 
    knitr %>% 
    asis_output()
  T
Predictors Estimates std. Error CI p
(Intercept) 330.50 133.14 50.77 – 610.22 0.023
Player: player 2 145.65 152.60 -174.95 – 466.26 0.352
Player: player 3 196.48 149.39 -117.38 – 510.33 0.205
Player: player 4 369.39 147.59 59.30 – 679.47 0.022
Player: player 5 -16.16 149.29 -329.81 – 297.49 0.915
Player: player 6 -167.59 231.68 -654.33 – 319.14 0.479
Playerplayer 1:diffELO 0.45 0.68 -0.98 – 1.89 0.515
Playerplayer 2:diffELO 1.57 0.55 0.41 – 2.74 0.011
Playerplayer 3:diffELO 0.80 0.54 -0.34 – 1.94 0.159
Playerplayer 4:diffELO 0.56 0.54 -0.57 – 1.68 0.313
Playerplayer 5:diffELO 0.79 0.54 -0.35 – 1.93 0.162
Playerplayer 6:diffELO 0.24 0.84 -1.53 – 2.00 0.781
Observations 30
R2 / R2 adjusted 0.775 / 0.638

Es interesante ver \(R^2\) ajustada del modelo es bastante alta. Una gran parte de la variabilidad observada en las mediciones de testosterona es explicada por unos niveles y diferencia de elo con rival, que dependen del individuo.

Vamos a ver qué pendiente tiene cada juhgador que se pueda atribuir a la diffELO con el rival, tras ajustar por la RONDA en la que se ha enfrentado a él.

####
#Estudio2
#Estudio1
#Queremos estudiar so T y C de un jugador es afectado por 
#la diferencia de ELO con el rival

df2=df %>% mutate(diffELO=ELO_RIVAL-ELO)

modelo21b=lm(formula = T ~  RONDA+Player + diffELO:Player,  data=df2)
summary(modelo21b)
## 
## Call:
## lm(formula = T ~ RONDA + Player + diffELO:Player, data = df2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -227.17  -39.18   25.13   61.27  131.41 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)   
## (Intercept)             187.1776   132.4408   1.413  0.17562   
## RONDA                    38.5193    15.9991   2.408  0.02770 * 
## Playerplayer 2          177.0537   136.2281   1.300  0.21106   
## Playerplayer 3          230.7029   133.5049   1.728  0.10210   
## Playerplayer 4          406.4312   132.0509   3.078  0.00682 **
## Playerplayer 5            8.9318   133.0687   0.067  0.94727   
## Playerplayer 6         -207.1440   206.5201  -1.003  0.32991   
## Playerplayer 1:diffELO    0.2944     0.6100   0.483  0.63557   
## Playerplayer 2:diffELO    1.6187     0.4921   3.289  0.00433 **
## Playerplayer 3:diffELO    0.9061     0.4836   1.874  0.07827 . 
## Playerplayer 4:diffELO    0.3325     0.4846   0.686  0.50187   
## Playerplayer 5:diffELO    0.8345     0.4815   1.733  0.10119   
## Playerplayer 6:diffELO    0.5504     0.7572   0.727  0.47716   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 118.6 on 17 degrees of freedom
## Multiple R-squared:  0.8325, Adjusted R-squared:  0.7143 
## F-statistic: 7.043 on 12 and 17 DF,  p-value: 0.000182
modelo21b %>% tab_model(show.se = TRUE) %>% return() %$% 
    knitr %>% 
    asis_output()
  T
Predictors Estimates std. Error CI p
(Intercept) 187.18 132.44 -92.25 – 466.60 0.176
RONDA 38.52 16.00 4.76 – 72.27 0.028
Player: player 2 177.05 136.23 -110.36 – 464.47 0.211
Player: player 3 230.70 133.50 -50.97 – 512.37 0.102
Player: player 4 406.43 132.05 127.83 – 685.03 0.007
Player: player 5 8.93 133.07 -271.82 – 289.68 0.947
Player: player 6 -207.14 206.52 -642.86 – 228.58 0.330
Playerplayer 1:diffELO 0.29 0.61 -0.99 – 1.58 0.636
Playerplayer 2:diffELO 1.62 0.49 0.58 – 2.66 0.004
Playerplayer 3:diffELO 0.91 0.48 -0.11 – 1.93 0.078
Playerplayer 4:diffELO 0.33 0.48 -0.69 – 1.36 0.502
Playerplayer 5:diffELO 0.83 0.48 -0.18 – 1.85 0.101
Playerplayer 6:diffELO 0.55 0.76 -1.05 – 2.15 0.477
Observations 30
R2 / R2 adjusted 0.833 / 0.714

Al tener en cuenta la ronda, ha mejorado aún más \(R^2\) y han mejorado las estimaciones de error de los coeficientes relevantes (las interacciones con diffELO).

modelo22=lmer(C ~ diffELO +(1|Player),  data=df2)
## boundary (singular) fit: see ?isSingular
summary(modelo22)
## Linear mixed model fit by REML ['lmerMod']
## Formula: C ~ diffELO + (1 | Player)
##    Data: df2
## 
## REML criterion at convergence: 150.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.2136 -0.6532 -0.1319  0.4076  3.9879 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Player   (Intercept) 0.000    0.000   
##  Residual             6.829    2.613   
## Number of obs: 30, groups:  Player, 6
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 3.942333   0.477122   8.263
## diffELO     0.006245   0.002992   2.087
## 
## Correlation of Fixed Effects:
##         (Intr)
## diffELO 0.000 
## convergence code: 0
## boundary (singular) fit: see ?isSingular

En formato para revista sería:

modelo22 %>% tab_model(show.se = TRUE) %>% return() %$% 
    knitr %>% 
    asis_output()
## Warning: Can't compute random effect variances. Some variance components equal zero.
##   Solution: Respecify random structure!
  C
Predictors Estimates std. Error CI p
(Intercept) 3.94 0.48 3.01 – 4.88 <0.001
diffELO 0.01 0.00 0.00 – 0.01 0.037
Random Effects
σ2 6.83
τ00 Player 0.00
N Player 6
Observations 30
Marginal R2 / Conditional R2 0.131 / NA 
ggplot(df2,aes(x=diffELO,y=C,col=Player,shape=Player),alpha=0.75)+geom_point()+geom_smooth(method="lm",se=TRUE,alpha=0.05)+theme_classic()+
  coord_cartesian(ylim=c(0,20))

Al ver la gráfica, observo la heterogeneidad de las pendientes, y aunque en el modelo aparezca que diffELO tiene efecto significativo, preferiría considerarlo como un falso positivo. No hablaría de la variable C.

ggplot(df2,aes(x=diffELO,y=C,lty=Player),alpha=0.75)+geom_point()+geom_smooth(method="lm",se = T)+theme_bw()

ggplot(df2,aes(x=diffELO,y=C,lty=Player),alpha=0.75)+geom_point()+geom_smooth(method="lm",se = F)+theme_bw()

ggplot(df2,aes(x=diffELO,y=T),alpha=0.75)+geom_point(aes(col=Player))+geom_smooth(method="lm")

####
#Estudio2
#Estudio1
#Queremos estudiar so T y C de un jugador es afectado por 
#la diferencia de ELO con el rival

df2=df %>% mutate(diffELO=ELO_RIVAL-ELO)

modelo21b=lm(formula = C ~ Player + diffELO:Player,  data=df2)
summary(modelo21b)
## 
## Call:
## lm(formula = C ~ Player + diffELO:Player, data = df2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.6392 -0.8665 -0.2109  0.7776  8.7227 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)  
## (Intercept)             2.0841039  2.7382239   0.761   0.4564  
## Playerplayer 2          0.4952781  3.1384728   0.158   0.8764  
## Playerplayer 3          1.1556519  3.0723613   0.376   0.7112  
## Playerplayer 4          0.8912676  3.0354594   0.294   0.7724  
## Playerplayer 5          2.3146183  3.0703706   0.754   0.4607  
## Playerplayer 6         -3.9005228  4.7646977  -0.819   0.4237  
## Playerplayer 1:diffELO -0.0063641  0.0140353  -0.453   0.6557  
## Playerplayer 2:diffELO -0.0083245  0.0113818  -0.731   0.4740  
## Playerplayer 3:diffELO  0.0017526  0.0111429   0.157   0.8768  
## Playerplayer 4:diffELO  0.0149190  0.0110096   1.355   0.1922  
## Playerplayer 5:diffELO  0.0009465  0.0111357   0.085   0.9332  
## Playerplayer 6:diffELO  0.0365598  0.0172626   2.118   0.0484 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.744 on 18 degrees of freedom
## Multiple R-squared:  0.3865, Adjusted R-squared:  0.01153 
## F-statistic: 1.031 on 11 and 18 DF,  p-value: 0.4606
modelo21b %>% tab_model(show.se = TRUE) %>% return() %$% 
    knitr %>% 
    asis_output()
  C
Predictors Estimates std. Error CI p
(Intercept) 2.08 2.74 -3.67 – 7.84 0.456
Player: player 2 0.50 3.14 -6.10 – 7.09 0.876
Player: player 3 1.16 3.07 -5.30 – 7.61 0.711
Player: player 4 0.89 3.04 -5.49 – 7.27 0.772
Player: player 5 2.31 3.07 -4.14 – 8.77 0.461
Player: player 6 -3.90 4.76 -13.91 – 6.11 0.424
Playerplayer 1:diffELO -0.01 0.01 -0.04 – 0.02 0.656
Playerplayer 2:diffELO -0.01 0.01 -0.03 – 0.02 0.474
Playerplayer 3:diffELO 0.00 0.01 -0.02 – 0.03 0.877
Playerplayer 4:diffELO 0.01 0.01 -0.01 – 0.04 0.192
Playerplayer 5:diffELO 0.00 0.01 -0.02 – 0.02 0.933
Playerplayer 6:diffELO 0.04 0.02 0.00 – 0.07 0.048
Observations 30
R2 / R2 adjusted 0.386 / 0.012

El valor de \(R^2\) ajustado y el de los coeficientes de la interaccion convencen aún más de que mejor nohablar mucho de C.