Startle response
Juan Steibel
6/17/2022
Analysis of startle response experiment
Read data
setwd("C:/Users/marti/OneDrive/Documents")
dts<-read_xlsx("Startle test data sheet(statistics)_1.xlsx",skip = 1)## New names:
## * `` -> ...11
## * `` -> ...12
## * `` -> ...13head(dts)## # A tibble: 6 x 13
## `Replicate#` `Pen#` `Pig#` `SR or SV` `Relaxed-\r\ntense\~ `Orientation\r\nsc~
## <dbl> <chr> <dbl> <chr> <dbl> <dbl>
## 1 1 ML5 3 SV 3 1
## 2 1 ML6 2 SV 3 3
## 3 1 ML6 3 SR 2 3
## 4 1 ML6 5 SV 2 3
## 5 1 ML6 6 SR 2 1
## 6 1 MR6 7 SR 2 3
## # ... with 7 more variables: Startle
## response <chr>,
## # Outcome behavior
## after end freeze <chr>,
## # Startle magnitude
## score
## (0-4) <dbl>, Time taken to resume
## (sec) <dbl>,
## # ...11 <lgl>, ...12 <chr>, ...13 <chr>#one row is used as title of spread sheet.
#avoid this in the future
#use skip for now
#column names too long, replace
#extra sumaries in additional columns cause problems, delete
#for now: use code to clean things up
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))
head(dts)## # A tibble: 6 x 11
## replicate pen pig stress RTS orientation startleresp OAF startlemag
## <dbl> <chr> <dbl> <chr> <dbl> <fct> <chr> <chr> <dbl>
## 1 1 ML5 3 SV 3 1 SI LD 1
## 2 1 ML6 2 SV 3 3 ST LD 1
## 3 1 ML6 3 SR 2 3 RA FD 3
## 4 1 ML6 5 SV 2 3 RA FD 3
## 5 1 ML6 6 SR 2 1 RA LD 3
## 6 1 MR6 7 SR 2 3 RA NIC 3
## # ... with 2 more variables: ttr <dbl>, unique_pen <fct>dts$pen_pos<-as.numeric(substr(dts$pen,3,3))Summary descriptives
These tables and graphic are included here for discussion. Do we see any confounding pattern? This does not need to be included in your m&m or results. I will comment it
#table(dts$pen,dts$stress)
#table(dts$stress,dts$RTS)
#table(dts$stress,dts$orientation)
#table(dts$stress,dts$OAF)
#table(dts$stress,dts$startleresp)
#table(dts$stress,dts$startlemag)
#table(dts$pen,dts$pen_pos)
#table(dts$pig,dts$stress,paste(dts$replicate,dts$pen,dts$pig))
boxplot(dts$ttr~dts$stress)
Analysis of Time to response (continuous variable)
This linear mixed effect model with the stress group (variable of interes=t), plus orientation, pen position and RTS (relaxed/tense status) as fixed covariates/ factors and pen as a random effect. Pen is represented by unique_pen and it represents the combination of replicate and pen.
lmm<-lmer(ttr~stress+orientation+as.factor(RTS)+pen_pos+(1|unique_pen),data=dts)
lmm## Linear mixed model fit by REML ['lmerModLmerTest']
## 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.771anova(lmm) ## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## stress 346.4 346.39 1 43.753 0.5092 0.47927
## orientation 4733.1 2366.56 2 43.135 3.4790 0.03974 *
## as.factor(RTS) 2373.5 1186.75 2 43.168 1.7446 0.18681
## pen_pos 1370.8 1370.78 1 24.025 2.0151 0.16858
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1This did not produce a significative difference in terms of stress/resilient status. We can discuss adding other correction factors
Magnitude of the startle response is an ordinal variable
We tried several models. First we fit a linear mixed model assuming the response as a Gaussian variable.
lmm2<-lmer(startlemag~stress+orientation+as.factor(RTS)+pen_pos+(1|unique_pen),data=dts)
lmm2## Linear mixed model fit by REML ['lmerModLmerTest']
## 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.29285anova(lmm2) ## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## stress 0.0475 0.0475 1 43.034 0.0914 0.763855
## orientation 0.5264 0.2632 2 44.198 0.5068 0.605895
## as.factor(RTS) 4.3526 2.1763 2 44.111 4.1902 0.021572 *
## pen_pos 4.2828 4.2828 1 22.473 8.2460 0.008747 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Second, we tried an ordinal regression WITHOUT 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 ' ' 1an ordinal regression WITH random effects! I recommend keeping these.
Finally, an ordinal logistic regression model with same random and fixed effects as ttr was fit. This is the model that should be replaced
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 occurredlmm4## 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.0999492 0.7377116 -0.4037286 -1.3280641 -4.5642544
## pen_pos
## -0.8185282
##
## No Scale coefficients
##
## Threshold coefficients:
## 1|2 2|3 3|4
## -8.633008 -6.742234 -4.322134
##
## log-likelihood: -53.17019
## AIC: 126.3404summary(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.87726369
## orientation2 0.7377 0.9968 0.7401 0.45923977
## orientation3 -0.4037 0.7470 -0.5405 0.58886889
## as.factor(RTS)2 -1.3281 1.1658 -1.1392 0.25461906
## as.factor(RTS)3 -4.5643 1.9385 -2.3546 0.01854400
## pen_pos -0.8185 0.3031 -2.7004 0.00692486
##
## No scale coefficients
##
## Threshold coefficients:
## Estimate Std. Error z value
## 1|2 -8.6330 2.2929 -3.7652
## 2|3 -6.7422 2.1020 -3.2075
## 3|4 -4.3221 1.8720 -2.3088
##
## log-likelihood: -53.17019
## AIC: 126.3404
## Condition number of Hessian: 3780.96Notice that the results vary minimally with different model assumptions.
Analysis of OAF
This is even trickier as it’s a multinomial response model. I tried two linear models and they did not converge. (all commented) Then tried a chi.square test of independence. Even though it was not significant, I added the results. Notice the expected vs. observed. Given the extreme frequencies, the p-value was obtained using MCMC resampling, thus, no d.f. are provided. This is an exact and not an approximate test.
# library(brms)
#dts$OAF<-as.factor(dts$OAF)
#m2 <- brm(OAF ~ stress+orientation+as.factor(RTS), data = dts, family = categorical(link #= "logit"))
#summary(m2)
##did not work well
#m3 <- brm(OAF ~ stress+orientation+as.factor(RTS), data = dts, family = multinomial(link #= "logit"))
#summary(m3)
indt<-chisq.test(dts$OAF,dts$stress,simulate.p.value = TRUE)
print(indt)##
## Pearson's Chi-squared test with simulated p-value (based on 2000
## replicates)
##
## data: dts$OAF and dts$stress
## X-squared = 7.2466, df = NA, p-value = 0.2454indt$observed## dts$stress
## dts$OAF SR SV
## CP 1 0
## E 1 0
## EX 7 10
## FD 4 8
## L 1 0
## LD 4 5
## NIC 8 3indt$expected## dts$stress
## dts$OAF SR SV
## CP 0.5 0.5
## E 0.5 0.5
## EX 8.5 8.5
## FD 6.0 6.0
## L 0.5 0.5
## LD 4.5 4.5
## NIC 5.5 5.5studentized residuals for time to response
For the continuous variable we can obtain the studentized residuals. These are the residual value (observed-predicted) divided by their estimated (model-based) standard deviation. They are useful to: 1) diagnose lack of normality, 2)flag observations for follow up.
In this case, the diagnostic plots does not point at any issues.
qqnorm(rstudent(lmm))
abline(0,1,col="red",lwd=2)