Checking Packages
library("lavaan")
library("effects")
library("psych")
library("lme4")
library("lmerTest")
library("sjstats")
library("sjmisc")
library("tidyverse")
library("tinytex")Attaching Resurgene Data and setting working directory
setwd("~/Dropbox/Vanderbilt/Old Courses/Spring 2019/Advanced SCD/Review Paper/Resurgence")
library(readxl)
Resurgence_Data_R <- read_excel("Resurgence_Data_R.xlsx")
View(Resurgence_Data_R)
attach(Resurgence_Data_R)Changing variable type to categorical
Resurgence_Data_R$Study_ID<- as.factor(Resurgence_Data_R$Study_ID)
Resurgence_Data_R$Resurgence<- as.factor(Resurgence_Data_R$Resurgence)
Resurgence_Data_R$SR_Type<- as.factor(Resurgence_Data_R$SR_Type)
vartable(Resurgence_Data_R)##Resurgence Models with random effects There’s an error message with these models, which I think means we don’t have enough data or variance in our data to run random effect models.
R1 = glmer(Resurgence ~ (1|Study_ID),data = Resurgence_Data_R, family = binomial)
summary(R1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Resurgence ~ (1 | Study_ID)
## Data: Resurgence_Data_R
##
## AIC BIC logLik deviance df.resid
## 24.4 26.8 -10.2 20.4 23
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.76038 -0.11772 0.01116 0.01160 0.56806
##
## Random effects:
## Groups Name Variance Std.Dev.
## Study_ID (Intercept) 320.7 17.91
## Number of obs: 25, groups: Study_ID, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.871 3.629 2.444 0.0145 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1R2 = glmer(Resurgence ~ SR_Type + (1|Study_ID),data = Resurgence_Data_R, family = binomial)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Hessian is numerically singular: parameters are not uniquely determinedsummary(R2)## Warning in vcov.merMod(object, use.hessian = use.hessian): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in vcov.merMod(object, correlation = correlation, sigm = sig): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Resurgence ~ SR_Type + (1 | Study_ID)
## Data: Resurgence_Data_R
##
## AIC BIC logLik deviance df.resid
## 26.0 29.6 -10.0 20.0 22
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.76544 -0.13242 0.00000 0.02184 0.56643
##
## Random effects:
## Groups Name Variance Std.Dev.
## Study_ID (Intercept) 222.1 14.9
## Number of obs: 25, groups: Study_ID, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.436e+00 7.098e+00 1.048 0.295
## SR_Type1 2.580e+01 6.476e+06 0.000 1.000
##
## Correlation of Fixed Effects:
## (Intr)
## SR_Type1 0.000
## convergence code: 0
## unable to evaluate scaled gradient
## Hessian is numerically singular: parameters are not uniquely determinedThis is what a random blog says about the error message "Factor z almost always takes value 2, and value 3 is achieved only once. Your response isn’t changing very much, so it’s not easy to estimate the parameters. The random effect also seems to load heavily on factor A. Your data set is not really informative about factor A and the response y. The optimizer knows this.
The Hessian is singular because the optimizer can’t distinguish some of your parameters. This makes sense, given that your data set does not really distinguish them either. So you need a more balanced data set. It would also be nice to have a larger sample size."
##Resurgence Model without random effects
R3 = glm(Resurgence ~ SR_Type ,data = Resurgence_Data_R, family = binomial)
summary(R3)##
## Call:
## glm(formula = Resurgence ~ SR_Type, family = binomial, data = Resurgence_Data_R)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.27352 -1.27352 0.00013 1.08424 1.08424
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.2231 0.4743 0.470 0.638
## SR_Type1 18.3429 2465.3257 0.007 0.994
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 31.343 on 24 degrees of freedom
## Residual deviance: 24.731 on 23 degrees of freedom
## AIC: 28.731
##
## Number of Fisher Scoring iterations: 17It’s not significant (p=.994), but the direction is positive. This means that escape only reinforcement leads to more resurgence than synthesized reinforcement.
Detaching resurgence data and attaching recurrence data
detach(Resurgence_Data_R)
library(readxl)
Recurrence_Data_R <- read_excel("Recurrence_Data_R.xlsx")
View(Recurrence_Data_R)
attach(Recurrence_Data_R)Changing variable type to categorical
Recurrence_Data_R$Study_ID<- as.factor(Recurrence_Data_R$Study_ID)
Recurrence_Data_R$Recurrence<- as.factor(Recurrence_Data_R$Recurrence)
Recurrence_Data_R$PB_Reinforced<- as.factor(Recurrence_Data_R$PB_Reinforced)
vartable(Recurrence_Data_R)##Recurrence Models with random effects This one actually ran with the random effect of study.
Rec1 = glmer(Recurrence ~ (1|Study_ID),data = Recurrence_Data_R, family = binomial)
summary(Rec1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Recurrence ~ (1 | Study_ID)
## Data: Recurrence_Data_R
##
## AIC BIC logLik deviance df.resid
## 32.4 35.4 -14.2 28.4 31
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.78437 -0.16119 0.01856 0.01998 0.95024
##
## Random effects:
## Groups Name Variance Std.Dev.
## Study_ID (Intercept) 150.3 12.26
## Number of obs: 33, groups: Study_ID, 15
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.766 3.159 2.458 0.014 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Rec2 = glmer(Recurrence ~ PB_Reinforced + (1|Study_ID),data = Recurrence_Data_R, family = binomial)
summary(Rec2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Recurrence ~ PB_Reinforced + (1 | Study_ID)
## Data: Recurrence_Data_R
##
## AIC BIC logLik deviance df.resid
## 34.3 38.8 -14.2 28.3 30
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.78313 -0.16069 0.02066 0.02257 0.94360
##
## Random effects:
## Groups Name Variance Std.Dev.
## Study_ID (Intercept) 147.8 12.16
## Number of obs: 33, groups: Study_ID, 15
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.507 3.228 2.326 0.020 *
## PB_Reinforced1 1.181 4.852 0.243 0.808
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## PB_Renfrcd1 -0.213If PB was reinforced (not on extinction) then it lead to more recurence of PB, but the difference is not significant (p=0.80).
##Recurrence Model without random effects
Rec3 = glm(Recurrence ~ PB_Reinforced ,data = Recurrence_Data_R, family = binomial)
summary(Rec3)##
## Call:
## glm(formula = Recurrence ~ PB_Reinforced, family = binomial,
## data = Recurrence_Data_R)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0393 -1.5096 0.8782 0.8782 0.8782
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.7538 0.4287 1.758 0.0787 .
## PB_Reinforced1 1.1921 1.1518 1.035 0.3007
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 38.673 on 32 degrees of freedom
## Residual deviance: 37.372 on 31 degrees of freedom
## AIC: 41.372
##
## Number of Fisher Scoring iterations: 4Same thing as above. If PB was reinforced (not on extinction) then it lead to more recurence of PB, but the difference is not significant.
Detaching recurrence data and attaching function data
detach(Recurrence_Data_R)
library(readxl)
Function_Data_R <- read_excel("Function_Data_R.xlsx")
View(Function_Data_R)
attach(Function_Data_R)Changing variable type to categorical
Function_Data_R$Study_ID<- as.factor(Function_Data_R$Study_ID)
Function_Data_R$Resurgence<- as.factor(Function_Data_R$Resurgence)
Function_Data_R$Mult_maintained<- as.factor(Function_Data_R$Mult_maintained)
vartable(Function_Data_R)##Function Models with random effects This one actually ran too with the random effect of study.
Fun1 = glmer(Resurgence ~ (1|Study_ID),data = Function_Data_R, family = binomial)
summary(Fun1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Resurgence ~ (1 | Study_ID)
## Data: Function_Data_R
##
## AIC BIC logLik deviance df.resid
## 24.4 26.8 -10.2 20.4 23
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.76038 -0.11772 0.01116 0.01160 0.56806
##
## Random effects:
## Groups Name Variance Std.Dev.
## Study_ID (Intercept) 320.7 17.91
## Number of obs: 25, groups: Study_ID, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.871 3.629 2.444 0.0145 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Fun2 = glmer(Resurgence ~ Mult_maintained + (1|Study_ID),data = Function_Data_R, family = binomial)
summary(Fun2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Resurgence ~ Mult_maintained + (1 | Study_ID)
## Data: Function_Data_R
##
## AIC BIC logLik deviance df.resid
## 26.3 30.0 -10.1 20.3 22
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.76370 -0.11851 0.00949 0.01423 0.56699
##
## Random effects:
## Groups Name Variance Std.Dev.
## Study_ID (Intercept) 302.9 17.4
## Number of obs: 25, groups: Study_ID, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.3202 4.3048 1.933 0.0533 .
## Mult_maintained1 0.9685 4.2525 0.228 0.8198
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Mult_mntnd1 -0.448Again not significant (p=0.819), but multiply maintained coded as 1 had more resurgence than whatever the 0 coding was.
##Function Model without random effects
Fun3 = glm(Resurgence ~ Mult_maintained ,data = Function_Data_R, family = binomial)
summary(Fun3)##
## Call:
## glm(formula = Resurgence ~ Mult_maintained, family = binomial,
## data = Function_Data_R)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8465 -1.3018 0.6335 1.0579 1.0579
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.2877 0.5401 0.533 0.594
## Mult_maintained1 1.2164 0.9501 1.280 0.200
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 31.343 on 24 degrees of freedom
## Residual deviance: 29.552 on 23 degrees of freedom
## AIC: 33.552
##
## Number of Fisher Scoring iterations: 4Not significant (p=0.20), but same direction as above. Multiply maintained has more resurgence.