Setup for data analysis

packages activation

library(dplyr)
library(readxl)
library(Matrix)
library(lme4)
library(ordinal)

Data Handling

dts<-read_xlsx("Startle test data sheet(statistics)_Dr.Steibel.xlsx",skip = 1)
## New names:
## • `` -> `...11`
## • `` -> `...12`
## • `` -> `...13`
head(dts)
## # A tibble: 6 × 13
##   Replic…¹ `Pen#` `Pig#` SR or…² Relax…³ Orien…⁴ Start…⁵ Outco…⁶ Start…⁷ Time …⁸
##      <dbl> <chr>   <dbl> <chr>     <dbl>   <dbl> <chr>   <chr>     <dbl>   <dbl>
## 1        1 ML5         3 SV            3       1 SI      LD            1      17
## 2        1 ML6         2 SV            3       3 ST      LD            1      15
## 3        1 ML6         3 SR            2       3 RA      FD            3      13
## 4        1 ML6         5 SV            2       3 RA      FD            3      31
## 5        1 ML6         6 SR            2       1 RA      LD            3      26
## 6        1 MR6         7 SR            2       3 RA      NIC           3      57
## # … with 3 more variables: ...11 <lgl>, ...12 <chr>, ...13 <chr>, and
## #   abbreviated variable names ¹​`Replicate#`, ²​`SR or SV\r\n(Stress)`,
## #   ³​`Relaxed-\r\ntense\r\nscore(1-3)\r\n(RTS)`,
## #   ⁴​`Orientation\r\nscore\r\n(1-3)\r\n(Orientation)`,
## #   ⁵​`Startle\r\nresponse\r\n(Startleresp)`,
## #   ⁶​`Outcome behavior\r\nafter end freeze\r\n(OAF)`,
## #   ⁷​`Startle magnitude\r\nscore\r\n(0-4)\r\n(Startlemag)`, …
dts<-dts[,1:10]
colnames(dts)<-c("replicate","pen","pig","stress","RTS","orientation","startleresp","OAF","startlemag","ttr")
dts<-mutate(dts,unique_pen=as.factor(paste(replicate,pen)),orientation=as.factor(orientation))
dts$pen_pos<-as.numeric(substr(dts$pen,3,3))
head(dts)
## # A tibble: 6 × 12
##   replicate pen     pig stress   RTS orien…¹ start…² OAF   start…³   ttr uniqu…⁴
##       <dbl> <chr> <dbl> <chr>  <dbl> <fct>   <chr>   <chr>   <dbl> <dbl> <fct>  
## 1         1 ML5       3 SV         3 1       SI      LD          1    17 1 ML5  
## 2         1 ML6       2 SV         3 3       ST      LD          1    15 1 ML6  
## 3         1 ML6       3 SR         2 3       RA      FD          3    13 1 ML6  
## 4         1 ML6       5 SV         2 3       RA      FD          3    31 1 ML6  
## 5         1 ML6       6 SR         2 1       RA      LD          3    26 1 ML6  
## 6         1 MR6       7 SR         2 3       RA      NIC         3    57 1 MR6  
## # … with 1 more variable: pen_pos <dbl>, and abbreviated variable names
## #   ¹​orientation, ²​startleresp, ³​startlemag, ⁴​unique_pen

Box plot checking for finding difference on ttr between stress levels

### difference in stress levels - no difference betwwen SR and SV
boxplot(dts$ttr~dts$stress)

aov1<-aov(dts$ttr~dts$stress)
summary(aov1)
##             Df Sum Sq Mean Sq F value Pr(>F)
## dts$stress   1     94    94.2   0.099  0.754
## Residuals   50  47646   952.9

Analysis of Time to response (continuous variable)

lmm<-lmer(ttr~stress+orientation+as.factor(RTS)+pen_pos+(1|unique_pen),data=dts)
lmm
## Linear mixed model fit by REML ['lmerMod']
## Formula: ttr ~ stress + orientation + as.factor(RTS) + pen_pos + (1 |  
##     unique_pen)
##    Data: dts
## REML criterion at convergence: 442.8785
## Random effects:
##  Groups     Name        Std.Dev.
##  unique_pen (Intercept)  8.813  
##  Residual               26.081  
## Number of obs: 52, groups:  unique_pen, 21
## Fixed Effects:
##     (Intercept)         stressSV     orientation2     orientation3  
##          85.216            5.844           27.105            1.887  
## as.factor(RTS)2  as.factor(RTS)3          pen_pos  
##         -21.950          -39.039           -4.771
summary(lmm)
## Linear mixed model fit by REML ['lmerMod']
## Formula: ttr ~ stress + orientation + as.factor(RTS) + pen_pos + (1 |  
##     unique_pen)
##    Data: dts
## 
## REML criterion at convergence: 442.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0787 -0.7016 -0.2009  0.6386  1.9318 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  unique_pen (Intercept)  77.67    8.813  
##  Residual               680.24   26.081  
## Number of obs: 52, groups:  unique_pen, 21
## 
## Fixed effects:
##                 Estimate Std. Error t value
## (Intercept)       85.216     19.727   4.320
## stressSV           5.844      8.189   0.714
## orientation2      27.105     11.031   2.457
## orientation3       1.887      9.038   0.209
## as.factor(RTS)2  -21.950     14.012  -1.566
## as.factor(RTS)3  -39.039     22.061  -1.770
## pen_pos           -4.771      3.361  -1.420
## 
## Correlation of Fixed Effects:
##             (Intr) strsSV orntt2 orntt3 a.(RTS)2 a.(RTS)3
## stressSV    -0.288                                       
## orientatin2 -0.211  0.273                                
## orientatin3 -0.293  0.159  0.436                         
## as.fc(RTS)2 -0.376 -0.159 -0.250 -0.057                  
## as.fc(RTS)3 -0.100 -0.312 -0.155 -0.157  0.620           
## pen_pos     -0.723  0.170  0.157  0.109 -0.241   -0.285
anova(lmm)
## Analysis of Variance Table
##                npar Sum Sq Mean Sq F value
## stress            1   25.5   25.55  0.0376
## orientation       2 4518.8 2259.38  3.3214
## as.factor(RTS)    2 3905.2 1952.58  2.8704
## pen_pos           1 1370.8 1370.78  2.0151

Magnitude of the startle response is an ordinal variable

(1) reponse follows the Gaussian assumption
lmm2<-lmer(startlemag~stress+orientation+as.factor(RTS)+pen_pos+(1|unique_pen),data=dts)
lmm2
## Linear mixed model fit by REML ['lmerMod']
## Formula: startlemag ~ stress + orientation + as.factor(RTS) + pen_pos +  
##     (1 | unique_pen)
##    Data: dts
## REML criterion at convergence: 123.4041
## Random effects:
##  Groups     Name        Std.Dev.
##  unique_pen (Intercept) 0.3542  
##  Residual               0.7207  
## Number of obs: 52, groups:  unique_pen, 21
## Fixed Effects:
##     (Intercept)         stressSV     orientation2     orientation3  
##         4.81559         -0.07038          0.13259         -0.16748  
## as.factor(RTS)2  as.factor(RTS)3          pen_pos  
##        -0.33252         -1.70314         -0.29285
summary(lmm2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: startlemag ~ stress + orientation + as.factor(RTS) + pen_pos +  
##     (1 | unique_pen)
##    Data: dts
## 
## REML criterion at convergence: 123.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9310 -0.6533  0.1134  0.5695  1.9619 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  unique_pen (Intercept) 0.1254   0.3542  
##  Residual               0.5194   0.7207  
## Number of obs: 52, groups:  unique_pen, 21
## 
## Fixed effects:
##                 Estimate Std. Error t value
## (Intercept)      4.81559    0.58636   8.213
## stressSV        -0.07038    0.23278  -0.302
## orientation2     0.13259    0.31963   0.415
## orientation3    -0.16748    0.26333  -0.636
## as.factor(RTS)2 -0.33252    0.41112  -0.809
## as.factor(RTS)3 -1.70314    0.63910  -2.665
## pen_pos         -0.29285    0.10198  -2.872
## 
## Correlation of Fixed Effects:
##             (Intr) strsSV orntt2 orntt3 a.(RTS)2 a.(RTS)3
## stressSV    -0.259                                       
## orientatin2 -0.193  0.258                                
## orientatin3 -0.281  0.124  0.446                         
## as.fc(RTS)2 -0.366 -0.186 -0.260 -0.032                  
## as.fc(RTS)3 -0.111 -0.331 -0.153 -0.131  0.632           
## pen_pos     -0.732  0.173  0.151  0.088 -0.242   -0.275
anova(lmm2)
## Analysis of Variance Table
##                npar Sum Sq Mean Sq F value
## stress            1 0.6964  0.6964  1.3408
## orientation       2 1.0372  0.5186  0.9985
## as.factor(RTS)    2 7.0309  3.5154  6.7686
## pen_pos           1 4.2828  4.2828  8.2460
(2) ordinal regression without any random effects
dts$startlemag<-factor(dts$startlemag,ordered = TRUE)
lmm3 = clm(startlemag ~ stress+orientation+as.factor(RTS)+pen_pos,data=dts)
anova(lmm3, type="III")
## Type III Analysis of Deviance Table with Wald chi-square tests
## 
##                Df  Chisq Pr(>Chisq)   
## stress          1 0.0265   0.870581   
## orientation     2 2.4672   0.291239   
## as.factor(RTS)  2 6.1596   0.045968 * 
## pen_pos         1 8.7219   0.003144 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(lmm3)
## formula: startlemag ~ stress + orientation + as.factor(RTS) + pen_pos
## data:    dts
## 
##  link  threshold nobs logLik AIC    niter max.grad cond.H 
##  logit flexible  52   -53.38 124.77 5(0)  5.85e-09 3.7e+03
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)   
## stressSV         -0.1004     0.6162  -0.163  0.87058   
## orientation2      1.0193     0.8195   1.244  0.21358   
## orientation3     -0.2185     0.6175  -0.354  0.72340   
## as.factor(RTS)2  -1.4023     1.0500  -1.335  0.18173   
## as.factor(RTS)3  -4.3629     1.7632  -2.474  0.01334 * 
## pen_pos          -0.7544     0.2555  -2.953  0.00314 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Threshold coefficients:
##     Estimate Std. Error z value
## 1|2   -8.078      1.925  -4.196
## 2|3   -6.281      1.793  -3.503
## 3|4   -4.035      1.647  -2.450
(3) ordinal with a random effect
lmm4 = clmm2(startlemag ~ stress+orientation+as.factor(RTS)+pen_pos,random=unique_pen,data=dts,Hess=TRUE)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
summary(lmm4)
## Cumulative Link Mixed Model fitted with the Laplace approximation
## 
## Call:
## clmm2(location = startlemag ~ stress + orientation + as.factor(RTS) + 
##     pen_pos, random = unique_pen, data = dts, Hess = TRUE)
## 
## Random effects:
##                 Var   Std.Dev
## unique_pen 0.403829 0.6354754
## 
## Location coefficients:
##                 Estimate Std. Error z value Pr(>|z|) 
## stressSV        -0.0999   0.6472    -0.1544 0.8772633
## orientation2     0.7377   0.9968     0.7401 0.4592364
## orientation3    -0.4037   0.7470    -0.5405 0.5888654
## as.factor(RTS)2 -1.3281   1.1658    -1.1392 0.2546121
## as.factor(RTS)3 -4.5643   1.9384    -2.3546 0.0185430
## pen_pos         -0.8185   0.3031    -2.7005 0.0069229
## 
## No scale coefficients
## 
## Threshold coefficients:
##     Estimate Std. Error z value
## 1|2 -8.6330   2.2927    -3.7654
## 2|3 -6.7422   2.1019    -3.2077
## 3|4 -4.3221   1.8719    -2.3090
## 
## log-likelihood: -53.17019 
## AIC: 126.3404 
## Condition number of Hessian: 3780.579
lmm4
## Cumulative Link Mixed Model fitted with the Laplace approximation
## 
## Call:
## clmm2(location = startlemag ~ stress + orientation + as.factor(RTS) + 
##     pen_pos, random = unique_pen, data = dts, Hess = TRUE)
## 
## Random effects:
##                 Var   Std.Dev
## unique_pen 0.403829 0.6354754
## 
## Location coefficients:
##        stressSV    orientation2    orientation3 as.factor(RTS)2 as.factor(RTS)3 
##     -0.09994919      0.73771158     -0.40372865     -1.32806410     -4.56425438 
##         pen_pos 
##     -0.81852824 
## 
## No Scale coefficients
## 
## Threshold coefficients:
##       1|2       2|3       3|4 
## -8.633008 -6.742234 -4.322134 
## 
## log-likelihood: -53.17019 
## AIC: 126.3404

23.01.07 update

stress vs. RTS - chisq test or fisher exact test (category-difference)
ct<-table(as.factor(dts$RTS),dts$stress)
t(ct)
##     
##       1  2  3
##   SR  3 23  0
##   SV  2 21  3
### expected value ratio check - <5 cells a lot
## ex<-chisq.test(ct)
## ex$expected
fisher.test(t(ct))
## 
##  Fisher's Exact Test for Count Data
## 
## data:  t(ct)
## p-value = 0.2736
## alternative hypothesis: two.sided