Using this code block in fresh install of R,
if (!require('readr')) install.packages('readr')
data <- as.data.frame(readr::read_delim("/Users/{le me}/Desktop/{path}/Histology Data for R.csv", delim = ","))
DataGroupLabel <- "3groupAnova"
subdata <- data[3:80,1:7] # Select a single trial of data to analyze for 3-way ANOVA
colnames(subdata) <- c("Tx","XCoordinate","Channel","VarianceLabel", "VarianceCode", "BrainRegion", "CellCount") # give the data simple labels
subdata$CellCount <- as.numeric(subdata$CellCount) # convert CellCount to numeric
# Calculate intraclass correlation (ICC)
ANOVAresult <- aov(CellCount ~ XCoordinate, data = subdata)
ANOVAtable <- summary(ANOVAresult)
names(ANOVAtable) <- "Ftable"
Atable <- as.data.frame(ANOVAtable$Ftable)
icc <- Atable[1,3]/(Atable[1,3] + Atable[2,3])
icc
# dummy code the 3 groups (sham = baseline = 0 in $CCI and 0 in $Xpro)
subdata$CCI <- 0
subdata$CCI[subdata$Tx=="CCI"]<- 1
subdata$Xpro <- 0
subdata$Xpro[subdata$Tx=="Xpro"]<- 1
subdata$CCI=as.factor(subdata$CCI)
subdata$Xpro=as.factor(subdata$Xpro)
#subdata$XCoordinate <- type_convert(subdata$XCoordinate)
head(subdata, 100) # Preview & inspect
if (!require('nlme')) install.packages('nlme') # fit multilevel model to run t-test
HCSresult <- nlme::gls(CellCount ~ CCI + Xpro,
# heterogeneous variances
weights=varIdent(form = ~1|VarianceCode),
correlation = corSymm(value=numeric(0),form = ~1|XCoordinate,fixed=T),
data = subdata)
residual <- HCSresult$residuals #setup plots for inspection
par(mfrow=c(1,2))
h <- hist(residual, xlab="Residual")
xfit <- seq(min(residual), max(residual), length = 40)
yfit <- dnorm(xfit, mean = mean(residual), sd = sd(residual))
yfit <- yfit * diff(h$mids[1:2]) * length(residual)
lines(xfit, yfit, col = "black", lwd = 2)
qqnorm(residual, frame = FALSE)
qqline(residual)
par(mfrow=c(1,1))
summary(HCSresult)
# fit null model for comparison
NULresult <- nlme::gls(CellCount ~ 1,
# heterogeneous variances
weights=varIdent(form = ~1|VarianceCode),
correlation = corSymm(value=numeric(0),form = ~1|XCoordinate,fixed=T),
data = subdata)
logLik(NULresult, REML=FALSE)
logLik(HCSresult, REML=FALSE)
chidiff <- -2*as.numeric(logLik(NULresult, REML=FALSE) - logLik(HCSresult, REML=FALSE))
chidiff # print delta chi-square value
qchisq(.95, df=2) # print critical value
round(pchisq(chidiff, df=2, lower.tail=FALSE), 4) # p-value
if (!require('lmeInfo')) install.packages('lmeInfo') # asymmetric conf interval for ES
ctrl <- lmeControl(opt='optim');
result <- lme(fixed = CellCount ~ CCI + Xpro, control=ctrl,
random = ~ 1 | XCoordinate,
correlation = corSymm(value=numeric(0),form = ~1|XCoordinate,fixed=T),
# heterogeneous variances
weights=varIdent(form = ~1|VarianceCode),
data = subdata)
#control=lmeControl(opt="optim") # Switching optimizers for occasional convergence errors
summary(result)
# extract_varcomp(result,separate_variances = TRUE)
ESresult <- lmeInfo::g_mlm(result, p_const = c(0,1,1), r_const = c(1,1,0,0,1,1,1),
infotype = "expected", separate_variances = TRUE)
summary(ESresult)
CI_g(ESresult, symmetric = FALSE)
I fixed a bunch of errors in this to get it to run all the way to the final ESresult <- lmeInfo ::gmlm function call. Unfortunately the lmeInfo module's GitHub page ( https://github.com/jepusto/lmeInfo ) doesn't seem to go into detail on how to specify the r_const = c() vector. My advisor said to use these modules to run this analysis, but is functionally unavailable to teach me how to use them. Specifically, how to understand the string of integers in r_const(). I think I'll need to adjust r_const for other blocks of data in this set, as many of the blocks are of unique sample sizes.
Got this error: Error in utils::combn(1:cor_q, 2) : n < m
...against the following example block of data:
CCI,2108,FITC,Ipsi Group 1,A,Cortex,12.2222222
CCI,1218,FITC,Ipsi Group 1,A,Cortex,19.4736842
CCI,2380,FITC,Ipsi Group 1,A,Cortex,15
CCI,4479,FITC,Ipsi Group 1,A,Cortex,9.47368421
CCI,949,FITC,Ipsi Group 1,A,Cortex,25.2631579
CCI,2129,FITC,Ipsi Group 1,A,Cortex,0
CCI,618,FITC,Ipsi Group 1,A,Cortex,15.5555556
CCI,5061,FITC,Ipsi Group 1,A,Cortex,28.3333333
CCI,2768,FITC,Ipsi Group 1,A,Cortex,25
CCI,2733,FITC,Ipsi Group 1,A,Cortex,31
CCI,2063,FITC,Ipsi Group 1,A,Cortex,25.5
CCI,5307,FITC,Ipsi Group 1,A,Cortex,36
CCI,1923,FITC,Ipsi Group 1,A,Cortex,16.3157895
CCI,9414,FITC,Contra Group 1,D,Cortex,1.05263158
CCI,9123,FITC,Contra Group 1,D,Cortex,4
CCI,10723,FITC,Contra Group 1,D,Cortex,5.2631579
CCI,9737,FITC,Contra Group 1,D,Cortex,6.5
CCI,7422,FITC,Contra Group 1,D,Cortex,8.42105263
CCI,9012,FITC,Contra Group 1,D,Cortex,5
CCI,8924,FITC,Contra Group 1,D,Cortex,5
CCI,9550,FITC,Contra Group 1,D,Cortex,9
CCI,8070,FITC,Contra Group 1,D,Cortex,6
CCI,8734,FITC,Contra Group 1,D,Cortex,7.5
CCI,10154,FITC,Contra Group 1,D,Cortex,9.5
CCI,10203,FITC,Contra Group 1,D,Cortex,38.5
CCI,11305,FITC,Contra Group 1,D,Cortex,24.4444444
Sham,803,FITC,Ipsi Group 2,B,Cortex,1.49342891
Sham,1937,FITC,Ipsi Group 2,B,Cortex,1.88101926
Sham,278,FITC,Ipsi Group 2,B,Cortex,1.38966092
Sham,1414,FITC,Ipsi Group 2,B,Cortex,0
Sham,2762,FITC,Ipsi Group 2,B,Cortex,2.86418056
Sham,3218,FITC,Ipsi Group 2,B,Cortex,1.06126694
Sham,157,FITC,Ipsi Group 2,B,Cortex,5.68860572
Sham,1174,FITC,Ipsi Group 2,B,Cortex,3.09463391
Sham,4481,FITC,Ipsi Group 2,B,Cortex,5.33253701
Sham,4011,FITC,Ipsi Group 2,B,Cortex,4.53607312
Sham,2784,FITC,Ipsi Group 2,B,Cortex,1.1154863
Sham,2379,FITC,Ipsi Group 2,B,Cortex,1.05037604
Sham,1417,FITC,Ipsi Group 2,B,Cortex,0
Sham,6739,FITC,Contra Group 2,E,Cortex,9.59452079
Sham,5166,FITC,Contra Group 2,E,Cortex,2.18957325
Sham,7252,FITC,Contra Group 2,E,Cortex,5.16175655
Sham,7115,FITC,Contra Group 2,E,Cortex,9.05843988
Sham,8000,FITC,Contra Group 2,E,Cortex,11.5952485
Sham,5843,FITC,Contra Group 2,E,Cortex,0
Sham,6361,FITC,Contra Group 2,E,Cortex,4.84261501
Sham,8444,FITC,Contra Group 2,E,Cortex,6.10209357
Sham,8247,FITC,Contra Group 2,E,Cortex,5.22677678
Sham,5260,FITC,Contra Group 2,E,Cortex,3.88888889
Sham,7526,FITC,Contra Group 2,E,Cortex,3.97645936
Sham,7495,FITC,Contra Group 2,E,Cortex,5.7390456
Sham,8052,FITC,Contra Group 2,E,Cortex,3.17302506
Xpro,2845,FITC,Ipsi Group 3,C,Cortex,33.6842105
Xpro,1379,FITC,Ipsi Group 3,C,Cortex,9.09090909
Xpro,1592,FITC,Ipsi Group 3,C,Cortex,21.6666667
Xpro,968,FITC,Ipsi Group 3,C,Cortex,6.25
Xpro,2982,FITC,Ipsi Group 3,C,Cortex,17.2222222
Xpro,1212,FITC,Ipsi Group 3,C,Cortex,6.31578947
Xpro,2563,FITC,Ipsi Group 3,C,Cortex,29
Xpro,4615,FITC,Ipsi Group 3,C,Cortex,30
Xpro,630,FITC,Ipsi Group 3,C,Cortex,23
Xpro,5014,FITC,Ipsi Group 3,C,Cortex,27
Xpro,1409,FITC,Ipsi Group 3,C,Cortex,16.5
Xpro,3969,FITC,Ipsi Group 3,C,Cortex,23
Xpro,3527,FITC,Ipsi Group 3,C,Cortex,15
Xpro,7253,FITC,Contra Group 3,F,Cortex,7
Xpro,7891,FITC,Contra Group 3,F,Cortex,9
Xpro,6711,FITC,Contra Group 3,F,Cortex,15.625
Xpro,7830,FITC,Contra Group 3,F,Cortex,11.1111111
Xpro,7218,FITC,Contra Group 3,F,Cortex,10.2857143
Xpro,7849,FITC,Contra Group 3,F,Cortex,4.57085395
Xpro,7303,FITC,Contra Group 3,F,Cortex,12.6019498
Xpro,6928,FITC,Contra Group 3,F,Cortex,33.8199036
Xpro,7991,FITC,Contra Group 3,F,Cortex,5.81483248
Xpro,8452,FITC,Contra Group 3,F,Cortex,1.28539661
Xpro,6717,FITC,Contra Group 3,F,Cortex,36.0048967
Xpro,7328,FITC,Contra Group 3,F,Cortex,23.0207289
Xpro,5772,FITC,Contra Group 3,F,Cortex,7
So, I've troubleshot all the way to the end and I'm stuck here on the last lines of code. Much appreciation to anyone for wisdom.
