Nociception
2025-06-24
data<-Behavioral.Trials.All.DataDoes hot predict tailfip?
#Binary Data
data$tailflip_binary <- ifelse(data$tailflip > 0, 1, 0)
#tailflip in trial? yes or no
tailfliphot <- glmer(tailflip_binary ~ hot + (1 | subject), family = binomial, data = data)
summary(tailfliphot)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: tailflip_binary ~ hot + (1 | subject)
## Data: data
##
## AIC BIC logLik -2*log(L) df.resid
## 111.1 119.9 -52.6 105.1 135
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.07949 -0.29656 0.06827 0.22654 1.69858
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 3.783 1.945
## Number of obs: 138, groups: subject, 23
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.0675 0.6604 -3.131 0.00174 **
## hoty 5.4006 1.0815 4.994 5.92e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## hoty -0.662model_null <- glmer(tailflip_binary ~ 1 + (1 | subject), family = binomial, data = data)
anova(model_null, tailfliphot)## Data: data
## Models:
## model_null: tailflip_binary ~ 1 + (1 | subject)
## tailfliphot: tailflip_binary ~ hot + (1 | subject)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## model_null 2 193.35 199.21 -94.676 189.35
## tailfliphot 3 111.14 119.92 -52.571 105.14 84.209 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Count data
#tailflip count
tailfliphot <-glmer(tailflip ~ hot + (1 | subject), family = poisson, data = data)
summary(tailfliphot)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: tailflip ~ hot + (1 | subject)
## Data: data
##
## AIC BIC logLik -2*log(L) df.resid
## 441.8 450.6 -217.9 435.8 135
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9539 -0.6873 -0.4916 0.4801 4.9277
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 0.2301 0.4796
## Number of obs: 138, groups: subject, 23
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.9476 0.2091 -4.532 5.85e-06 ***
## hoty 2.0237 0.1922 10.531 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## hoty -0.812model_null <- glmer(tailflip ~ 1 + (1 | subject), family = poisson, data = data)
anova(model_null, tailfliphot)## Data: data
## Models:
## model_null: tailflip ~ 1 + (1 | subject)
## tailfliphot: tailflip ~ hot + (1 | subject)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## model_null 2 610.82 616.68 -303.41 606.82
## tailfliphot 3 441.77 450.55 -217.89 435.77 171.05 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#hot predicts tailflip significantly, both for binary (tailflip or no tailflip) and for counts.
#Trial had no effect - took it out of the model.
Count data
library(scales) # for hue_pal()
# Get all unique subjects in the dataset (regardless of treatment)
unique_subjects <- unique(data$subject)
# Generate a large enough palette (here using ggplot2's hue_pal across subjects)
subject_colors <- setNames(hue_pal()(length(unique_subjects)), unique_subjects)
df_summary <- data %>%
group_by(hot) %>%
summarise(
mean_tailflip = mean(tailflip, na.rm = TRUE),
se_tailflip = sd(tailflip, na.rm = TRUE) / sqrt(n()),
.groups = "drop"
)
p <- ggplot(data, aes(x = hot, y = tailflip, color = subject)) +
geom_jitter(width = 0.3, height= 0, alpha = 0.7, size = 2) +
scale_color_manual(values = subject_colors) +
scale_x_discrete(labels = c("n" = "control", "y" = "hot")) +
geom_point(data = df_summary, aes(x = hot, y = mean_tailflip),
color = "black", size = 3, inherit.aes = FALSE) +
geom_errorbar(data = df_summary, aes(x = hot, ymin = mean_tailflip - se_tailflip, ymax = mean_tailflip + se_tailflip),
width = 0.2, color = "black", inherit.aes = FALSE) +
labs(title = "Tailflip Counts by Subject with Mean ± SE",
x = "Treatment",
y = "Number of Tail Flips") +
theme_minimal() +
guides(color = guide_legend(title = "Subject"))
p +
scale_y_continuous(expand = expansion(mult = c(0, 0.05))) +
annotate("text", x = 2.36, y = -Inf, label = "\nz = 10.531\np < 2e-16\nIRR ≈ 7.6",
hjust = 0, vjust = -0.1, size = 4) +
guides(color = "none")
Binary data
library(scales) # for hue_pal()
# Make sure tailflip_binary is 0/1 in data
data <- data %>%
mutate(tailflip_binary = ifelse(tailflip > 0, 1, 0))
# Get all unique subjects in the dataset (regardless of treatment)
unique_subjects <- unique(data$subject)
# Generate a large enough palette (here using ggplot2's hue_pal across subjects)
subject_colors <- setNames(hue_pal()(length(unique_subjects)), unique_subjects)
p_binary_alt <- ggplot(data, aes(x = hot, y = tailflip_binary, color = subject)) +
geom_jitter(width = 0.3, height = 0.05, alpha = 0.7, size = 2) +
scale_color_manual(values = subject_colors) +
stat_summary(aes(x = hot, y = tailflip_binary), fun = mean, geom = "point",
size = 3, color = "black", inherit.aes = FALSE) +
stat_summary(aes(x = hot, y = tailflip_binary), fun.data = mean_cl_normal, geom = "errorbar",
width = 0.2, color = "black", inherit.aes = FALSE)+
scale_x_discrete(labels = c("n" = "control", "y" = "hot")) +
scale_y_continuous(breaks = c(0, 1), labels = c("No Tailflip", "Tailflip"),
limits = c(-0.1, 1.0)) +
labs(title = "Binary Tail Flip Outcomes by Subject",
x = "Treatment",
y = "Tailflip - Yes (1) or No (0)") +
theme_minimal() +
guides(color = "none") + theme_minimal()
p_binary_alt +
scale_y_continuous(expand = expansion(mult = c(0, 0.05))) +
annotate("text", x = 2.30, y = -Inf, label = "z = 4.994\np < 5.92e-07\nOR ≈ 222",
hjust = 0, vjust = -0.1, size = 4) + guides(color = "none") + theme_minimal()## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
# Make sure tailflip_binary is 0/1 in data
data <- data %>%
mutate(tailflip_binary = ifelse(tailflip > 0, 1, 0))
# Summarize proportion of trials with at least one tailflip
binary_summary <- data %>%
group_by(hot) %>%
summarise(
prop_tailflip = mean(tailflip_binary, na.rm = TRUE),
se_prop = sqrt(prop_tailflip * (1 - prop_tailflip) / n()),
.groups = "drop"
)
p_binary <- ggplot(binary_summary, aes(x = hot, y = prop_tailflip, fill = hot)) +
geom_col(alpha = 0.8, width = 0.5) +
geom_errorbar(aes(ymin = prop_tailflip - se_prop, ymax = prop_tailflip + se_prop),
width = 0.2, color = "black") +
scale_fill_manual(values = c("n" = "gray60", "y" = "red2")) +
scale_x_discrete(labels = c("n" = "control", "y" = "hot")) +
labs(title = "Probability of Tailflip Occurrence by Treatment",
x = NULL,
y = "Proportion with Tailflip") +
theme_minimal() +
guides(fill = "none") # remove legend if desired
print(p_binary)
#residuals
# Tailflip count model (Poisson)
model_count <- glmer(tailflip ~ hot + (1 | subject), data = data, family = poisson(link = "log"))
# Tailflip binary model (Poisson GLMM)
model_binary <- glmer(tailflip_binary ~ hot + (1 | subject), data = data, family = poisson(link = "log"))## boundary (singular) fit: see help('isSingular')# Pearson residuals
pearson_resid_count <- resid(model_count, type = "pearson")
fitted_count <- fitted(model_count)
# Residual vs. fitted plot
plot(fitted_count, pearson_resid_count,
xlab = "Fitted values", ylab = "Pearson residuals",
main = "Residuals vs Fitted (Count Model)")
abline(h = 0, col = "red")
# QQ plot of residuals
qqnorm(pearson_resid_count, main = "QQ Plot of Residuals (Count Model)")
qqline(pearson_resid_count, col = "red")