Startle response
Juan Steibel
6/17/2022
Analysis of startle response experiment
Read data
setwd("c:/Documents and Settings/marti/Documents/")
dts<-read_xlsx("Startle test data sheet(statistics)_1.xlsx",skip = 1)## New names:
## * `` -> ...11
## * `` -> ...12
## * `` -> ...13head(dts)## Warning in fansi::strwrap_ctl(x, width = max(width, 0), indent = indent, :
## Encountered a C0 control character, see `?unhandled_ctl`; you can use
## `warn=FALSE` to turn off these warnings.## # 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.
dts$RTS<-as.factor(dts$RTS)
lmm<-lmer(ttr~stress+orientation+RTS+pen_pos+(1|unique_pen),data=dts)
lmm## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: ttr ~ stress + orientation + 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 RTS2
## 85.216 5.844 27.105 1.887 -21.950
## RTS3 pen_pos
## -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 *
## 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 ' ' 1ls_means(lmm,which = NULL,level = 0.95,ddf = c("Kenward-Roger"),pairwise = TRUE)## Least Squares Means table:
##
## Estimate Std. Error df t value lower upper
## orientation1 - orientation2 -27.1048 11.5683 43.0 -2.3430 -50.4347 -3.7749
## orientation1 - orientation3 -1.8869 9.4338 41.2 -0.2000 -20.9363 17.1624
## orientation2 - orientation3 25.2179 11.2524 44.7 2.2411 2.5505 47.8852
## RTS1 - RTS2 21.9500 14.8048 37.1 1.4826 -8.0449 51.9450
## RTS1 - RTS3 39.0390 23.0141 43.3 1.6963 -7.3632 85.4411
## RTS2 - RTS3 17.0890 17.8251 44.8 0.9587 -18.8162 52.9942
## Pr(>|t|)
## orientation1 - orientation2 0.02382 *
## orientation1 - orientation3 0.84245
## orientation2 - orientation3 0.03003 *
## RTS1 - RTS2 0.14662
## RTS1 - RTS3 0.09700 .
## RTS2 - RTS3 0.34285
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Confidence level: 95%
## Degrees of freedom method: Kenward-Rogerls_means(lmm,which = NULL,level = 0.95,ddf = c("Kenward-Roger"),pairwise = FALSE)## Least Squares Means table:
##
## Estimate Std. Error df t value lower upper Pr(>|t|)
## orientation1 44.86895 8.76024 39.7 5.1219 27.15991 62.57800 8.158e-06
## orientation2 71.97375 12.14900 34.5 5.9243 47.29727 96.65022 1.023e-06
## orientation3 46.75588 8.64042 38.3 5.4113 29.26878 64.24299 3.572e-06
## RTS1 74.86253 14.62662 30.7 5.1182 45.01843 104.70663 1.573e-05
## RTS2 52.91251 4.71934 16.3 11.2119 42.92329 62.90174 4.451e-09
## RTS3 35.82354 17.37321 45.0 2.0620 0.83152 70.81557 0.04501
##
## orientation1 ***
## orientation2 ***
## orientation3 ***
## RTS1 ***
## RTS2 ***
## RTS3 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Confidence level: 95%
## Degrees of freedom method: Kenward-RogerThis 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+RTS+pen_pos+(1|unique_pen),data=dts)
lmm2## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: startlemag ~ stress + orientation + 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 RTS2
## 4.81559 -0.07038 0.13259 -0.16748 -0.33252
## RTS3 pen_pos
## -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
## 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+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
## 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+RTS, data = dts, family = categorical(link = #"logit"))
#did not work well
#m1 <- brm(prog2 ~ ses + write, data = ml, family = multinomial(link = "logit"))
#summary(m1)
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.2494indt$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)