#LMM
library(lme4)
m2 <- clm(categorical_FOD_FODs ~ 0 + Condition_FODs:Language_used_FODs - 1, random = ~1|subject, data = indvar_FODs)
summary(m2)
#LM
m2 <- clm(categorical_FOD_FODs ~ Condition_FODs:Language_used_FODs - 1, data = indvar_FODs)
summary(m2)
#post-hoc analysis
library(lsmeans)
contrast(m2, "pairwise", list(Condition_FODs = c("A", "B", "C", "D")))
library(emmeans)
emmeans(m2, pairwise~Condition_FODs, adjust = "bonferroni")
# dependent variable
y <- as.factor(indvar_FODs$categorical_FOD_FODs)
# independent variables
x1 <- indvar_FODs$Condition_FODs
x2 <- indvar_FODs$Language_used_FODs
# model with optimizer argument
mod_FODs <- glmer(y ~ x1 * x2 + (1|subject_FODs), family = binomial, data = indvar_FODs,
control = glmerControl(optimizer = "bobyqa"))
# print model
summary(mod_FODs)
fix invalid codeMon, 01 May 2023 Model_DD_1 <- glmer(cat_SNB_new ~ Condition_SNB + Gender_SNB +
(1+Condition_SNB|subject_SNB) + (1+Condition_SNB|Gender_SNB),
data = s_complete,
family = binomial)
# model fitting
model_FODs <- multinom(categorical_FOD_FODs ~ ., data = indvar_FODs)
summary(model_FODs)
library(multcomp)
anova_result <- aov(Trait ~ Condition, data = data)
posthoc <- glht(anova_result, linfct = mcp(Condition_SNB = "Tukey"))
summary(posthoc)
# this is the code for a latex table
library(texreg)
screenreg(anova_result, posthoc)
library(lme4)
library(optimx)
# dependent variable
y <- as.factor(indvar_FODs$categorical_FOD_FODs)
# independent variables
x1 <- indvar_FODs$Condition_FODs
x2 <- indvar_FODs$Language_used_FODs
# model
mod_FODs <- glmer(y ~ x1*x2 + (1|subject_FODs), family=binomial, data=indvar_FODs)
# print model
mod_FODs
fix invalid codeSat, 06 May 2023
fix invalid codeFri, 05 May 2023 r
s_complete$Handedness_SNB <- factor(s_complete$Handedness_SNB)
table(s_complete$Handedness_SNB)
levels(s_complete$Handedness_SNB)
# R script
library(nlme)
library(lme4)
# analysis on the effect of Cond on direction
lmer(direction ~ Cond + (1 | colours), data = indvar_FODs_new)
# analysis on the interaction effect of Cond and colours on direction
lmer(direction ~ Cond * colours + (1 | colours), data = indvar_FODs_new)
I want to run a likelihood ratio test for cumulative link models. For that, I need to specify two models (m1 and m2).
library(multcomp)
summary(glht(model, linfct = mcp (Condition_FODs*Language_used_FODs="Tukey")))
fix invalid codeFri, 05 May 2023 > s_complete$developmentaldisorder_SNB <- factor(s_complete$developmentaldisorder_SNB)
> table(s_complete$developmentaldisorder_SNB)
No Yes
2712 796
> levels(s_complete$developmentaldisorder_SNB)
NULL
fix invalid codeFri, 12 May 2023 ### Getting and Cleaning Data
#### Installation
Download and install the latest version of [R](https://cran.r-project.org/), and [RStudio](https://www.rstudio.com/products/rstudio/download/).
Clone the [GitHub repo](https://github.com/bkkaggle/getting-cleaning-data) or download the zip file and extract.
#### Usage
Open the `run_analysis.R` file in RStudio.
Set the working directory to the extracted directory.
Source the script
indvar_FODs_new
Cond colours direction
1 Ay blue morning
2 Ay pink morning
3 Ay pink morning
4 Ay pink morning
5 Ay pink morning
6 Ay pink morning
7 Ay pink morning
8 Ay blue morning
9 Ay blue morning
10 Ay pink morning
11 Ay pink morning
12 Ay blue morning
13 Ay blue morning
14 Ay blue morning
15 Ay blue morning
16 Ay blue morning
17 Ay blue morning
18 Ay blue morning
19 Ay pink morning
20 Ay blue morning
21 Ay blue morning
22 Ay blue morning
23 Ay blue morning
24 Ay blue morning
25 Ay blue morning
26 Ay blue morning
27 Ay blue morning
28 Ay blue morning
29 Ay blue morning
30
# Group by Gender_SNB and num_SNB, and calculate summary statistics
av_data1_all <- s_complete %>%
group_by(Gender_SNB, num_SNB) %>%
summarise(mean_num_SNB = mean(num_SNB),
sum_synch = sum(synch),
count = n()) %>%
ungroup()
# Your code here...
### Problem 3.2
fix invalid codeFri, 12 May 2023 Fixing invalid data:
mod0 <- glmer(categorical_FOD_FODs~Condition_FODs*Gender_FODs+
(1|subject_FODs),
na.action=na.omit,
data=indvar_FODs,
control = lmeControl(opt = "optim"))
## create a design matrix
FODs_design <- model.matrix(~Language_used_FODs*Condition_FODs , data=indvar_FODs)
FODs_design <- FODs_design[,-1]
head(FODs_design)
attach(indvar_FODs)
FODs_model <- multinom(categorical_FOD_FODs ~ Language_used_FODs*Condition_FODs , data=indvar_FODs, Hess=TRUE)
summary(FODs_model)
summary(FODs_model)$coefficients
aic(FODs_model)
anova(FODs_model, type="Chisq")
Analysis of Deviance Table
Model: multinom(formula = categorical_FOD_FODs ~ Language_used_FODs * Condition_FODs,
data = indvar_FODs, Hess = TRUE)
Terms added sequentially (
fix invalid codeMon, 24 Apr 2023
fix invalid codeFri, 12 May 2023
fix invalid codeFri, 12 May 2023
> s_complete$cat_SNB_new[s_complete$cat_SNB_new == "0"] <- "No"
> s_complete$cat_SNB_new[s_complete$cat_SNB_new == "1"] <- "Yes"
> model <- lmer(cat_SNB_new ~ Condition_SNB * Dotcounting_SNB + (1 | subject_SNB), data = na.omit(s_complete))
> summary(model)
Linear mixed model fit by REML ['lmerMod']
Formula: cat_SNB_new ~ Condition_SNB * Dotcounting_SNB + (1 | subject_SNB)
Data: na.omit(s_complete)
REML criterion at convergence: 611.2
Scaled residuals:
Min 1Q Median 3Q Max
-3.2218 -0.8447 0.0724 0.8313 3.1351
Random effects:
Groups Name Variance Std
glm(categorical_FOD_FODs ~ 0 + Condition_FODs:Language_used_FODs, data = indvar_FODs, family=binomial(link="logit"))
get_colours <- function(ggplot2_object){
ggplot2_object +
theme(panel.background = element_rect(fill = NA)) +
theme(panel.grid.major = element_blank()) +
theme(panel.grid.minor = element_blank()) +
theme(panel.border = element_rect(fill = NA)) +
theme(axis.line = element_blank()) +
theme(axis.text.x = element_blank()) +
theme(axis.text.y = element_blank()) +
theme(axis.ticks = element_blank()) +
theme(axis.title.x = element_blank()) +
theme(axis.title.y = element_blank()) +
theme(plot.background = element_blank()) +
theme(plot.title = element_blank()) +
theme(legend.position = "none") +
theme(legend.background = element_blank()) +
theme(legend.key = element_blank()) +
theme(legend
mod0 <- glmer(categorical_FOD_FODs~Condition_FODs*Gender_FODs+
(1|subject_FODs),
na.action=na.omit,
data=indvar_FODs,
control = lmeControl(opt = "optim"))
The ideal test to conduct is the likelihood ratio test.
The group means of the FOD_FODs are not significantly different from zero (M = -0.0644, SE = 0.0058, t = -11.045, p < .001) and the variance is also significantly different from zero (M = 1.10, SE = 0.009, t = 120.818, p < .001).
The model is not significant (M = 0.4038, SE = 0.0406, t = 9.983, p < .001) and the condition has no effect (M = 0.3213, SE = 0.3236, t = 1.000, p = 0.324).
The model is not significant (M = 0.0147, SE = 0.0074, t = 1.997, p = 0.048) and the condition has no effect (M = 0.98738, SE = 0.0155, t = 63.904, p < .001).
fix invalid codeFri, 12 May 2023
fix invalid codeWed, 24 May 2023 ### `tidy.anova` (with `print.tidy.anova`)
Tidy the result of `anova` into a table.
fix invalid codeSat, 06 May 2023 > anova_table <- Anova(model, type = "III")
> library(report)
> report(anova_table)
Error in `[[<-.data.frame`(`*tmp*`, "Mean_Square", value = numeric(0)) :
Ersetzung hat 0 Zeilen, Daten haben 3
T <- chisq.test(cont_table)
T
fix invalid codeMon, 24 Apr 2023 ## Regression
Regression analysis is a form of predictive modelling technique which investigates the relationship between a dependent (target) and independent variable (s) (predictor). This technique is used for forecasting, time series modelling and finding the causal effect relationship between the variables. For example, relationship between rash driving and number of road accidents by a driver is best studied through regression.
Regression analysis is an important tool for modelling and analyzing data. Here, we fit a regression line to observed data. One sample of observations is taken and a line is fitted to the data. The best fit line is known as regression line and represented by a linear equation Y=a+bX.
Linear regression is a statistical approach for modelling relationship between a dependent variable with a given set of independent variables.
Mathematically a linear relationship represents a straight line when plotted as a graph. And a non-linear relationship is represented by a curve in a graph.
The term “linearity” in algebra refers to a linear relationship between two or more variables. If we draw this relationship in a two-dimensional space (between two variables, in this case), we get a straight line.
library(ordinal)
m2 <- clm(categorical_FOD_FODs ~ 0 + Condition_FODs:Gender_quantised_FODs - 1,
random = ~1|subject_FODs, data = indvar_FODs)
summary(m2)
fix invalid codeFri, 12 May 2023
# Add your function here
library(nlme)
model1 <- glmer(categorical_FOD_FODs~Condition_FODs*Gender_FODs,random=~1|subject_FODs,
na.action=na.omit,
method = "ML",
data=indvar_FODs,
control = lmeControl(opt = "optim"))
fix invalid codeWed, 24 May 2023 ### Answer:
TukeyHSD(aov(L1 ~ Condition*Language, data = indvar))
## R code for the 2nd and 3rd question:
# model with optimizer argument
mod_FODs <- glmer(y ~ x1 * x2 + (1|subject_FODs), family = binomial, data = indvar_FODs,
control = glmerControl(optimizer = "bobyqa"))
# print model
summary(mod_FODs)
fix invalid codeWed, 24 May 2023
# Write your code here!
library(car)
TukeyHSD(aov(categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + Condition_FODs*Language_used_FODs, indvar_FODs))
pvtable(m2)
The results are reported in APA style in the table below.
The last example is a good one. The formula is categorical_FOD_FODs ~ 0 + Condition_FODs:Gender_quantised_FODs - 1
and this is the formula for the model. The m2 is the model with the formula. However, What does m1 stand for? What does 0 stand for? what does the : stand for? What does the - 1 stand for?
fix invalid codeMon, 01 May 2023 ## Fix order
def add(a, b):
return a + b
add(1, 2)
The stepwise AIC algorithm selects the model with the lowest AIC value. According to this test we do not reject the null hypothesis, and therefore, we do not select the model m2.
fix invalid codeFri, 12 May 2023
> my_data <- indvar_FODs %>%
+ group_by(Condition_FODs, Language_used_FODs) %>%
+ summarise(mean = mean(categorical_FOD_FODs))
> my_data
# A tibble: 20 x 3
# Groups: Condition_FODs [4]
Condition_FODs Language_used_FODs mean
<fct> <fct>
1 A English 0.7477272
2 A German 0.6437500
3 A Hungarian 0.8125000
4 A Italian 0.5375000
5 A Turkish 0.8689873
6 B English 0.7738636
7 B German 0.8875000
8 B Hungarian 0.9062500
9 B Italian 0.7125000
10 B Turkish 0.9528302
11 C English 0.7397727
12 C German 0.618
chisq.test(
cont_table,
sim = TRUE,
B = 10000,
correct = FALSE
)
simulate this in python
fix invalid codeMon, 01 May 2023 > Model_DD <- glmer(cat_SNB_new ~ Condition_SNB * Gender_SNB + (1 + Condition_SNB | subject_SNB),
+ data = s_complete,
+ family = binomial,
+ control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 20000)))
boundary (singular) fit: see help('isSingular')
avfourfods <- ggplot(meanFODs_indvar_FODs, aes(x = Condition_FODs, y = mean, fill = Condition_FODs)) +
geom_bar(stat = "identity", position = "dodge") +
geom_errorbar(aes(ymin = mean - sd/sqrt(count), ymax = mean + sd/sqrt(count)),
width = 0.2, position = position_dodge(width = 0.9)) +
xlab("Condition") +
ylab("Mean of FODs") +
ggtitle("Average of first order differences (FODs): Four conditions") +
theme_minimal() +
theme(panel.border = element_rect(color = "black", fill = NA, size = 1)) +
labs(fill = "Conditions") +
geom_hline(yintercept = 0, color = "black", linetype = "solid", size = 2)
#model1
model1<-aov(categorical_FOD_FODs ~ Condition_FODs + Gender_FODs, data=indvar_FODs)
summary(model1)
#model2
model2<-aov(categorical_FOD_FODs ~ Condition_FODs*Gender_FODs, data=indvar_FODs)
summary(model2)
According to the Shapiro-Wilk test, the residual is not normally distributed (p=0). However, the model is not singular (FALSE) and the effect is not significant (p=1).
fix invalid codeWed, 24 May 2023 ## Fixing Errors
### Fixing Errors
- You can use `\` to wrap code across lines
- Use `{}` to wrap code across multiple lines
- Use `#` to comment out lines
fix invalid codeWed, 24 May 2023 ## Recommandations
- Please use an R package to interface with markdown.
- Please use an R package to interface with HTML.
- Please use an R package to interface with LaTeX.
- Please use an R package to interface with `.docx` files.
- Please use an R package to interface with `.xlsx` files.
- Please use an R package to interface with `.pptx` files.
- Please use an R package to interface with `.pdf` files.
- Please use an R package to interface with `.mp4` files.
- Please use an R package to interface with `.mp3` files.
- Please use an R package to interface with `.png` files.
- Please use an R package to interface with `.jpg` files.
- Please use an R package to interface with `.gif` files.
- Please use an R package to interface with .csv files.
- Please use an R package to interface with `.xls` files.
- Please use an R package to interface with `.ppt` files.
- Please use an
These are code and results. My question is how to report these results in APA style.
# fit the model
fit <- glm(categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + Condition_FODs:Language_used_FODs, data = indvar_FODs, family = "binomial")
# optimizer
# modify the optimizer
fit2 <- glm(categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + Condition_FODs:Language_used_FODs, data = indvar_FODs, family = "binomial",
control = glm.control(epsilon = 1e-4, maxit = 100))
# if singularity
# the model is singular, therefore the model is not valid
fit3 <- glm(categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + Condition_FODs:Language_used_FODs, data = indvar_FODs, family = "binomial",
control = gl
fix invalid codeMon, 01 May 2023 # How to contribute
Please read [CONTRIBUTING.md](https://gist.github.com/PurpleBooth/b24679402957c63ec426) for details on our code of conduct, and the process for submitting pull requests to us.
# Versioning
We use [SemVer](http://semver.org/) for versioning. For the versions available, see the [tags on this repository](https://github.com/your/project/tags).
# Authors
- **Billie Thompson** - _Initial work_ - [PurpleBooth](https://github.com/PurpleBooth)
See also the list of [contributors](https://github.com/your/project/contributors) who participated in this project.
# License
This project is licensed under the MIT License - see the [LICENSE.md](LICENSE.md) file for details
# Acknowledgments
- Hat tip to anyone whose code was used
- Inspiration
- etc
res.aov1 <- aov(categorical_FOD_FODs ~ Language_used_FODs + Error(Condition_FODs / Language_used_FODs), indvar_FODs)
summary(res.aov1)
library(MASS)
library(lsmeans)
model <- multinom(categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + Condition_FODs*Language_used_FODs, indvar_FODs)
summary(model)
post_hoc <- lsmeans(model, ~ Condition_FODs + Language_used_FODs + Condition_FODs*Language_used_FODs, adjust = "tukey")
post_hoc
fix invalid codeFri, 12 May 2023 ## Summary
[table](https://www.tablesgenerator.com/markdown_tables)
| Study | Frequencies |
| ----------------- | ----------- |
| Dotcounting_SNB | A-B-C-D |
| Dotcounting_SNB | A-C-B-D |
| Dotcounting_SNB | A-D-B-C |
| Dotcounting_SNB | B-A-C-D |
| Dotcounting_SNB | B-A-D-C |
| Dotcounting_SNB | B-C-D-A |
| Dotcounting_SNB | D-C-B-A |
| Dotcounting_SNB | missing |
| Flanker_SNB | Congruent |
| Flanker_SNB | Incongruent|
| Flanker_SNB | missing |
| Stroop_
# define a function
myfun <- function(x){
a <- mean(x)
b <- sd(x)
c <- table(x)
return(c(a, b, c))
}
# use apply over the dataframe
apply(s_complete[c("num_SNB", "Condition_SNB")], 2, myfun)
model1 <- glmer(categorical_FOD_FODs ~ Condition_FODs + (1|subject_FODs),
data = indvar_FODs, family = "binomial", control = glmerControl(optimizer = "bobyqa"))
fix invalid codeWed, 10 May 2023 **Question**: How do you interpret the results of the model (i.e. beta weights, p-values)?
fix invalid codeMon, 24 Apr 2023 What is the problem here?
## 8.0 References
1. https://en.wikipedia.org/wiki/Statistical_inference
2. https://en.wikipedia.org/wiki/Statistical_hypothesis_testing
3. http://www.stat.yale.edu/Courses/1997-98/101/meancomp.htm
4. https://cran.r-project.org/web/packages/lme4/vignettes/lmer.pdf
5. https://cran.r-project.org/web/packages/lmerTest/vignettes/lmerTest.pdf
6. https://neuroconductor.org/packages/release/bioc/html/jmv.html
7. https://neuroconductor.org/packages/release/bioc/html/lme4.html
8. https://neuroconductor.org/packages/release/bioc/html/lmerTest.html
9. https://neuroconductor.org/packages/release/bioc/html/pbkrtest.html
fix invalid codeFri, 05 May 2023 [转载](https://www.jianshu.com/p/7e9eb09e8dc7)
#run model using lme4 package
model <- lmer(FOD_FODs ~ Condition_FODs + (1|subject_FODs), data = indvar_FODs)
fix invalid codeSat, 06 May 2023 model_full2 <- glmer(cat_SNB_new ~ 0 + developmentaldisorder_SNB * Condition_SNB -1 + (1+developmentaldisorder_SNB|subject_SNB),
data = s_complete, family = binomial)
levels(s_complete$Gender_SNB)
[1] "Female" "Male" "Non-Binary / Genderqueer"
anova(Model_DD_1)
Analysis of Variance Table
npar Sum Sq Mean Sq F value
Condition_SNB 4 136.301 34.075 34.0754
Gender_SNB 2 1.684 0.842 0.8421
Condition_SNB:Gender_SNB 6 7.063 1.177 1.1771
fix invalid codeMon, 24 Apr 2023 ---
### indexing in R
- `(x)`: `x` is an object, in R, it is a vectore or data frame
- `x[i]`: `x` is a vector, `i` is an integer, which is the index of `x`
- `x[[i]]`: `x` is a data frame, `i` is an integer, which is the sequence number of `x`
- `x$i`: `x` is a data frame, `i` is a string, which is the variable name
- `x[i,j]`: `x` is a matrix, `i` and `j` are integers, which are the index of `x`
- `x[[i,j]]`: `x` is a data frame, `i` and `j` are integers, which are the sequence number of `x`
- `x$i[j]`: `x` is a data frame, `i` is a string, and `j` is an integer, which is the variable name and index
- `
fix invalid codeMon, 01 May 2023
fix invalid codeWed, 24 May 2023 > eta_squared <- eta_squared(model, anova_table = tidy(model))
> eta_squared
analytic numeric
0.3075934 0.3075899
fix invalid codeWed, 10 May 2023 ***NOTE:***
- This is just a draft
- We are not sure if this will work
The function add takes two arguments a and b and returns the sum of a and b.
The model1 variable is the output of applying the linear mixed-effect model analysis to the data in the variable long1.RSPL using the linear mixed-effect model analysis package lme4. The linear mixed-effect model analysis takes as input the formula categorical_FOD_FODs~Condition_FODs*Gender_FODs, (the dependent variable categorical_FOD_FODs is predicted by the predictor variables Condition_FODs and Gender_FODs and their interaction Condition_FODs*Gender_FODs), random=~1|ID/Ch, (the random effects are the intercepts nested in the variables ID and Ch), na.action=na.omit (missing data points are removed), method = "ML" (maximum likelihood estimation is used), data=long1.RSPL (the data is the data in the variable long1.RSPL), control = lmeControl(opt = "optim") (the control is the optimization process).
write the APA style report here
fix invalid codeMon, 24 Apr 2023 ggplot(melt(av_data1_all, id.vars=c("ID","Gender_SNB", "Age_SNB")), aes(x=num_SNB, y=value, color=Age_SNB, group=Age_SNB))+
geom_smooth(method=glm, method.args= list(family = binomial(logit)),
se = FALSE) +
geom_point() +
#geom_point(aes(shape = block)) +
xlab('Delay (ms)') + ylab('Prop. of Synch Responses') +
#ggtitle("Session 2")+
theme_bw()+ # removes background
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+ #removes grid
facet_wrap(~ID)+
scale_colour_manual(values = c("red", "blue", "black"))+
theme(axis.text.x = element_text(angle =
library(lme4)
library(ggplot2)
library(MASS)
s_complete = read.csv('s_complete.csv')
s_complete$Condition_SNB = factor(s_complete$Condition_SNB)
s_complete$Dotcounting_SNB = factor(s_complete$Dotcounting_SNB)
s_complete$cat_SNB_new = factor(s_complete$cat_SNB_new)
#create dataframe with interactions
s_complete$Dotcounting_SNB_Condition_SNB = interaction(s_complete$Dotcounting_SNB, s_complete$Condition_SNB)
Model_MSG <- glmer(cat_SNB_new ~ 0 + Condition_SNB * Dotcounting_SNB - 1 + (1 | subject_SNB),
data = s_complete,
family = binomial,
control = glmerControl(optimizer = "bobyqa", tolPwrss = 1e-08))
library(dplyr)
s_complete %>%
group_by(Condition_SNB) %>%
summarise(mean = mean(num_SNB), sd = sd(num_SNB), count = length(num_SNB), "max_value" = max(num_SNB), "min_value" = min(num_SNB))
library(nlme)
model <- lme(FOD_FODs ~ Condition_FODs, random = ~1|subject_FODs, data = indvar_FODs)
fix invalid codeSat, 22 Apr 2023 model <- lme(FOD_FODs ~ Condition_FODs, random = ~1|subject_FODs, data = indvar_FODs)
anova(model)
glmer('categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + (1 | subject_FODs)',
+ data = indvar_FODs, family = binomial, method = "", data= ,control = (opt = ""))
glm(categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + Condition_FODs*Language_used_FODs,
family = multinomial("identity"), data = indvar_FODs)
library(multcomp)
glht(model, linfct = mcp(Condition_FODs = "Tukey"))
# load lme4 package
library(lme4)
# load multicomp package
library(multicomp)
# load the data
IndVar_FODs <- read.csv("/Users/sarahszalay/Dropbox/Thesis/Data/IndVar_FODs.csv", header = TRUE)
IndVar_FODs$categorical_FOD_FODs <- factor(IndVar_FODs$categorical_FOD_FODs)
IndVar_FODs$Condition_FODs <- factor(IndVar_FODs$Condition_FODs)
IndVar_FODs$Language_used_FODs <- factor(IndVar_FODs$Language_used_FODs)
IndVar_FODs$subject_FODs <- factor(IndVar_FODs$subject_FODs)
IndVar_FODs$Language_used_FODs <- relevel(IndVar_FODs$Language_used_FODs, "English")
IndVar_F
fix invalid codeTue, 18 Apr 2023 ## Installing
You can download the package in R with:
fix invalid codeWed, 24 May 2023 r
# 3.3
lm.fit = lm(Word.Accuracy~Condition*Language_used, data = data)
anova(lm.fit)
m1 <- glht(model, linfct = mcp(Condition_FODs = "Tukey"))
summary(m1)
indvar_FODs <- ifelse(data$nativelang_FODs == "Danish", "Danish",
ifelse(data$nativelang_FODs == "English", "English",
ifelse(data$nativelang_FODs == "German", "German",
ifelse(data$nativelang_FODs == "Hungarian", "Hungarian",
ifelse(data$nativelang_FODs == "Italian", "Italian",
ifelse(data$nativelang_FODs == "Polish", "Polish",
ifelse(data$nativelang_FODs == "Russian", "Russian",
ifelse(data$nativelang_FODs == "Turkish", "Turkish",
ifelse(data$nativelang_FODs == "Ukrainian-Russian", "Ukrainian-Russian",
ifelse(data$nativelang_FODs == "Turkish-German", "Turkish-German",
> indvar_FODs$categorical_FOD_FODs
ascending descending identity
1683 1432 206
> indvar_FODs$categorical_FOD_FODs
ascending descending identity
1683 1432 206
> indvar_FODs$categorical_FOD_FODs
ascending descending identity
1683 1432 206
> indvar_FODs$categorical_FOD_FODs
ascending descending identity
1683 1432 206
> indvar_FODs$categorical_FOD_FODs
ascending descending identity
1683 1432 206
> indvar_FODs$categorical_FOD_FODs
ascending descending identity
1683 1432 206
> indvar_FODs$categorical_FOD_FODs
ascending descending identity
fix invalid codeFri, 12 May 2023
# dependent variable
y <- as.factor(indvar_FODs$categorical_FOD_FODs)
# independent variables
x1 <- indvar_FODs$Condition_FODs
x2 <- indvar_FODs$Language_used_FODs
# model with optimizer argument
mod_FODs <- glmer(y ~ x1 * x2 + (1|subject_FODs), family = binomial, data = indvar_FODs,
control = glmerControl(optimizer = "bobyqa"))
# print model
summary(mod_FODs)
glm_indvar_FODs<- glm(categorical_FOD_FODs ~ ., data=indvar_FODs, family = "multinomial")
summary(glm_indvar_FODs)
library(knitr)
library(kableExtra)
library(tidyverse)
data
%>%
select(.matrix)
%>%
filter(p.value < 0.05)
%>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
# This is a function that will take four inputs:
# 1) matrix: a matrix of data
# 2) outcome: the column name of the outcome/dependent variable
# 3) predictors: a vector of column names, indicating which columns are the predictors/independent variables
# 4) family: the family object, e.g. binomial()
# This function will return a summary of a logistic regression model.
# The function requires that you have the 'dplyr' package installed.
log_reg <- function(matrix, outcome, predictors, family) {
# Install and load the 'dplyr' package.
#install.packages("dplyr")
#library(dplyr)
# Fit a logistic regression model, using the variables in 'predictors'
# as the predictors and the variable in 'outcome' as the outcome
# The model should use the family indicated by 'family'.
# Make sure to save your model in the 'fit' object!
#fit <- glm
library(lme4)
# dependent variable
y <- as.factor(indvar_FODs$categorical_FOD_FODs)
# independent variables
x1 <- indvar_FODs$Condition_FODs
x2 <- indvar_FODs$Language_used_FODs
# model
mod_FODs <- glmer(y ~ x1*x2 + (1|subject_FODs), family=binomial, data=indvar_FODs)
# print model
mod_FODs
> multinom(categorical_FOD_FODs~Condition_FODs*Language_used_FODs, data = indvar_FODs)
The function add takes two arguments a and b and returns the sum of a and b.
print(pairwise.t.test(dataset$HDL.C, dataset$DIABETES, p.adjust.method = "bonferroni"))
condition_FODs: F(4, 3139)=0.35, p=0.846, η2=4.42e-4
Age_FODs: F(15, 3139)=0.084, p=1.000, η2=4.02e-4
Condition_FODs:Age_FODs: F(45, 3139)=0.046, p=1.000, η2=6.52e-4
library(knitr)
library(kableExtra)
crosstab <- kable(
mtcars[, c(1,3,5,6)], "html",
caption = "Cross-tabulation of gears and carb",
booktabs = TRUE
) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
crosstab
fix invalid codeWed, 24 May 2023 Or you can provide a url:
The BIC test is ideal to conduct.
We used linear mixed effect model to analyse the interaction effect of condition and gender on the categorical FODs data.
posthoc <- glht(anova_result, linfct = mcp(Condition_SNB = "Tukey"))
# create a table of the results
tab <- as.data.frame(summary(posthoc))
tab <- tab[-1,]
#rename the columns
names(tab)[names(tab) == "p values"] <- "p"
names(tab)[names(tab) == "Estimate"] <- "Est"
names(tab)[names(tab) == "Estimate[Standard Error]"] <- "Std.Err"
names(tab)[names(tab) == "Std.Err"] <- "Std.Error"
names(tab)[names(tab) == "Lower"] <- "Lower.CI"
names(tab)[names(tab) == "Upper"] <- "Upper.CI"
names(tab)[names(tab) == "Lower.CI"] <- "Lower.95CI"
names(tab)[names(tab) == "Upper.CI"] <- "Upper.95CI"
#create the table
print(
library(nnet)
set.seed(1234)
# split dataset into train and test
train <- sample(1:nrow(indvar_FODs), nrow(indvar_FODs)*.7)
test <- -train
# this function will return a predictions where the sum of probabilities of all
# possible values of the categorical_FOD_FODs is equal to 1.
mn_logit_model <- multinom(categorical_FOD_FODs ~ Condition_FODs*Language_used_FODs,
data = indvar_FODs[train,], maxit = 50)
# prediction
mn_logit_model_pred <- predict(mn_logit_model, indvar_FODs[test,])
# accuracy
mean(max.col(mn_logit_model_pred)==max.col(indvar_FODs[test,]))
fix invalid codeSat, 06 May 2023 ### Python
fix invalid codeWed, 10 May 2023 <details>
<summary>Click to see the solution</summary>
<p>
- The model formula is not correct. There is no random effect in your model.
- The model has not converged. Change the optimizer to NAG, and it runs really good.
library(nnet)
m1 <- multinom(categorical_FOD_FODs ~ Condition_FODs*Language_used_FODs, data = indvar_FODs)
summary(m1)
Call:
multinom(formula = categorical_FOD_FODs ~ Condition_FODs * Language_used_FODs,
data = indvar_FODs)
Coefficients:
(Intercept) Condition_FODsA Condition_FODsB
-5.56764 0.16689 -0.01620
Condition_FODsC Condition_FODsD Language_used_FODsEnglish
0.09732 -0.16689 0.36466
Language_used_FODsGerman Language_used_FODsHungarian Language_used_FODsItalian
-0.36466 -0.36466 -0.03095
Language_used_F
# dependent variable
y <- as.factor(indvar_FODs$categorical_FOD_FODs)
# independent variables
x1 <- indvar_FODs$Condition_FODs
x2 <- indvar_FODs$Language_used_FODs
# model
mod_FODs <- glmer(y ~ x1*x2 + (1|subject_FODs), family=binomial, data=indvar_FODs)
# print model
mod_FODs
fix invalid codeFri, 12 May 2023
def latex_table(table_header, table_data):
table = table_header
...
chisq.test(cont_table)
# multinomial logistic regression
glm(categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs,
family = ,
data = indvar_FODs)
fix invalid codeSat, 06 May 2023
# for categorical_FOD_FODs
MultinomialLogisticRegression_FODs <- multinom(categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + Condition_FODs:Language_used_FODs, data = indvar_FODs)
summary(MultinomialLogisticRegression_FODs)
# Call:
# multinom(formula = categorical_FOD_FODs ~ Condition_FODs +
# Language_used_FODs + Condition_FODs:Language_used_FODs,
# data = indvar_FODs)
#
# Coefficients:
# (Intercept) Condition_FODsB Condition_FODsC
# 1 -0.87860 0.41173
# Condition_FODsD Language_used_FODsGerman
# 1.0 -0.12083
#
model_FODs_1 = glmer(categorical_FOD_FODs ~ Condition_FODs * Language_used_FODs + (1|subject_FODs), data = indvar_FODs, family = 'binomial')
fix invalid codeWed, 24 May 2023 {r, warning=FALSE, message=FALSE, echo=FALSE}
eta_squared(model)
## #2: Plot a function
plot(x, y)
The function `plot` takes two arguments, `x` and `y`, and plots an x-y graph of the two.
## #3: Plot different types of functions
plot(x, y1); plot(x, y2)
The function `plot` takes two arguments, `x` and `y1`, and plots an x-y graph of the two. The function `plot` takes two arguments, `x` and `y2`, and plots an x-y graph of the two.
## #4: Combine multiple lines of code
def add(a, b):
return a + b
add(2, 5)
The function `add` takes two arguments, `a` and `b`, and returns the sum of a and b. The function `add` takes two arguments, `2` and `5`, and returns the sum of `2` and `5`.
## #5: Put a function inside a function
def add(a, b):
return a + b
def model_snb(s_complete):
s_complete <- read.csv("s_complete_snb_with_model.csv")
Model_MSG <- glmer(cat_SNB_new ~ 0 + Condition_SNB * Dotcounting_SNB - 1 + (1 | subject_SNB),
data = s_complete,
family = binomial,
control = glmerControl(optimizer = "bobyqa", tolPwrss = 1e-08))
return(Model_MSG)
model_snb(s_complete)
data(iris)
iris <- within(iris, {Species <- factor(Species)})
model <- lm(Sepal.Length ~ Species * Sepal.Width, data = iris)
TukeyHSD(model, "Species")
TukeyHSD(model, "Sepal.Width")
TukeyHSD(model, "Species:Sepal.Width")
posthoc <- glht(anova_result, linfct = mcp(Condition_SNB = "Tukey"))
posthoc_summ <- summary(posthoc)
posthoc_summ
avfourfods <- ggplot(meanFODs_indvar_FODs, aes(x = Condition_FODs, y = mean, fill = Condition_FODs)) +
geom_bar(stat = "identity", position = "dodge") +
geom_errorbar(aes(ymin = mean - sd/sqrt(count), ymax = mean + sd/sqrt(count)),
width = 0.2, position = position_dodge(width = 0.9)) +
xlab("Condition") +
ylab("Mean of FODs") +
ggtitle("Average of first order differences (FODs): Four conditions") +
theme_minimal() +
theme(panel.border = element_rect(color = "black", fill = NA, size = 1)) +
labs(fill = "Conditions") +
geom_hline(yintercept = 0, color = "black", linetype = "solid", size = 2)
#install.packages("lme4")
#install.packages("optimx")
#install.packages("minqa")
#install.packages("nloptr")
#install.packages("RcppEigen")
#install.packages("MatrixModels")
library(lme4)
library(optimx)
library(minqa)
library(nloptr)
library(RcppEigen)
library(MatrixModels)
lmerControl(optCtrl = list(maxfun = 2e5))
model_logistic_regression_factor_FODs <- glmer(categorical_FOD_FODs ~ 1 + Condition_FODs * Language_used_FODs +
(1|subject_FODs),
data=indvar_FODs, family = binomial(link = "logit"))
summary(model_logistic_regression_factor_FODs)
> summary(model_logistic_regression_factor_FODs)
Generalized linear mixed model fit by maximum likelihood (La
model <- multinom(categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + Condition_FODs*Language_used_FODs, indvar_FODs)
summary(model)
model1 <- lmer(categorical_FOD_FODs~Condition_FODs*Gender_FODs + (1|subject_FODs),
data=indvar_FODs,
control = lmerControl(
optimizer ='bobyqa', optCtrl=list(method='nlminb', maxfun=6000000), calc.derivs = FALSE))
summary(model1)
coef(model1)
fixef(model1)
sigma(model1)
sigma(model1)
sigma(model1)
sigma(model1)
sigma(model1)
sigma(model1)
sigma(model1)
sigma(model1)
sigma(model1)
sigma(model1)
sigma(model1)
sigma(model1)
sigma(model1)
sigma(model1)
sigma(model1)
sigma(model1)
sigma
fit <- glm(categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + Condition_FODs*Language_used_FODs,
+ family = multinomial("identity"), data = indvar_FODs)
glmer('categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + (1 | subject_FODs)',
+ data = indvar_FODs, family = binomial)
boundary (singular) fit: see help('isSingular')
model1 <- lme(FOD_FODs ~ 0 + Condition_FODs*Age_FODs, random=~1|subject_FODs,
na.action=na.omit,
method = "ML", data=indvar_FODs,control = lmeControl(opt = "optim"))
ano <- anova(model1)
round(ano, 3)
numDF denDF F-value p-value
Condition_FODs 4 3139 0.347 0.846
Age_FODs 15 3139 0.084 1.000
Condition_FODs:Age_FODs 45 3139 0.046 1.000
eta_squared(anova(model1), partial=TRUE)
# Effect Size for ANOVA
Parameter | Eta2 (partial) | 95% CI
-------------------------------------------------------
Condition_FODs | 4.42e-04 | [0.00, 1.00]
Age_FODs
model1_singular <- try(glmer(categorical_FOD_FODs ~ Condition_FODs + developmentaldisorder_FODs + Condition_FODs:developmentaldisorder_FODs + (1|subject_FODs),
data = indvar_FODs,
family = "binomial",
control = glmerControl(optimizer = "bobyqa")), TRUE)
if (class(model1_singular) == "try-error") {
model1 <- glmer(categorical_FOD_FODs ~ Condition_FODs + developmentaldisorder_FODs + Condition_FODs:developmentaldisorder_FODs + (1|subject_FODs),
data = indvar_FODs,
family = "binomial",
control = glmerControl(optimizer = "Nelder_Mead"))
}
else {
model1 <- model1_singular
}
fix invalid codeWed, 24 May 2023 ## Simple Math
try(m1 <- clm(categorical_FOD_FODs ~ 0 + Condition_FODs:Language_used_FODs, data = indvar_FODs), silent = TRUE)
indvar_FODs$categorical_FOD_FODs = factor(indvar_FODs$categorical_FOD_FODs)
m1 <- clm(categorical_FOD_FODs ~ 0 + Condition_FODs:Language_used_FODs - 1, data = indvar_FODs)
> head(data)
categorical_FOD_FODs Condition_FODs Language_used_FODs
1 D B Turkish
2 D D German
3 D D German
4 D D German
5 D A Turkish
6 D C Turkish
APA style is a style that is used for writing papers or journal articles in the social sciences.
For the chi-square test, we get $X^{2}(4, N = 3321) = 1.1716, p = .8828$.
For the multinomial logit model, we find that $A$ has a significant effect on the outcome, $B$ and $C$ do not ($p < .05$).
fix invalid codeWed, 24 May 2023 ## Credits
The package author is [Benjamin Wolfe](https://github.com/benjamin-wolfe). The package was inspired by the [Academic](https://github.com/gcushen/hugo-academic) theme for Hugo.
library(nnet)
fit.mlr_FODs <- multinom(
categorical_FOD_FODs~1,
data=indvar_FODs,
family = "binomial"(link = 'logit')
)
summary(fit.mlr_FODs)
# http://www.ats.ucla.edu/stat/r/dae/mlogit.htm
# http://www.ats.ucla.edu/stat/r/faq/mlogit.htm
library(mlogit)
mlogit(categorical_FOD_FODs ~ 1 | Condition_FODs + Language_used_FODs , data = indvar_FODs, method = "nr", print.level = 0)
#> Warning: The following arguments have been renamed: to --> data
#> Warning: The following arguments have been renamed: to --> data
#>
#> Call:
#> mlogit(formula = categorical_FOD_FODs ~ 1 | Condition_FODs +
#> Language_used_FODs, data = indvar_FODs, method = "nr",
#> print.level = 0)
#>
#> Coefficients:
#> (Intercept) Condition_FODsA Condition_FODsB
#
x <- 1
if (x) {
x <- 2
}
x
m1 <- clm(categorical_FOD_FODs ~ 0 + Condition_FODs:Language_used_FODs - 1, data = indvar_FODs)
fix invalid codeWed, 24 May 2023 ## Citation
To cite `statistics` in publications use:
## The results section
This section is usually the most difficult one, because you have to report the results in a clear and concise manner. The results section is not just a list of the results and the statistics, it is also the place where you explain them.
**So, what are the main elements of a results section?**
- You start by reporting the main results of your study, using descriptive statistics and/or figures.
- You report the statistics that you used to assess the assumptions of your model.
- Following, you report the main statistics of your analyses (e.g., regression weights, f-values, p-values, etc.).
- Next, you explain, in the text, the main results and their interpretation.
- You use tables and figures to clarify the results and to make it easier for the reader to understand them.
**Example**
fix invalid codeWed, 24 May 2023 ### Reference
[1]
levels(indvar_FODs$Condition_FODs)
levels(indvar_FODs$Language_used_FODs)
levels(indvar_FODs$categorical_FOD_FODs)
fix invalid codeFri, 12 May 2023 More information about syntax can be found [here](https://help.github.com/articles/basic-writing-and-formatting-syntax/)
# Notes on the data
## Format
The data is in the form of a long table, with each observation representing a
subject with a row. Observations are grouped into "experiments", with each
observation belonging to one experiment. Each experiment has a certain number
of columns, and each column is a different question or measurement.
The data file has three columns containing information about the experiment, and
some number of columns containing information about the subjects. The three
columns are:
1. `experiment`: The name of the experiment. All the rows that belong to the
same experiment have the same value in this column.
2. `experiment_name`: The name of the experiment in the form of a sentence.
3. `experiment_description`: A description of the experiment.
## Missing data
If a subject does not answer a question, the corresponding cell in the table
contains the value `missing`.
avfourfods +
geom_hline(yintercept = 0, color = '#000000', linetype = 'solid', size = 2) +
scale_fill_manual(values = c('#82c4fa','#fda02a','#00d2ff','#ff0033'))
mod0 <- glmer(categorical_FOD_FODs~Condition_FODs*Gender_FODs+
(1|subject_FODs),
na.action=na.omit,
data=indvar_FODs,
control = lmeControl(opt = "optim"))
mod <- glmer('categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + (1 | subject_FODs)',
data = indvar_FODs, family = binomial)
glmer(categorical_FOD_FODs~Condition_FODs*Gender_FODs+
(1|subject_FODs),
na.action=na.omit,
data=indvar_FODs,
control = lmeControl(opt = "optim"))
data(diamonds)
lm(price ~ carat + clarity + cut, data = diamonds)
# Call:
# lm(formula = price ~ carat + clarity + cut, data = diamonds)
#
# Coefficients:
# (Intercept) carat clarity1 clarity2 clarity3 clarity4 clarity5 clarity6 clarity7 cut1 cut2 cut3 cut4
# -8.001e+03 -6.351e+03 -8.578e+03 -9.923e+03 -8.724e+03 -7.864e+03 -7.909e+03 -7.929e+03 -7.925e+03 -7.658e+03 -7.871e+03 -6.890e+03 -8.318e+03
# cut5
# -7.932e+03
posthoc <- glht(anova_result, linfct = mcp(Condition_SNB = "Tukey"))
summary(posthoc)
> summary(mod)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: binomial ( logit )
Formula: categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + (1 | subject_FODs)
Data: indvar_FODs
AIC BIC logLik deviance df.resid
9493.3 9617.2 -4739.7 9479.3 4304
Scaled residuals:
Min 1Q Median 3Q Max
-1.8235 -0.7290 -0.6139 0.7590 1.8492
Random effects:
Groups Name Variance Std.Dev.
subject_FODs (Intercept) 0.07844 0.2815
Number of obs: 4320, groups: subject_FODs, 160
Fixed effects:
Estimate Std. Error z value Pr(
python
# The function add takes two arguments a and b and returns the sum of a and b.
def add(a, b):
return a + b
The ideal test to conduct is likelihood ratio tests of cumulative link models.
fix invalid codeFri, 05 May 2023 ## Testing
Examples of testing include:
> indvar_FODs$Condition_FODs<-factor(indvar_FODs$Condition_FODs, levels=c("A", "B", "C", "D"), labels=c("A", "B", "C", "D"))
> indvar_FODs$developmentaldisorder_FODs<-factor(indvar_FODs$developmentaldisorder_FODs, levels=c("No", "Yes"), labels=c("No", "Yes"))
> indvar_FODs$categorical_FOD_FODs<-factor(indvar_FODs$categorical_FOD_FODs, levels=c("ascending", "identity", "descending"), labels=c("ascending", "identity", "descending"))
> model1 <- glmer(categorical_FOD_FODs ~ Condition_FODs + developmentaldisorder_FODs + Condition_FODs:developmentaldisorder_FODs + (1|subject_FODs),
+
library(ggplot2)
# I have a dataset that has four columns, x, y, fill, and order.
# The x and y are are categorical, the fill is continious, and the order is a integer
# I want to change the order of the bars, manually as I wish. How can I do this?
ggplot(df, aes(x = x, y = y, fill = fill)) +
geom_bar(stat = 'identity') +
scale_fill_gradient(low = 'blue', high = 'red') +
guides(fill = guide_colorbar(order = 1)) #<- I want to change the order here
fix invalid codeFri, 21 Apr 2023 test code
the model, its results, and the normality test
fix invalid codeMon, 01 May 2023 Model_DD <- glmer(cat_SNB_new ~ Condition_SNB * Gender_SNB + (1 + Condition_SNB | subject_SNB),
data = s_complete,
family = binomial,
control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 20000)))
indvar_FODs_no_missing <- indvar_FODs[complete.cases(indvar_FODs), ]
## [3.2] Multinomial Logistic Regression
Condition_FODs
L1_or_not_FODs
Condition_FODs:L1_or_not_FODs
mod <- glmer('categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + (1 | subject_FODs)', data = indvar_FODs, family=binomial)
mod
all.equal(model$coefficients, coef(model))
tidy_fod_mult <- tidy(model, conf.int = TRUE)
tidy_fod_mult
mod <- glmer('categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + (1 | subject_FODs)',
data = indvar_FODs, family = binomial)
## https://stats.idre.ucla.edu/r/dae/multinomial-logistic-regression/
# Example of a multinomial logistic regression
#
# The data are from the example in Agresti (2013), Categorical Data Analysis, 3rd ed.,
# Table 10.2, p. 497.
#
# The data are from a survey to determine the preferred vacation site
# of 100 people, who were asked to pick their first, second, and third
# choice among the following four sites. The data are coded as
# 0 = no preference, 1 = first choice, 2 = second choice, 3 = third choice.
#
# The response is categorical, with four categories:
# 1 = Grand Canyon, 2 = Redwood Forest, 3 = Yosemite, 4 = Yellowstone.
# There are four predictors: age, sex, education, and income.
# The predictors are continuous.
#
# The model is
# logit(p[i,j]) = beta[0,j] +
fix invalid codeFri, 05 May 2023
fix invalid codeSat, 06 May 2023 ## @todo
- [x] Add R code
- [x] Add R error code
- [x] Add R output code
- [x] Add math
- [x] Add R math code
- [x] Add math R code
- [x] Add table
- [x] Add R table code
- [x] Add table R code
- [x] Add image
- [x] Add R image code
- [x] Add image R code
- [x] Add package import
- [x] Add R package import code
- [x] Add package import R code
- [x] Add package install
- [x] Add R package install code
- [x] Add package install R code
- [x] Add R code output
- [x] Add R code hidden output
- [ ] Add rmd to html
- [ ] Add rmd to pdf
- [ ] Add rmd to word
- [ ] Add rmd to markdown
m2 <- clm(categorical_FOD_FODs ~ 0 + Condition_FODs:Language_used_FODs - 1, random = ~1|subject, data = indvar_FODs)
summary(m2)
fix invalid codeFri, 05 May 2023 > table(s_complete$is_SNB)
False True
2892 716
> table(s_complete$diagnosis_SNB)
False True
2970 638
fix invalid codeSat, 06 May 2023 ## 移除文件夹中的文件，不包括子文件夹中的文件
The p-values of the ANOVA tests show that the effect of the main factors Condition_FODs and language_used_FODs_ind are not signifficant. The p-value of the interaction term is also not significant. This means that the model does not fit the data well. This is confirmed by the low partial eta squared values of the main factors and the interaction.
library("knitr")
library("kableExtra")
kable(pairwise, caption = "Pairwise Comparisons") %>%
kable_styling("striped") %>%
add_footnote("P value adjustment: tukey method for comparing a family of 4 estimates")
contrasts(df$condition) <- contr.sum
model <- lm(data ~ condition, data = df)
anova(model, type = "III")
fix invalid codeMon, 01 May 2023 Model_DD_1 <- glmer(cat_SNB_new ~ Condition_SNB * Gender_SNB + (1 | subject_SNB) + (1 | Gender_SNB:Condition_SNB),
data = s_complete,
family = binomial)
The results of the analysis are presented in table 1. According to the table, the effect of condition and age on FODs was not significant. In addition, there was also no significant interaction between condition and age. The partial eta-square indicated that condition and age have a trivial effect on FODs (eta-squared < 0.01) and the R2GLMM indicated that condition and age explained less than 1% of the variance of the data.
Table 1. Results of ANOVA
fix invalid codeMon, 01 May 2023 Model_DD <- glmer(cat_SNB_new ~ Condition_SNB * Gender_SNB + (1 + Condition_SNB | subject_SNB),
data = s_complete,
family = binomial,
control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 20000)))
print(pairwise, p.adj = "tukey")
My_Data <- read.csv('~/Desktop/My_Data.csv')
head(My_Data)
## Coding for a Mixed ANOVA
## Set up the model
My_Data$direction <- relevel(My_Data$direction, ref = "morning")
My_Data$Cond <- relevel(My_Data$Cond, ref = "Ay")
My_Data$colours <- relevel(My_Data$colours, ref = "blue")
My_Data.lm <- lm(direction ~ Cond * colours, data = My_Data)
## Model fit
anova(My_Data.lm)
fix invalid codeFri, 12 May 2023 ## expand_more
fix invalid codeMon, 01 May 2023 Model_DD_1 <- glmer(cat_SNB_new ~ Condition_SNB*Gender_SNB + (1+Condition_SNB|subject_SNB) + (1+Condition_SNB|Gender_SNB),
data = s_complete,
family = binomial)
fix invalid codeMon, 01 May 2023 > Model_DD_1 <- glmer(cat_SNB_new ~ Condition_SNB * Gender_SNB + (1 | subject_SNB) + (1 | Gender_SNB:Condition_SNB),
+ data = s_complete,
+ family = binomial)
boundary (singular) fit: see help('isSingular')
fix invalid codeWed, 10 May 2023 # Modeling - Mixed Effects
In this section we will try to predict whether a participant was classified as a native speaker of Dutch or English.
## Predicting Participant Native Language
> model <- glmer(cat_SNB_new ~ 0 + Condition_SNB * nativelang_SNB - 1+ (1 | subject_SNB),
+ data = s_complete, family = binomial)
> summary(model)
> anova(model)
> model_interactions <- glmer(cat_SNB_new ~ 0 + Condition_SNB * nativelang_SNB + (1 | subject_SNB),
+ data = s_complete, family = binomial)
> summary(model_interactions)
> anova(model_interactions)
## Model Fitting Diagnostics
> x <- cbind(s_complete$cat_SNB_new, s_complete$subject_SNB)
> x <- cbind(s_complete$cat_SNB_new, s_complete$subject_SNB)
fix invalid codeWed, 24 May 2023 > eta_squared <- eta_squared(model)
fix invalid codeWed, 24 May 2023 {r}
library(lsr)
etaSquared(model)
indvar_FODs<-read.csv(file="indvar_FODs.csv", header=TRUE, sep=",")
head(indvar_FODs)
#transforming categorical_FOD_FODs variable to factor variable
indvar_FODs$categorical_FOD_FODs<-as.factor(indvar_FODs$categorical_FOD_FODs)
#Treating categorical_FOD_FODs variable as dependent variable
mydata<-indvar_FODs[,c("categorical_FOD_FODs","Condition_FODs","Language_used_FODs","alltogether_FODs","sentence_FODs","verb_FODs","num_FODs")]
#logistic regression
logistic_regression = glm(categorical_FOD_FODs ~ Condition_FODs * Language_used_FODs * alltogether_FODs * sentence_FODs * verb_FODs *
library(multcomp)
plot(glht(model, linfct = mcp(Condition_FODs = "Tukey")))
fix invalid codeWed, 10 May 2023 str(s_complete)
# install.packages("lme4")
library(lme4)
logistic.model <- glmer(categorical_FOD_FODs~Condition_FODs + Language_used_FODs + Condition_FODs*Language_used_FODs + (1|subject_FODs), data=indvar_FODs, family=binomial)
summary(logistic.model)
fix invalid codeSat, 06 May 2023 1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, ... :
Model failed to converge: degenerate Hessian with 1 negative eigenvalues
2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, ... :
Model is nearly unidentifiable: large eigenvalue ratio
- Rescale variables?
fix invalid codeSat, 06 May 2023 fix invalid code:
> model_full2 <- glmer(cat_SNB_new ~ 0 + developmentaldisorder_SNB * Condition_SNB -1 + (1+developmentaldisorder_SNB|subject_SNB),
+ data = s_complete, family = binomial)
Warning messages:
1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
unable to evaluate scaled gradient
2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge: degenerate Hessian with 1 negative eigenvalues
fix invalid codeWed, 24 May 2023 library(car)
library(reshape2)
options(contrasts = c("contr.sum", "contr.poly"))
chisq.test(cont_table)
Pearson's Chi-squared test
data: cont_table
X-squared = 105.87, df = 6, p-value < 2.2e-16
please conduct post-hoc test in this data
m2 <- clm(categorical_FOD_FODs ~ 0 + Condition_FODs:Gender_quantised_FODs - 1,
+ random = ~1|subject_FODs, data = indvar_FODs)
> summary(m2)
formula: categorical_FOD_FODs ~ 0 + Condition_FODs:Gender_quantised_FODs - 1
data: indvar_FODs
link threshold nobs logLik AIC niter max.grad cond.H
logit flexible 3321 -2913.62 5839.25 6(0) 8.61e-13 4.7e+01
Coefficients:
Estimate Std. Error z value Pr(>|z|)
Condition_FODsA:Gender_quantised_FODs 0.12910 0.06843 1.887 0.0592 .
Condition_FODsB:Gender_quantised_FODs -0.06120
require(nnet)
data(iris)
m1 <- multinom(Species ~ Sepal.Length + Sepal.Width, data = iris)
m1
summary(m1)
m2 <- multinom(Species ~ Sepal.Length*Sepal.Width, data = iris)
m2
summary(m2)
fix invalid codeWed, 10 May 2023 R
model_full2 <- glmer(cat_SNB_new ~ 0 + developmentaldisorder_SNB * Condition_SNB -1 + (1+developmentaldisorder_SNB|subject_SNB),
data = s_complete, family = binomial)
I am writing a result section of a research paper. Can you please report these results in APA style?
model <- lmer(num_SNB ~ 0 + Condition_SNB*Language_used_SNB + (1|subject_SNB), data = s_complete)
anova(model)
Type III Analysis of Variance Table with Satterthwaite's method Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
Condition_SNB 10826.0 2706.50 4 145.6 340.1028 <2e-16 ***
Language_used_SNB 42.5 10.62 4 39.2 1.3345 0.2743
Condition_SNB:Language_used_SNB 57.7 4.81 12 3449.2 0.6039 0.8407
--- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> m2 <- clm(categorical_FOD_FODs ~ 0 + Condition_FODs:Gender_quantised_FODs - 1,
+ random = ~1|subject_FODs, data = indvar_FODs)
> summary(m2)
formula: categorical_FOD_FODs ~ 0 + Condition_FODs:Gender_quantised_FODs - 1
data: indvar_FODs
link threshold nobs logLik AIC niter max.grad cond.H
logit flexible 3321 -2913.62 5839.25 6(0) 8.61e-13 4.7e+01
Coefficients:
Estimate Std. Error z value Pr(>|z|)
Condition_FODsA:Gender_quantised_FODs 0.12910 0.06843 1.887 0.0592 .
Condition_FODsB:Gender_quantised_FODs -0
#creating the model
mod <- glmer('categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + (1 | subject_FODs)',
data = indvar_FODs, family = binomial)
summary(mod)
lmer(FOD_FODs ~ Condition_FODs + (1|subject_FODs), data = indvar_FODs)
fix invalid codeSat, 06 May 2023 ### R Markdown
R Markdown is an authoring format that enables easy creation of dynamic documents, presentations, and reports from R. It combines the core syntax of markdown (an easy-to-write plain text format) with embedded R code chunks that are run so their output can be included in the final document.
- [R Markdown](RMarkdown.ipynb)
- [R Markdown Cheat Sheet](RMarkdown_Cheat_Sheet.ipynb)
### R Notebooks
The R Notebook is an R Markdown document with chunks that can be executed independently and interactively, with output visible immediately beneath the input. R Notebooks are a new feature in RStudio 1.0, and are currently available as a preview release.
- [R Notebooks](RNotebooks.ipynb)
- [R Notebooks Cheat Sheet](RNotebooks_Cheat_Sheet.ipynb)
### Shiny
Shiny is an R package that makes it easy to build interactive web apps straight from R. You can host standalone apps on a webpage or embed them in R Mark
fix invalid codeMon, 24 Apr 2023 av_data1_all<-s_complete %>%
group_by(Gender_SNB, num_SNB)%>%
dplyr::summarize(prop=mean(num_SNB), synch = sum(num_SNB),
n = n(), .groups='drop')
The difference between these two tests is that the first is a fixed effect test, and the second is a mixed effects test (random effect).
The first test is a multinomial logistic regression. We are testing whether there is a significant interaction between condition and gender on a categorical DV that is a factor of FOD. The results show that the interaction between condition and gender is not significant.
The second test is a mixed effects multinomial regression. We are testing whether there is a significant interaction between condition and handedness on a categorical DV that is a factor of FOD. The results show that the interaction between condition and handedness is significant.
Since both are multinomial logistic regressions, we have to use a multinomial test to determine whether there is a significant main effect of condition on FOD. To do this we tested the effect of condition on FOD using a multinomial logistic regression. The results of this test showed that there was not a significant main effect of condition on FOD.
> p_value <- shapiro.test(resid(model1))$p.value
> round(p_value, 3)
[1] 0
> check_singularity(model1)
[1] FALSE
> ano <- anova(model1)
> round(ano, 3)
numDF denDF F-value p-value
Condition_FODs 4 3262 0.402 0.807
language_used_FODs_ind 4 40 0.008 1.000
Condition_FODs:language_used_FODs_ind 12 3262 0.053 1.000
> library(predictmeans)
> eta_squared(anova(model1), partial=TRUE)
# Effect Size for ANOVA
Parameter | Eta2 (partial) | 95% CI
---------------------------------------------------------------------
Condition_FODs | 4.93e-04 | [0.00, 1.00]
language_used_FODs_ind
s <- read.csv('s.csv', header = TRUE, sep = ',')
s_complete <- s[complete.cases(s$num_SNB), ]
summary(s_complete$num_SNB)
aggregate(num_SNB ~ Condition_SNB, data = s_complete, FUN = summary)
sd(s_complete$num_SNB)
count <- table(s_complete$num_SNB)
fix invalid codeFri, 05 May 2023 ### `X:`
>
fix invalid codeSat, 06 May 2023 anova_table <- Anova(model, type = "III")
library(report)
report(anova_table)
model1 <- glmer(categorical_FOD_FODs~Condition_FODs*Gender_FODs,random=~1|subject_FODs,
na.action=na.omit,
method = "ML",
data=indvar_FODs,
control = lmeControl(opt = "optim"))
mod <- glmer('categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + (1 | subject_FODs)',
+ data = indvar_FODs, family = binomial)
fix invalid codeMon, 24 Apr 2023 r
av_data1_all<-s_complete %>%
group_by(Gender_SNB, num_SNB)%>%
dplyr::summarize(prop=mean(num_SNB), synch = sum(num_SNB),
n = n(), .groups='drop')
d<-
av_data1_all %>% ggplot( aes(x=num_SNB, y=prop, color=Gender_SNB, group=Gender_SNB))+
geom_smooth(av_data1_all,method=glm, method.args= list(family = binomial(logit)),
mapping = aes(y=prop,group=Gender_SNB),linetype="dashed", position ="identity",se = FALSE)+ #mapping = aes(y=prop,group=Age, col='average'),linetype="dashed", position ="identity",se = FALSE)+
geom_smooth(method=glm, method.args
fit4 <- glm(categorical_FOD_FODs ~ 0 + Condition_FODs:nativelang_FODs - 1,
family = binomial(link = 'logit'), data = indvar_FODs)
summary(fit4)
cont_table <- data.frame(
table(gender, voting)
)
chisq_result <- chisq.test(cont_table)
chisq_result
fix invalid codeSat, 06 May 2023 ## Installation
You can install the released version of fixr from [CRAN](https://CRAN.R-project.org) with:
fix invalid codeSat, 06 May 2023 > anova_table <- Anova(model, type = "III")
> library(report)
> report(anova_table)
<table>
<tr>
<th colspan="4">Anova Table (Type III tests)</th>
</tr>
<tr>
<td></td>
<th>SS</th>
<th>DF</th>
<th>MS</th>
</tr>
<tr>
<th>Girth</th>
<td>206.9</td>
<td>31</td>
<td>6.68</td>
</tr>
<tr>
<th>Volume</th>
<td>12305.1</td>
<td>1</td>
<td>12305.1</td>
</tr>
<tr>
<th>Residuals</th>
<td>2051.6</td>
<td>29</
if (summary(model1)$isSingular) {
model1 <- glmer(categorical_FOD_FODs ~ developmentaldisorder_FODs + Condition_FODs + developmentaldisorder_FODs: Condition_FODs + (1|subject_FODs),
data = indvar_FODs,
family = "binomial",
control = glmerControl(optimizer = "bobyqa"))
}
fix invalid codeMon, 01 May 2023 Model_DD <- glmer(cat_SNB_new ~ Condition_SNB * Gender_SNB + (1 + Condition_SNB | subject_SNB),
data = s_complete,
family = binomial,
control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 20000)))
fix invalid codeSat, 06 May 2023 model_full2 <- glmer(cat_SNB_new ~ 0 + developmentaldisorder_SNB * Condition_SNB (1+developmentaldisorder_SNB|subject_SNB),
data = s_complete, family = binomial)
fix invalid codeMon, 24 Apr 2023
lmer_model <- lmer(FOD_FODs ~ Condition_FODs + (1|subject_FODs), data = indvar_FODs)
AIC is ideal for a model with a smaller number of parameters as it penalizes complexity.
fix invalid codeMon, 01 May 2023 [1] "Model_DD"
[1] "Model_DD:Cat_SNB ~ Condition_SNB + Gender_SNB + Condition_SNB:Gender_SNB (1|subject_SNB)"
Warning message:
In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.172378 (tol = 0.002, component 1)
> fix invalid code:
def add(a, b):
return a - b
fix invalid codeTue, 18 Apr 2023
library(lsmeans)
library(emmeans)
clm1 <- clm(master_data$categorical_FODs ~ master_data$Condition_FODs:master_data$Language_used_FODs)
emmeans(clm1, infer = c(TRUE, TRUE), pairwise ~ Condition_FODs:Language_used_FODs)
fix invalid codeMon, 01 May 2023 Model_DD <- glmer(cat_SNB_new ~ Condition_SNB * Gender_SNB + (1 + Condition_SNB | subject_SNB),
data = s_complete,
family = binomial,
control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 20000)))
# You can change the data to this
df <- data.frame(categorical_FOD_FODs = rep(c("ascending", "descending"),
each = 8),
Condition_FODs = rep(c("Verbal", "Vm", "Manual", "Mv"),
each = 2),
Frequency = c(100.0,
98.4,
94.4,
97.7,
0.0,
1.6,
5.6,
2.3))
df
mod <- glmer('categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + (1 | subject_FODs)', data = indvar_FODs,
family=binomial)
boundary (singular) fit: see help('isSingular')
mod
mod <- glmer('categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + (1 | subject_FODs)',
data = indvar_FODs, family = binomial)
The test is used to compare two models: the null model, and the model with changes. The p-value is used to determine whether to keep the changes in the model.
A p-value of less than 0.05 indicates that the changes improve the model fit.
A p-value of more than 0.05 indicates that the changes do not improve the model fit.
indvar_FODs$Condition_FODs <- as.factor(indvar_FODs$Condition_FODs)
indvar_FODs$Language_used_FODs <- as.factor(indvar_FODs$Language_used_FODs)
m1 <- clm(categorical_FOD_FODs ~ 0 + Condition_FODs:Language_used_FODs - 1, data = indvar_FODs)
summary(m1)
library(lme4)
df$Condition_FODs<-factor(df$Condition_FODs, levels=c("A", "B", "C", "D"), labels=c("A", "B", "C", "D"))
df$developmentaldisorder_FODs<-factor(df$developmentaldisorder_FODs, levels=c("No", "Yes"), labels=c("No", "Yes"))
df$categorical_FOD_FODs<-factor(df$categorical_FOD_FODs, levels=c("ascending", "identity", "descending"), labels=c("ascending", "identity", "descending"))
model1 <- glmer(categorical_FOD_FODs ~ Condition_FODs + developmentaldisorder_FODs + Condition_FODs:developmentaldisorder_FODs + (1|subject_FODs), data = df, family = "binomial", control = glmerControl(optimizer = "bobyqa"))
summary(model1)
model_comparison <- emmeans(model, "Language_used_FODs")
summary(model_comparison, test = "bonferroni")
This is for code blocks
### Task 1
Write a program that takes a number from the user and prints the result to check if it is a prime number.
What is a prime number?
A prime number is a natural number greater than 1 that has no positive divisors other than 1 and itself.
fix invalid codeMon, 01 May 2023 Model_DD_1 <- glmer(cat_SNB_new ~ Condition_SNB*Gender_SNB + (1+Condition_SNB|subject_SNB) + (1+Condition_SNB|Gender_SNB),
data = s_complete,
family = binomial)
fix invalid codeSat, 06 May 2023 > anova_table <- Anova(model, type = "III")
> report(anova_table)
fix invalid codeMon, 01 May 2023 fix invalid code:
Model_DD <- glmer(cat_SNB_new ~ Condition_SNB * Gender_SNB + (1 + Condition_SNB | subject_SNB),
data = s_complete,
family = binomial,
control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 20000)))
model <- glmmTMB(num_SNB ~ Condition_SNB * Gender_SNB + (1 | subject_SNB),
data = s_complete,
family = poisson(),
control = glmmTMBControl(optimizer = "nlminb", optArgs = list(maxiter = 2000)),
verbose = FALSE)
# Solution 1
# install.packages('lme4')
library(lme4)
# Solution 2
m1 <- clm(categorical_FOD_FODs ~ 0 + Condition_FODs:Language_used_FODs, data = indvar_FODs)
# Solution 3
m1 <- clm(categorical_FOD_FODs ~ 0 + Condition_FODs:Language_used_FODs, data = indvar_FODs)
fix invalid codeFri, 12 May 2023 > dotcounting_a <- table(s_complete$Dotcounting_SNB)
> dotcounting_a
A-B-C-D A-C-B-D A-D-B-C B-A-C-D B-A-D-C B-C B-C-D-A D-C-B-A
2711 80 80 159 80 80 80 159
> dotcounting_a <- dotcounting_a[!is.na(dotcounting_a)]
> dotcounting_a
A-B-C-D A-C-B-D A-D-B-C B-A-C-D B-A-D-C B-C B-C-D-A D-C-B-A
2711 80 80 159 80 80 80 159
> dotcounting_a <- as.numeric(dotcounting_a)
> dotcounting_a
[1] 2711 80
# multinomial logistic regression with interaction effect
model_mlr_interaction <- mlr(categorical_FOD_FODs ~ Condition_FODs * Language_used_FODs, data=indvar_FODs, family=multinomial("identity"))
summary(model_mlr_interaction)
fix invalid codeMon, 01 May 2023 Model_DD <- glmer(cat_SNB_new ~ Condition_SNB * Gender_SNB + (1 + Condition_SNB | subject_SNB),
+ data = s_complete,
+ family = binomial,
+ control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 20000)))
boundary (singular) fit: see help('isSingular')
posthoc_FODs <- glht(model, linfct = mcp(Condition_FODs = "Tukey", Language_used_FODs = "Tukey", Condition_FODs*Language_used_FODs = "Tukey"))
summary(posthoc_FODs)
cont_table <- table(data$Sex, data$Survived)
chisq_result <- chisq.test(cont_table)
chisq_result
### ANOVA
# One-way ANOVA
anova(lm(Performance ~ Condition, data=data))
# Two-way ANOVA
anova(lm(Performance ~ Condition*Language_used, data=data))
# > library(car)
# > with(mydata, Anova(lm(categorical ~ Condition*Language_used, data=data), type=3))
# Analysis of Deviance Table
#
# Model 1: categorical ~ Condition + Language_used + Condition:Language_used
# Model 2: categorical ~ Condition * Language_used
# Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F)
# 1 336 464.24
# 2 330 463.75 6 0.4935 0.288 0.917
#
# Analysis of Deviance Table
#
# Model 1: categorical ~ Condition + Language_used + Condition:Language_used
# Model 2: categorical ~ Condition * Language_used
# Res.Df Res.
sum(!is.na(indvar_FODs$categorical_FOD_FODs) & !is.na(indvar_FODs$Gender_quantised_FODs))
fix invalid codeMon, 24 Apr 2023 Gender_SNB num_SNB prop synch n
<fct> <dbl> <dbl> <int> <int>
1 male 7 7 0.35 49 1
2 male 42 42 0.35 49 1
3 male 42 42 0.35 49 1
4 male 42 42 0.35 49 1
5 male 42 42 0.35 49 1
6 male 42 42 0.35 49 1
7 male 42 42 0.35 49 1
8 male 42 42 0.35 49 1
9 male 42 42 0.35 49 1
10 male 42 42 0.35 49 1
# … with 21 more rows
The model1 is the name of the model created using the lme function. The lme function takes a number of arguments. The first two arguments are the response variable and the predictor variables. In this model the response variable is categorical_FOD_FODs and the predictor variables are Condition_FODs*Gender_FODs. The random argument represents the random effects of the model and in this case, the random effects are a random intercept for each participant and a random slope for each participant and each chunk condition. The random effects are nested within the participant and the chunk conditions are nested within each participant. The na.action argument specifies how to deal with missing values, in this case missing values will be removed from the model. The method argument specifies the fitting method to be used, in this case maximum likelihood estimation (ML). The data argument specifies the dataframe from which the variables are taken, in this case the dataframe is called long1.RSPL. The control argument specifies the options to be used in the model, in this case the optim argument specifies the optimization algorithm to be used, in this case a quasi-Newton method.
FODs_Logistic <- multinom(categorical_FOD_FODs ~ Condition_FODs*Language_used_FODs, indvar_FODs)
aov.m <- aov(direction ~ Cond + colours + Cond:colours + Error(colours), data=indvar_FODs_new)
summary(aov.m)
def model_msg(data, family, control):
return glmer(cat_SNB_new ~ 0 + Condition_SNB * Dotcounting_SNB - 1 + (1 | subject_SNB),
data = data,
family = family,
control = control)
# model with optimizer argument
mod_FODs_optimizer <- glmer(y ~ x1 * x2 + (1|subject_FODs), family = binomial, data = indvar_FODs,
control = glmerControl(optimizer = "bobyqa"))
# print model
summary(mod_FODs_optimizer)
mod_FODs_optimizer <- glmer(y ~ x1 * x2 + (1|subject_FODs), family = binomial, data = indvar_FODs,
control = glmerControl(optimizer = "bobyqa"))
multinomial <- function(link)
{
switch(link, identity = "multinomial.identity", logit = "multinomial.logit", probit = "multinomial.probit",
cauchit = "multinomial.cauchit")
}
glm(categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + Condition_FODs*Language_used_FODs,
family = multinomial("identity"), data = indvar_FODs)
The model uses the lme function from the nlme package to construct a mixed model. The function will output a model with a random intercept for the subject ID and a random slope for the subject ID and the condition.
The model has a residual deviance of 4598.5 and 3313 degrees of freedom, which is greater than the null deviance of 4603.3 and 3320 degrees of freedom.
The intercept is 0.02715 with a standard error of 0.09178.
The condition A is 0 (or omitted).
The condition B is -0.17148 with a standard error of 0.9832, which is significantly different from 0 (z = -1.744, p < .05).
The condition C is -0.14466 with a standard error of 0.9814, which is significantly different from 0 (z = -1.474, p < .05).
The condition D is -0.12076 with a standard error of 0.9831, which is significantly different from 0 (z = -1.228, p < .05).
The language German is 0.02854 with a standard error of 0.17759, which is not significantly different from 0 (z = 0.161, p > .05).
The language Hungarian is 0.16117 with a standard error of 0
function that:
library(nortest)
ad.test(x, alternative = c("two.sided", "less", "greater"),
exact = FALSE, tol = 1e-08, maxpts = 500000,
abseps = 1e-05, releps = 1e-08)
model <- glm(categorical_FOD_FODs ~ Condition_FODs * Language_used_FODs, family=multinomial("identity"), data=mydata)
The "> >" is a mistake.
model1 <- glmer(categorical_FOD_FODs ~ Condition_FODs + developmentaldisorder_FODs + Condition_FODs:developmentaldisorder_FODs + (1|subject_FODs),
data = indvar_FODs, family = "binomial", control = glmerControl(optimizer = "bobyqa"))
boundary (singular) fit: see help('isSingular')
fix invalid codeMon, 24 Apr 2023 s_data <- s_data %>%
group_by(Gender_SNB, num_SNB) %>%
dplyr::summarize(prop=mean(num_SNB), synch = sum(num_SNB),
n = n(), .groups='drop')
{r}
# first calculate the effect size
library(predictmeans)
model1_eta <- eta_squared(anova(model1), partial = TRUE)
round(model1_eta, 3)
# now report in APA style
apa.table(model1_eta,
p.value = ano$p.value,
stars = c("ns","ns","ns"),
digits = c(3,3,3),
custom.stars = TRUE,
custom.stars.table = c("", "*", "**", "***"),
star.symbol = "",
star.fontface = "plain",
star.cex = 1,
border = "topbottom",
caption = "",
align = "c",
booktab = TRUE,
space = "",
only.means = FALSE,
omit.empty.rows = TRUE,
print.results = "none",
file = NULL,
table.placement = NULL,
library(lme4)
mod0 <- glmer(categorical_FOD_FODs~Condition_FODs*Gender_FODs+
(1|subject_FODs),
na.action=na.omit,
data=indvar_FODs,
control = lmeControl(opt = "optim"))
Warning in glmer(categorical_FOD_FODs ~ Condition_FODs * Gender_FODs + (1 | :
Use control=glmerControl(..) instead of passing a list of class “list”
Error in (function (optimizer = c("bobyqa", "Nelder_Mead"), restart_edge = FALSE, :
argument 12 matches multiple formal arguments
library(ordinal)
m2 <- clm(categorical_FOD_FODs ~ 0 + Condition_FODs:Gender_FODs - 1,
random = ~1|subject, data = indvar_FODs)
> summary(m2)
formula: categorical_FOD_FODs ~ 0 + Condition_FODs:Gender_FODs - 1
data: indvar_FODs
link threshold nobs logLik AIC niter max.grad cond.H
logit flexible 3321 -2910.89 5847.79 6(0) 1.75e-12 1.4e+03
Coefficients: (1 not defined because of singularities)
Estimate Std. Error z value Pr(>|z|)
Condition_FODsA:Gender_FODsFemale 0.32677 0.33910 0.964 0.335
Condition_FODsB:Gender_FODsFemale 0.03530 0
fix invalid codeMon, 24 Apr 2023 {r, eval=FALSE,echo=TRUE}
d<-
av_data1_all %>% ggplot( aes(x=num_SNB, y=prop, color=Age_SNB, group=Age_SNB))+
geom_smooth(av_data1_all,method=glm, method.args= list(family = binomial(logit)),
mapping = aes(y=prop,group=Age_SNB),linetype="dashed", position ="identity",se = FALSE)+ #mapping = aes(y=prop,group=Age, col='average'),linetype="dashed", position ="identity",se = FALSE)+
geom_smooth(method=glm, method.args= list(family = binomial(logit)),
se = FALSE) +
geom_point() +
#geom_point(aes(shape = block)) +
xlab('Delay (ms)') + ylab('Prop. of Synch Responses
model1 fits a mixed effects model with the categorical_FOD_FODs as the dependent variable and Condition_FODs and Genders_FODs as fixed effects. Random intercepts were fitted at the level of the Participant and the
Variable (Ch)
class(indvar_FODs$categorical_FOD_FODs)
[1] "integer"
fix invalid codeSat, 06 May 2023 r
> report(anova_table)
fix invalid codeSat, 06 May 2023 > model_full2 <- glmer(cat_SNB_new ~ 0 + developmentaldisorder_SNB * Condition_SNB -1 + (1|subject_SNB),
+ data = s_complete, family = binomial)
Warning messages:
1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
unable to evaluate scaled gradient
2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge: degenerate Hessian with 1 negative eigenvalues
3: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model is nearly unidentifiable: very large eigenvalue
- Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
- Rescale variables?
fix invalid codeFri, 05 May 2023 def add(a, b):
return a + b
multcomp::glht(model1, linfct = mcp(Condition_FODs = "Tukey"))
options(warn=-1)
library(car)
library(lme4)
library(lmerTest)
library(lmtest)
library(multcomp)
library(lsmeans)
library(plyr)
library(nlme)
library(afex)
library(emmeans)
library(effsize)
library(compute.es)
library(multcompView)
library(ggplot2)
library(sjPlot)
library(texreg)
fix invalid codeWed, 24 May 2023 library(model.summary)
eta_squared <- etaSquared(model)
fix invalid codeWed, 24 May 2023 You can also use it to show an error:
model <- lm(categorical_FOD_FODs ~ Condition_FODs*Gender_FODs, data=indvar_FODs)
anova(model)
fix invalid codeWed, 24 May 2023 r
# Now, we can visualise the results of the ANOVA with eta2:
library(ggpubr)
# Create data frame with eta2 values
df1 <- as.data.frame(eta_squared(model))
df1
glmer(cat_SNB_new ~ 0 + Condition_SNB * Dotcounting_SNB - 1 + (1 | subject_SNB),
data = s_complete,
family = binomial,
control = glmerControl(optimizer = "bobyqa", tolPwrss = 1e-08))
optCtrl=list(method="bobyqa", optCtrl=list(maxfun=2e5, calc.derivs=TRUE))
mod_FODs_optimizer <- glmer(y ~ x1 * x2 + (1|subject_FODs), family = binomial, data = indvar_FODs,
control = glmerControl(optimizer = "bobyqa"))
fix invalid codeMon, 01 May 2023 Model_DD <- glmer(cat_SNB_new ~ Condition_SNB * Gender_SNB + (1 + Condition_SNB | subject_SNB),
data = s_complete,
family = binomial,
control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 20000)))
mod0 <- glmer(categorical_FOD_FODs~Condition_FODs*Gender_FODs+
(1|subject_FODs),
na.action=na.omit,
data=indvar_FODs,
control = lmeControl(opt = "optim"))
fix invalid codeFri, 12 May 2023 ### Exercise 4
ds.test <- shapiro.test(resid(model))
fix invalid codeWed, 24 May 2023 </details>
<details><summary>week 2</summary>
mod_FODs_optimizer <- glmer(y ~ x1 * x2 + (1|subject_FODs), family = binomial, data = indvar_FODs,
control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun=2e5)))
m1 <- clm(categorical_FOD_FODs ~ 0 + Condition_FODs:Language_used_FODs, data = indvar_FODs)
lme4::clm(categorical_FOD_FODs ~ Condition_FODs:Language_used_FODs, data = indvar_FODs)
> install.packages("lme4")
> library(lme4)
> logistic_model <- glmer(categorical_FOD_FODs ~ Condition_FODs * Language_used_FODs + (1|subject_FODs), family=binomial, data=indvar_FODs)
> summary(logistic_model)
> optimizer <- optimizerControl(maxit = 10000)
> logistic_model <- glmer(categorical_FOD_FODs ~ Condition_FODs * Language_used_FODs + (1|subject_FODs), family=binomial, data=indvar_FODs, control=optimizer)
> summary(logistic_model)
> fixef(logistic_model)
(Intercept) Condition_FODsB Condition_FODsC Condition_FODsD Language_used_FODsEnglish Language_used_FODsGerman Language_used_FODsHungarian Language_used_FOD
fix invalid codeFri, 12 May 2023 *Quick Tip:*
- To print a backtick (`), use `\``.
- To print a backtick (
fix invalid codeSat, 06 May 2023 # Reproducible Research
What is Reproducible Research?
- The ability to reproduce results from a research article
- Refers to both the computational workflows and the data
Why is Reproducible Research Important?
- Ensures the results are trustworthy
- Allows others to extend your work
- Allows others to check your work
How do I do Reproducible Research?
- Use tools such as Jupyter Notebooks
<img src="images/jupyter.png" width=300>
- Store your data in a version control system
<img src="images/github.png" width=300>
- Use version control to track changes to code
<img src="images/git.png" width=300>
# Jupyter Notebooks
<img src="images/jupyter.png" width=300>
- Notebooks are a means to combine code with descriptions and output
- Notebooks are a means to combine Markdown and code
The main benefits of Jupyter Notebooks are:
- The ability to combine code, documentation and
av_data1_all <- s_complete %>%
group_by(Gender_SNB, num_SNB) %>%
summarise(mean_num_SNB = mean(num_SNB),
sum_synch = sum(synch),
count = n()) %>%
ungroup()
mod0 <- glmer(categorical_FOD_FODs~Condition_FODs*Gender_FODs+
(1|subject_FODs),
na.action=na.omit,
data=indvar_FODs,
control = lmeControl(opt = "optim"))
fix invalid codeFri, 05 May 2023 {% endcomment %}
fix invalid codeFri, 05 May 2023 Model_MSG <- glmer(cat_SNB_new ~ 0 + Condition_SNB * Dotcounting_SNB - 1 + (1 | subject_SNB),
data = s_complete,
family = binomial,
control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e5)))
fix invalid codeWed, 24 May 2023 ### Step 5. Conduct pairwise comparisons (optional)
fix invalid codeWed, 24 May 2023 fix invalid code:
> library(car)
> etaSquared(model)
fix invalid codeMon, 24 Apr 2023 ### 3.2.2 Data visualization
> m2 <- clm(categorical_FOD_FODs ~ 0 + Condition_FODs:Language_used_FODs - 1, random = ~1|subject, data = indvar_FODs)
> summary(m2)
formula: categorical_FOD_FODs ~ 0 + Condition_FODs:Language_used_FODs - 1
data: indvar_FODs
link threshold nobs logLik AIC niter max.grad cond.H
logit flexible 3321 -2909.91 5861.83 6(0) 1.02e-12 2.3e+02
Coefficients: (1 not defined because of singularities)
Estimate Std. Error z value Pr(>|z|)
Condition_FODsA:Language_used_FODsEnglish 0.15985 0.16285 0.982 0.3263
Condition_FODsB:Language_used_FODsEnglish -0
indvar_FODs_new <- indvar_FODs %>%
filter(direction %in% c("morning", "noon", "night"))
glmer(direction ~ Cond + colours + Cond:colours +
(1|colours) + (1|Cond), data = indvar_FODs_new, family = "binomial")
multinomial <- function(link)
{
switch(link, identity = "multinomial.identity", logit = "multinomial.logit", probit = "multinomial.probit",
cauchit = "multinomial.cauchit")
}
fix invalid codeMon, 01 May 2023 r
## Model DD
Model_DD_1 <- glmer(cat_SNB_new ~ Condition_SNB*Gender_SNB + (1+Condition_SNB|subject_SNB) + (1+Condition_SNB|Gender_SNB),
data = s_complete,
family = binomial)
fix invalid codeSat, 06 May 2023 model_full2 <- glmer(cat_SNB_new ~ 0 + developmentaldisorder_SNB * Condition_SNB -1 + (1+developmentaldisorder_SNB|subject_SNB),
y <- as.numeric(s_complete$cat_SNB_new)
model <- lmer(y ~ Condition_SNB * Dotcounting_SNB + (1 | subject_SNB), data = na.omit(s_complete))
### this is the logistic regression
mod_logistic <- glmer(categorical_FOD_FODs ~ (1|subject_FODs) +
Condition_FODs * Language_used_FODs,
data = indvar_FODs, family = binomial)
mod_logistic2 <- glmer(categorical_FOD_FODs ~ (1|subject_FODs) +
Condition_FODs * Language_used_FODs +
(1|subject_FODs:Condition_FODs) +
(1|subject_FODs:Language_used_FODs),
data = indvar_FODs, family = binomial)
mod_logistic3 <- glmer(categorical_FOD_FODs ~ (1|subject_FODs) +
Condition_FODs * Language_used_FODs +
(1|subject_FODs:
I want to run in R the following command: plot(residuals(lm(y~x)) ~ fitted(lm(y~x)))
the result is:
### [Test 2]
fix invalid codeFri, 05 May 2023 # How to fix:
> Model_MSG <- glmer(cat_SNB_new ~ 0 + Condition_SNB * Dotcounting_SNB - 1 + (1 | subject_SNB),
+ data = s_complete,
+ family = binomial,
+ control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 5000)))
## example2:
> Model_MSG <- glmer(cat_SNB_new ~ 0 + Condition_SNB * Dotcounting_SNB - 1 + (1 | subject_SNB),
+ data = s_complete,
+ family = binomial,
+ control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 5000)))
## example3:
> Model_MSG <- glmer(cat_SNB_new ~ 0 + Condition_SNB * Dotcounting_SNB - 1 + (1 | subject_SNB),
+ data = s_complete,
+
fix invalid codeWed, 24 May 2023 fix invalid code:
> library(car)
> eta_squared <- etaSquared(model)
Error in etaSquared(model) : could not find function "etaSquared"
fix invalid codeWed, 24 May 2023 eta_squared <- eta_squared(model)
ggplot2:::ggplot2_colors
# get the colour codes used in ggplot
ggplot2::ggplot_build(avfourfods)$data[[1]]$colour
# get the colour codes in the same order to use in other places
avfourfods <- avfourfods + scale_fill_manual(values=c("#CC3333", "#3366FF", "#FF9900", "#33CC33"))
## This is a list of questions
1. This is a question
2. This is another question
3. This is a question
library(multcomp)
mod_glht <- glht(mod, linfct = mcp(Condition_FODs = "Tukey"))
summary(mod_glht)
library(lsmeans)
lsmeans(mod, ~Condition_FODs)
fix invalid codeFri, 05 May 2023 Model_MSG <- glmer(cat_SNB_new ~ 0 + Condition_SNB * Dotcounting_SNB - 1 + (1 | subject_SNB),
data = s_complete,
family = binomial,
control = glmerControl(optimizer = "Nelder_Mead"))
library(multcomp)
library(multcompView)
# load data
data(pwlf)
posthoc <- glht(pwlf$glht, linfct = mcp(factor(pwlf$condition) = "Tukey"))
# write table
multcompView::mctable(posthoc)
# Creating a model
model <- lmer(formula = direction ~ Cond * colours + (1|colours), data = indvar_FODs_new)
# I will add the confidence intervals in my plot to show the variability in the data.
ggplot(df, aes(x = categorical_FOD_FODs, y = Frequency_FODs, fill = Condition_FODs)) +
geom_bar(stat = "identity", position = position_dodge()) +
geom_errorbar(aes(ymin = Frequency_FODs - sd, ymax = Frequency_FODs + sd), position = position_dodge(.9), width = .25) +
scale_fill_manual(values = colors,
labels = c("Verbal", "Vm", "Manual", "Mv")) +
labs(
x = "FODs (direction)",
y = "Frequency of FODs",
title = "FODs values (ascending/descending) in four condition",
fill = "Conditions (renamed)"
) +
theme_minimal() +
theme(
fix invalid codeFri, 12 May 2023 ### Reporting Guidelines
> Please see [CONTRIBUTING](https://github.com/IBM/MAX-Object-Detector/blob/master/CONTRIBUTING.md) for more details on how to report issues and the process for submitting pull requests.
### LICENSE
[Apache 2.0](LICENSE)
[travis-image]: https://travis-ci.org/IBM/MAX-Object-Detector.svg?branch=master
[travis-url]: https://travis-ci.org/IBM/MAX-Object-Detector
[license-image]: https://img.shields.io/pypi/l/pipenv.svg
[license-url]: https://github.com/IBM/MAX-Object-Detector/blob/master/LICENSE
[codacy-image]: https://api.codacy.com/project/badge/Grade/0b6f950d5f2045b98b1e35a7e2a4b8d2
[cod
indvar_FODs$categorical_FOD_FODs <- as.factor(indvar_FODs$categorical_FOD_FODs)
m1 <- clm(categorical_FOD_FODs ~ 0 + Condition_FODs:Language_used_FODs - 1, data = indvar_FODs)
rs_complete <- read.csv("rs_complete.csv")
mean(rs_complete$num_RS)
summary(rs_complete$num_RS)
aggregate(num_RS ~ Condition_RS, data = rs_complete, FUN = summary)
sd(rs_complete$num_RS)
count <- table(rs_complete$num_RS)
fix invalid codeMon, 24 Apr 2023 av_data1_all<-s_complete %>%
group_by(Gender_SNB, num_SNB)%>%
dplyr::summarize(prop=mean(num_SNB), synch = sum(num_SNB),
n = n(), .groups='drop')
library(lme4)
mod <- glmer('categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + (1 | subject_FODs)', data = indvar_FODs,
family=binomial)
s_complete <- read.csv("C:/Users/simge/Desktop/R/s_complete.csv")
s_complete$num_SNB <- ifelse(s_complete$Condition_SNB == "SNB_positive", 1, 0)
summary(s_complete$num_SNB)
aggregate(num_SNB ~ Condition_SNB, data = s_complete, FUN = summary)
sd(s_complete$num_SNB)
count <- table(s_complete$num_SNB)
model1 <- glmer(categorical_FOD_FODs ~ Condition_FODs + developmentaldisorder_FODs + Condition_FODs:developmentaldisorder_FODs +
(1|subject_FODs), data = indvar_FODs, family = "binomial", control = glmerControl(optimizer = "bobyqa"))
library(lme4)
library(lmerTest)
# dependent variable
y <- as.factor(indvar_FODs$categorical_FOD_FODs)
# independent variables
x1 <- indvar_FODs$Condition_FODs
x2 <- indvar_FODs$Language_used_FODs
# model with optimizer argument
mod_FODs <- glmer(y ~ x1 * x2 + (1|subject_FODs), family = binomial, data = indvar_FODs,
control = glmerControl(optimizer = "bobyqa"))
# print model
summary(mod_FODs)
fix invalid codeSat, 06 May 2023 r
> anova_table <- Anova(model, type = "III")
> library(ez)
> ezANOVA(model)
$ANOVA
Effect DFn DFd F p p<.05 ges
1 (Intercept) 1 6 563.547681 0.000000e+00 * 0.9137616
2 Group2 2 12 1.763823 1.484842e-01 0.0322126
3 Group2:Group3 2 12 0.064462 9.308931e-01 0.0012488
4 Group3 2 12 0.066804 9.057730e-01 0.0012984
5 Residuals 18 96
$`Sphericity Corrections`
GGe HF eps
1 1.00000 0.924
2 1.02369 0.923 0.01
3 1.02444 0.922 0.01
4 1.02469 0
multinomial <- function(link) {
switch(link, identity = "multinomial.identity", logit = "multinomial.logit", probit = "multinomial.probit",
cauchit = "multinomial.cauchit")
}
glm(categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + Condition_FODs*Language_used_FODs,
family = multinomial("identity"), data = indvar_FODs)
fix invalid codeFri, 12 May 2023 ## Manuals
* [R Markdown: The Definitive Guide](https://bookdown.org/yihui/rmarkdown/)
* [R Markdown Cookbook](https://bookdown.org/yihui/rmarkdown-cookbook/)
* [R Markdown Reference Guide](https://bookdown.org/yihui/rmarkdown-reference/)
* [R Markdown: The Big Split](https://rmarkdown.rstudio.com/articles_split_files.html)
* [Package rmarkdown](https://cran.r-project.org/web/packages/rmarkdown/rmarkdown.pdf)
* [R Markdown Cheatsheet](https://www.rstudio.com/wp-content/uploads/2016/03/rmarkdown-cheatsheet-2.0.pdf)
* [R Markdown Cheat Sheet](https://www.rstudio.com/wp-content/uploads/2015/02/rmarkdown-cheatsheet.pdf)
## Examples
* [R Markdown Examples](https://
fix invalid codeWed, 24 May 2023 library(car)
# path to your data, and name of data file
getwd()
chisq_test <- function(cont_table) {
chisq_result <- chisq.test(cont_table)
return(chisq_result)
}
chisq_result <- chisq_test(cont_table)
print(chisq_result)
library(tidyverse)
library(stringr)
library(car)
library(stargazer)
library(lme4)
library(lmerTest)
library(lattice)
library(lsmeans)
library(multcomp)
library(MASS)
library(xtable)
library(psych)
library(corrplot)
library(ggpubr)
fix invalid codeSat, 22 Apr 2023 Error in `summary.lme`(model, cor = cor) :object 'model' not found
library(multcomp)
library(lsmeans)
#Generate data
data <- data.frame(
ID = 1:60,
Condition_SNB = c(rep("A", 30), rep("B", 30)),
SNB_Score = c(rnorm(30, 110, 10), rnorm(30, 100, 10))
)
#Perform the ANOVA
anova_result <- lm(SNB_Score ~ Condition_SNB, data)
#Perform the post hoc test
posthoc <- glht(anova_result, linfct = mcp(Condition_SNB = "Tukey"))
#Print out the table
print(summary(posthoc))
mod <- glmer('categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + (1 | subject_FODs)', data = indvar_FODs,
+ family=binomial)
boundary (singular) fit: see help('isSingular')
mod
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: binomial ( logit )
Formula: categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + (1 | subject_FODs)
Data: indvar_FODs
AIC BIC logLik deviance df.resid
4616.503 4671.476 -2299.252 4598.503 3312
Random effects:
Groups Name Std.Dev.
subject_FODs (Intercept) 1.403e-07
Number of obs: 3321, groups: subject_FODs,
fix invalid codeSat, 06 May 2023 model_full2 <- glmer(cat_SNB_new ~ 0 + developmentaldisorder_SNB * Condition_SNB -1 + (1+developmentaldisorder_SNB|subject_SNB),
data = s_complete, family = binomial)
In the function glm, the argument categorical_FOD_FODs is the dependent variable of the model, the argument data is the data frame of the independent variables, the argument family is the family of the dependent variable, the argument control is how the model is optimized, and the argument penalty is the method of regularization.
fix invalid codeSat, 06 May 2023 > anova_table <- Anova(model, type = "III")
> report(anova_table)
Error in `[[<-.data.frame`(`*tmp*`, "Mean_Square", value = numeric(0)) :
Ersetzung hat 0 Zeilen, Daten haben 3
# import the datafile
data <- read.csv('C:/Users/jessicasa/Documents/Courses/Spring_2019/STA_6414/Final_project/data/exp1_FODs.csv')
# dependent variable
y <- as.factor(data$categorical_FOD_FODs)
# independent variables
x1 <- data$Condition_FODs
x2 <- data$Language_used_FODs
# model with optimizer argument
mod_FODs <- glmer(y ~ x1 * x2 + (1|subject_FODs), family = binomial, data = data,
control = glmerControl(optimizer = "bobyqa"))
# print model
summary(mod_FODs)
library(tidyverse)
library(lme4)
FODs_direction <- read.csv("FODs_direction.csv", header = T)
ggplot(FODs_direction, aes(x = categorical_FOD_FODs, y = Frequency, fill = Condition_FODs)) +
geom_bar(stat = "identity", position = position_dodge()) +
scale_fill_manual(values = c("#F8766D", "#7CAE00", "#00BFC4", "#C77CFF"),
labels = c("Verbal", "Vm", "Manual", "Mv")) +
labs(
x = "FODs (direction)",
y = "Frequency of FODs",
title = "FODs values (ascending/descending) in four condition",
fill = "Conditions (renamed)"
) +
theme_minimal() +
theme(
plot.title = element_
m2 <- glmer('FOD_categorical ~ Condition + Language_used + (1 | subject)',
+ data = indvar, family = binomial)
fit <- lmer('FOD_categorical ~ Condition + Language_used + (1 | subject)',
+ data = indvar)
# R code:
writeTable(coefs, caption = "Linear regression coefficients for the relationship between Condition_FODs and FOD_FODs")
## 6.2.2 FODs_clm: raw data and z-scored data
model <- glmer(categorical_FOD_FODs ~ Condition_FODs + developmentaldisorder_FODs + (1 | subject_FODs),
family = binomial(link = 'logit'), data = indvar_FODs)
mod_FODs <- glmer(y ~ x1 * x2 + (1|subject_FODs), family = binomial, data = indvar_FODs,
control = glmerControl(optimizer = "bobyqa"))
fix invalid codeWed, 24 May 2023 Fix invalid code:
def add(a, b):
return a - b
fix invalid codeWed, 10 May 2023 model_full2 <- glmer(cat_SNB_new ~ 0 + developmentaldisorder_SNB * Condition_SNB -1 + (1+developmentaldisorder_SNB+Condition_SNB|subject_SNB),
data = s_complete, family = binomial)
Model_MSG <- glmer(cat_SNB_new ~ 0 + Condition_SNB * Dotcounting_SNB - 1 + (1 | subject_SNB),
data = s_complete,
family = binomial,
control = glmerControl(optimizer = "bobyqa", tolPwrss = 1e-08))
fit4 <- glm(categorical_FOD_FODs ~ 0 + Condition_FODs:nativelang_FODs - 1,
family = binomial(link = 'logit'), data = indvar_FODs)
summary(fit4)
The LRT test tells us that the first model is better than the second because the p-value is larger than 0.05.
fix invalid codeFri, 12 May 2023
library(lme4)
# create a dataframe that only needs the columns of interest
indvar_FODs_new <- indvar_FODs[,c("direction","Cond","colours")]
# run the mixed model
mixed_model <- lmer(direction ~ Cond * colours + (1|colours), data = indvar_FODs_new)
f <- function(x) {
if (x < 0) {
return(0)
} else {
return(x)
}
}
f(-1)
You can use [this link](https://stats.idre.ucla.edu/other/mult-pkg/faq/general/faq-how-do-i-report-results-from-a-logistic-regression/) to check how to report the results in APA style.
</div>
<!--
###
The function add takes two arguments a and b and returns the sum of a and b.
my.data <- data.frame(means, sd, se)
colnames(my.data) <- c("condition", "mean", "sd", "se")
#means, sd, se are vectors of the appropriate length with the respective values
kable(my.data, caption = "Reaction times", booktabs = T)
> model <- lmer(SNB_new ~ Condition_SNB * Dotcounting_SNB + (1 | subject_SNB), data = na.omit(s_complete))
> summary(model)
Linear mixed model fit by REML ['lmerMod']
Formula: SNB_new ~ Condition_SNB * Dotcounting_SNB + (1 | subject_SNB)
Data: na.omit(s_complete)
REML criterion at convergence: 628.8
Scaled residuals:
Min 1Q Median 3Q Max
-3.3291 -0.8359 -0.2435 0.7277 2.9290
Random effects:
Groups Name Variance Std.Dev.
subject_SNB (Intercept) 0.2609 0.5107
Residual 1.4618 1.2081
Number of obs: 3506, groups: subject_SNB, 209
Fixed effects:
Estimate Std. Error
library(lme4)
model <- lmer(FOD_FODs ~ Condition_FODs + (1|subject_FODs), data = indvar_FODs)
The output shows a main effect of condition F(4,3139)=0.347, p=0.846 and a main effect of age F(15,3139)=0.84, p=1.000, no significant interaction effect F(45,3139)=0.046, p=1.000. The partial eta squared values were 0.0004 for condition, 0.0004 for age and 0.0006 for the interaction term. The marginal and conditional R squared values were 0.0012 and 0.0011 respectively.
#No-test case:
#Test case:
mod1 <- glm(categorical_FOD_FODs ~ 0 + Condition_FODs:Language_used_FODs, data = indvar_FODs, family = "binomial")
mod1
summary(mod1)
#If the code for my condition is 00 and the code for my language_used is 01, how can I calculate the probability of my categorical_FOD_FODs given my condition and language_used?
#If I am interested in the odds ratio (OR), how can I calculate it?
fix invalid codeWed, 10 May 2023 I did not check the code by running, try to check it by running.
library(multcompView)
posthoc <- glht(anova_result, linfct = mcp(Condition_SNB = "Tukey"))
library(multcomp)
posthoc <- glht(anova_result, linfct = mcp(Condition_SNB = "Tukey"))
summary(posthoc)
we can add the L1(0.001) penalty to the glmer function in the R package penalized. The penalty is added to the log-likelihood function and is minimized along with the rest of the function
indvar_FODs$categorical_FOD_FODs[which(!is.na(indvar_FODs$categorical_FOD_FODs))]
indvar_FODs$ID[which(!is.na(indvar_FODs$categorical_FOD_FODs))]
indvar_FODs$Gender_quantised_FODs[which(!is.na(indvar_FODs$categorical_FOD_FODs))]
> pairwise
contrast estimate SE df t.ratio p.value
A - B 0.0747 0.190 3317 0.392 0.9795
A - C -0.0476 0.190 3317 -0.250 0.9945
A - D 0.0215 0.191 3317 0.113 0.9995
B - C -0.1222 0.190 3317 -0.643 0.9182
B - D -0.0532 0.191 3317 -0.279 0.9924
C - D 0.0691 0.190 3317 0.363 0.9836
diamonds %>%
ggplot() +
aes(x = price, y = cut) +
geom_bin2d()
model_FODs<-multinom(categorical_FOD_FODs~Condition_FODs*Language_used_FODs, data=indvar_FODs)
> summary(mod)
Generalized linear mixed model fit by the Laplace approximation
Formula: categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + (1 | subject_FODs)
Data: indvar_FODs
AIC BIC logLik deviance
1805.1 1828.0 -898.6 1797.1
Random effects:
Groups Name Variance Std.Dev.
subject_FODs (Intercept) 0.000981 0.03110
Number of obs: 3360, groups: subject_FODs, 240
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.3914 0.1146 -3.417 0.000633 ***
Condition_FODsB 0.2631 0.1572 1.674 0.094595 .
Condition_FODsC 0.1414 0.1572 0
fix invalid codeMon, 01 May 2023 Model_DD <- glmer(cat_SNB_new ~ Condition_SNB * Gender_SNB + (1 + Condition_SNB | subject_SNB),
data = s_complete,
family = binomial,
control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 20000)))
fix invalid codeWed, 24 May 2023 ### Results
In terms of accuracy, the model predicted the result for each subject correctly 84.6% of the time. This is slightly above the chance level of 80%.
m2 <- clm(categorical_FODs ~ 0 + Condition_FODs:Language_used_FODs - 1, random = ~1|subject, data = indvar_FODs)
summary(m2)
> mod <- glmer('categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + (1 | subject_FODs)',
+ data = indvar_FODs, family = binomial)
> summary(mod)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: binomial ( logit )
Formula: categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + (1 |
subject_FODs)
Data: indvar_FODs
AIC BIC logLik deviance df.resid
8289.4 8332.8 -4139.7 8279.4 4571
Scaled residuals:
Min 1Q Median 3Q Max
-2.5844 -0.7938 -0.1721 0.7055 28.6594
Random effects:
Groups Name
library(lme4)
# data
data <- read.table("http://www.stat.ufl.edu/~aa/social/data/nlswork.dat", sep="",
header=T, na.strings=".")
# logistic regression model
model <- glmer(union ~ educ + (1|id), data=data, family="binomial")
summary(model)
library(multcomp)
#comparison of all three conditions
TukeyHSD(mod_FODs, "x1")
#comparison of languages
TukeyHSD(mod_FODs, "x2")
#comparison for conditions and languages
TukeyHSD(mod_FODs, "x1:x2")
Edit this according to the data
#get odds ratios
exp(cbind(OR = coef(mod_FODs), confint(mod_FODs)))
#get odds ratios with robust standard errors
exp(cbind(OR = coef(mod_FODs), confint(mod_FODs, method = "Wald")))
fix invalid codeFri, 05 May 2023 > Model_MSG <- glmer(cat_SNB_new ~ 0 + Condition_SNB * Dotcounting_SNB - 1 + (1 | subject_SNB),
+ data = s_complete,
+ family = binomial,
+ control = glmerControl(optimizer = "bobyqa", tolPwrss = 1e-08))
Warning message:
In commonArgs(par, fn, control, environment()) :
maxfun < 10 * length(par)^2 is not recommended.
mod <- glmer('categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + (1 | subject_FODs)',
data = indvar_FODs, family = binomial)
avfourfods$theme$colors$fill
fix invalid codeFri, 12 May 2023 ### Compute the number of errors and number of error-free sentences
fix invalid codeMon, 01 May 2023 {r}
summary(Model_DD_1)
Model_MSG <- glmer(cat_SNB_new ~ 0 + Condition_SNB * Dotcounting_SNB - 1 + (1 | subject_SNB),
data = s_complete,
family = binomial,
control = glmerControl(optimizer = "bobyqa", tolPwrss = 1e-08))
chisq.test(cont_table)
# Perform chi-squared test
chisq_result <- chisq.test(cont_table)
chisq_result
Pearson's Chi-squared test
data: cont_table
X-squared = 105.87, df = 6, p-value < 2.2e-16
> glm(categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + Condition_FODs:Language_used_FODs, data=indvar_FODs, family=multinomial(link="logit"))
Warning message:
glm.fit: fitted probabilities numerically 0 or 1 occurred
Call: glm(formula = categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + Condition_FODs:Language_used_FODs, family = multinomial(link = "logit"), data = indvar_FODs)
Coefficients:
(Intercept)
ascending -0.002
descending 1.511
identity 3.368
Condition_FODsB 0.079
Condition_FODsC 0.243
Condition_FODsD
fix invalid codeWed, 24 May 2023 # eta_squared(model)
fix invalid codeTue, 18 Apr 2023 ``
### Use `#` to create an H2 header:
This will show up in the table of contents as well.
### Use `##` to create an H3 header:
This will show up in the table of contents as well.
This will show up in the table of contents as well.
### Use `###` to create an H4 header:
This will show up in the table of contents as well.
### Use `####` to create an H5 header:
This will show up in the table of contents as well.
### Use `#####` to create an H6 header:
This will show up in the table of contents as well.
### Use `*` or `-` or `+` to create a bulleted list:
- Item 1
- Item 2
- Item 2a
- Item 2b
- Item 2c
- Item 3
- Item 3a
- Item 3b
- Item 4
### Use numbers followed by periods to create a numbered list:
1. Item 1
1. Item
fix invalid codeSat, 22 Apr 2023 ### oops
i did a bad thing:
The function add takes two arguments a and b and returns the sum of a and b.
The Chi-square test for independence is used to determine if there is a significant association between two categorical variables. Results show a P-value of 0.8828 for Pearson's Chi-squared test of independence. As the P-value is greater that the threshold of 0.05, we fail to reject the null hypothesis and conclude that there is no statistically significant association between categorical_FOD_FODs and Gender_FODs.
library(car)
library(nnet)
library(MASS)
set.seed(123)
n = length(indvar_FODs$Condition_FODs)
FODs.train <- sample(1:n, 0.7*n)
FODs.test <- setdiff(1:n, FODs.train)
train.features <- model.matrix(categorical_FOD_FODs ~ -1 + Condition_FODs * Language_used_FODs, data=indvar_FODs)
test.features <- model.matrix(categorical_FOD_FODs ~ -1 + Condition_FODs * Language_used_FODs, data=indvar_FODs)[FODs.test,]
FODs.pred.1 <- multinom(categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs, data=indvar_FODs[FODs.train,])
FODs.pred.
aggregate(s_complete$num_SNB~s_complete$Condition_SNB, FUN = function(x){
c(mean = mean(x), sd = sd(x), count = length(x))
})
mod <- glmer('categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + (1 | subject_FODs)', data = indvar_FODs,
family=binomial)
## Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
## Family: binomial ( logit )
## Formula: categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs +
## (1 | subject_FODs)
## Data: indvar_FODs
##
## AIC BIC logLik deviance df.resid
## 4616.503 4671.476 -2299.252 4598.503 3312
##
## Random effects:
## Groups Name Std.Dev.
## subject_FODs (Intercept) 1.403e-07
## Number of obs: 3321, groups: subject_
Pairwise t-test
$`A - B`
Welch Two Sample t-test
data: mean_age by group
t = 0.3924, df = 3316.7, p-value = 0.97952
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.3199497 0.4645497
sample estimates:
mean in group A mean in group B
53.065 52.991
$`A - C`
Welch Two Sample t-test
data: mean_age by group
t = -0.25048, df = 3316.7, p-value = 0.9945
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.3165497 0.2445497
sample estimates:
mean in group A mean in group C
53.065 53.123
$`A - D`
Welch Two Sample t-test
data
Dear Dr. Kanai,
The results of the lme model are reported in APA style in the table below:
| Effect | F-value | p-value | eta-squared |
|--------|---------|---------|-------------|
| Condition | 0.347 | 0.846 | 0.004 |
| Age | 0.084 | 1.000 | 0.004 |
| Interaction | 0.046 | 1.000 | 0.007 |
Thank you for your time,
Best wishes,
[your name]
## Task 8. (1 point)
Why do we use the command `summary(model1)` to obtain the estimates of the fixed effects?
Your answer should be given in the following format:
model <- multinom(categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + Condition_FODs*Language_used_FODs, indvar_FODs)
summary(model)
fix invalid codeSat, 06 May 2023 1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, ... :
Model failed to converge with max|grad| = 0.000191281 (tol = 0.001, component 1)
2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, ... :
Model failed to converge with max|grad| = 0.000372457 (tol = 0.001, component 1)
3: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, ... :
Model is nearly unidentifiable: very large eigenvalue
- Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
- Rescale variables?
> model_full2 <- glmer(cat_SNB_new ~ 0 + developmentaldisorder_SNB * Condition_SNB -1 + (1|subject_SNB),
+ data =
mod <- glmer('categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + (1 | subject_FODs)',
data = indvar_FODs, family = binomial)
fix invalid codeSat, 06 May 2023 model_full2 <- glmer(cat_SNB_new ~ 0 + developmentaldisorder_SNB * Condition_SNB -1 + (1+developmentaldisorder_SNB|subject_SNB),
data = s_complete, family = binomial)
fix invalid codeSat, 06 May 2023
# First, we would like to look at the factors that influence the individual variable choice
# i.e. what variable was chosen in what condition for a given participant
indvar_FODs <- FOD_FODs %>%
group_by(subject_FODs, Condition_FODs, Language_used_FODs, categorical_FOD_FODs) %>%
summarise(n = n())
indvar_FODs$subject_FODs <- as.factor(indvar_FODs$subject_FODs)
indvar_FODs$Condition_FODs <- as.factor(indvar_FODs$Condition_FODs)
indvar_FODs$Language_used_FODs <- as.factor(indvar_FODs$Language_used_FODs)
indvar_FODs$categorical_FOD_FODs <- as.factor(indvar_FODs$categorical_FOD_FODs)
# we need to make
cont_table
chisq.test(cont_table)
fix invalid codeSat, 06 May 2023 > anova_table <- Anova(model, type = "III")
> library(report)
> report(anova_table)
Error in `[[<-.data.frame`(`*tmp*`, "Mean_Square", value = numeric(0)) :
Ersetzung hat 0 Zeilen, Daten haben 3
multinom(categorical_FOD_FODs ~ Condition_FODs * Language_used_FODs, indvar_FODs)
avfourfods +
scale_fill_manual(values = c("#00AFBB", "#E7B800", "#FC4E07", "#56B1F7"))
> summary(s_complete$num_SNB)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.000 0.000 0.000 3.854 6.000 12.000
> aggregate(num_SNB ~ Condition_SNB, data = s_complete, FUN = summary)
Condition_SNB num_SNB.Min. num_SNB.1st Qu. num_SNB.Median num_SNB.Mean
1 1 0.0 0.5 2.5 4.5
2 2 0.0 0.0 0.0 2.0
3 3 0.0 0.0 0.0 2.0
num_SNB.3rd Qu. num_SNB.Max.
1 6.0 12
2 0.0 4
3 0.0 6
> sd(s_complete$num_SNB)
[1] 4.553355
>
F(4, 3139) = 0.34, p = 0.846, η2 = 0.000442.
F(15, 3139) = 0.08, p = 1.000, η2 = 0.000402.
F(45, 3139) = 0.05, p = 1.000, η2 = 0.000652.
R2m = 0.00117
R2c = 0.00118
fix invalid codeWed, 24 May 2023 {r}
# condition_SNB
# language_used_SNB
# eta_squared(model)
# fix_SNB
# condition_SNB * fix_SNB
# eta_squared(model)
# condition_SNB
# language_used_SNB
# eta_squared(model)
The function add takes two arguments a and b and returns the sum of a and b.
model1 <- lme(FOD_FODs ~ 0 + Condition_FODs*language_used_FODs_ind, random=~1|subject_FODs,
na.action=na.omit,
method = "ML", data=indvar_FODs,control = lmeControl(opt = "optim"))
p_value <- shapiro.test(resid(model1))$p.value
round(p_value, 3)
[1] 0
check_singularity(model1)
[1] FALSE
ano <- anova(model1)
round(ano, 3)
numDF denDF F-value p-value
Condition_FODs 4 3262 0.402 0.807
language_used_FODs_ind 4 40 0.008 1.000
Condition_FODs:language_used_FODs_ind 12 3262 0.
## 4.3.3.3
> fit3 <- glm(categorical_FOD_FODs ~ 0 + Condition_FODs + Gender_FODs - 1,
+ family = binomial(link = 'logit'), data = indvar_FODs)
> summary(fit3)
Call:
glm(formula = categorical_FOD_FODs ~ 0 + Condition_FODs + Gender_FODs -
1, family = binomial(link = "logit"), data = indvar_FODs)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.315 -1.159 -1.111 1.191 1.268
Coefficients:
Estimate Std. Error z value Pr(>|z|)
Condition_FODsA 0.06420 0.09171 0.699 0.484
Condition_FODsB
library(nnet)
#without interaction effect
#m1 <- multinom(categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs, data=indvar_FODs)
#summary(m1)
#with interaction effect
m2 <- multinom(categorical_FOD_FODs ~ Condition_FODs * Language_used_FODs, data=indvar_FODs)
summary(m2)
#in this case, the interaction effect is not significant
#the softmax regression with interaction effect is worse than that without interaction effect
#cross validation with bootstrap
library(boot)
#train the model with all data
#without interaction effect
model1 <- multinom(categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs, data=indvar_FODs)
#with interaction effect
model2 <- multinom(categorical_FOD_FODs ~ Condition_FODs * Language_used_
mod <- glmer('categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + (1 | subject_FODs)',
+ data = indvar_FODs, family = binomial)
ranef(mod)
fixef(mod)
# Use anova to compare the models
anova(mod, mod_new)
### Contingency table
cont_table <- table(iris$Species, iris$Petal.Width)
cont_table
### Perform chi-squared test
chisq_result <- chisq.test(cont_table)
chisq_result
### Post-hoc test
# step 1: install the package
install.packages("vcd")
# step2: load the package
library(vcd)
# step3: test
assocstats(cont_table)
library(openxlsx)
read.xlsx2('C:/Users/Nathan/Desktop/indvar_FODs.xlsx', sheet = 1, as.data.frame = T)
model1 <- lme(FOD_FODs ~ 0 + Condition_FODs * developmentaldisorder_FODs, random = ~1|subject_FODs, na.action = na.omit, method = "ML", data = indvar_FODs, control = lmeControl(opt = "optim")) anova(model1) numDF denDF F-value p-value Condition_FODs 4 3271 0.4038075 0.8060 developmentaldisorder_FODs 1 43 0.3213196 0.5738 Condition_FODs:developmentaldisorder_FODs 3 3271 0.0147161 0.9976 shapiro.test(resid(model1)) Shapiro-Wilk normality test data: resid(model1) W = 0.98738, p-value < 2.2e-16
fix invalid codeFri, 12 May 2023 * 运行代码
fix invalid codeSat, 06 May 2023 ### level 2
fix invalid codeFri, 12 May 2023
fix invalid codeFri, 12 May 2023 ### This is how you insert image
### This is how you insert video
<iframe width="420" height="315" src="https://www.youtube.com/embed/O72FWNeO-xY" frameborder="0" allowfullscreen></iframe>
### How to add a hyperlink
<a href="https://www.rstudio.com/">Visit RStudio </a>
### How to insert a footnote
Please check the footnote[^footnote].
[^footnote]: Please note that this is the footnote.
### How to highlight some text
** Important note**: Please check the readme file before proceeding.
### How to add a table
| Column 1 | Column 2 |
|----------|----------|
| A1 | B1 |
| A2 | B2 |
| A3 | B3 |
### How to add a table with column
fix invalid codeMon, 01 May 2023 Model_DD <- glmer(cat_SNB_new ~ Condition_SNB * Gender_SNB + (1 | subject_SNB) + (1 | session_SNB),
data = s_complete,
family = binomial,
control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 20000)))
head(indvar_FODs)
indvar_FODs$Condition_FODs <- relevel(indvar_FODs$Condition_FODs, ref = "Not")
indvar_FODs$Language_used_FODs <- relevel(indvar_FODs$Language_used_FODs, ref = "English")
m1 <- clm(categorical_FOD_FODs ~ 0 + Condition_FODs:Language_used_FODs, data = indvar_FODs)
summary(m1)
m1 <- clm(categorical_FOD_FODs ~ 0 + Condition_FODs:Language_used_FODs + I(Age_FODs^2), data = indvar_FODs)
summary(m1)
fix invalid codeWed, 24 May 2023 ### Use \*\* to create superscripts and subscripts
I've found that using `^` and `_` to create superscripts and subscripts can make the code hard to read. I like to use `\*\*` instead.
set.seed(123)
FODs <- with(indvar_FODs,
data.frame(
categorical_FOD_FODs = categorical_FOD_FODs,
Condition_FODs = factor(Condition_FODs),
Language_used_FODs = factor(Language_used_FODs),
stringsAsFactors = FALSE
)
)
n <- nrow(FODs)
train <- sample(1:n, n * 0.8)
chisq.test(cont_table)
logit<-glm(categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs,
data = indvar_FODs,
family = binomial(link = "logit"))
summary(logit)
#install.packages("multcomp")
library("multcomp")
categorical_FOD_FODs_analysis <- glht(model_categorical_FOD_FODs, linfct = mcp(Condition_FODs = "Tukey", Gender_FODs = "Tukey"))
summary(categorical_FOD_FODs_analysis)
install.packages("xtable")
library(xtable)
print(xtable(mytable, digits = 2),type="latex",file="mytable.tex",caption.placement = "top")
fix invalid codeMon, 01 May 2023
fix invalid codeSat, 06 May 2023 r
x <- matrix(1:3, ncol=3, nrow=3)
x
anova(lm(categorical_FOD_FODs ~ Condition_FODs*Gender_FODs, data=indvar_FODs))
[R code for multinomial logistic regression](https://github.com/zgols/open-science-course/blob/master/code/code_multinomial_logistic_regression.R)
[R code for multinomial logistic regression - FODs](https://github.com/zgols/open-science-course/blob/master/code/code_multinomial_logistic_regression_FODs.R)
[R code for multinomial logistic regression - FODs - hierarchical](https://github.com/zgols/open-science-course/blob/master/code/code_multinomial_logistic_regression_FODs_hierarchical.R)
[R code for multinomial logistic regression - FODs - hierarchical - intercepts](https://github.com/zgols/open-science-course/blob/master/code/code_multinomial_logistic_regression_FODs_hierarchical_intercepts.R)
fit3 <- glm(categorical_FOD_FODs ~ 0 + Condition_FODs + Language_used_FODs + Condition_FODs:Language_used_FODs - 1, family = binomial(link = 'logit'), data = indvar_FODs)
summary(fit3)
library(multcomp)
glht1 <- glht(model, linfct = mcp(Language_used_FODs="Tukey"))
summary(glht1)
confint(glht1)
fix invalid codeWed, 24 May 2023 We can see that this does what we expect, but it's not perfect. Let's try
tweaking the code a bit:
fix invalid codeFri, 12 May 2023 > s_complete[s_complete$Dotcounting_SNB == "missing", names(s_complete) == "Dotcounting_SNB", drop = FALSE]
m1 <- glm(categorical_FOD_FODs ~ Condition_FODs*Language_used_FODs, data = indvar_FODs, family = "multinomial")
summary(m1)
function that: replace 0 with NA when there are missing data points
library(knitr)
kable(fit4)
fix invalid codeWed, 24 May 2023 Anova Table (Type II tests)
Response: Sepal.Length
Sum Sq Df F value Pr(>F)
Petal.Width 0.448 1 0.06633 0.7974
Residuals 233.924 98
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##The analysis is a multilevel model with a categorical dependent variable and a categorical independent variable. This is a multilevel generalized linear model, so glmer must be used.
model <- glmer(categorical_FOD_FODs ~ Condition_FODs*Gender_quantised_FODs + (1|subject_FODs), data = indvar_FODs, family = 'binomial')
summary(model)
The function add takes two arguments a and b and returns the sum of a and b.
#randomised factorial design
#https://stats.stackexchange.com/questions/83585/two-way-anova-with-random-factors
#https://stats.stackexchange.com/questions/7384/how-to-specify-random-effects-in-a-mixed-model-using-the-lmer-function?rq=1
#https://www.thertrader.com/2016/01/30/repeated-measures-anova-and-linear-mixed-effects-models-in-r/
#https://stats.stackexchange.com/questions/255888/how-to-compute-two-way-interaction-effect-in-lmer
#https://rstudio-pubs-static.s3.amazonaws.com/200194_fc1c5f5ab7a34c99b9d7d0c0faa936e1.html
#https://stackoverflow.com/questions/43821727/how-to-test-
fix invalid codeWed, 24 May 2023 eta_squared <- eta_squared(model, measure = "RMSEA", type = "robust")
fix invalid codeWed, 24 May 2023 ##### Effect sizes
Hey, I thought we talked about this.
fix invalid codeWed, 24 May 2023 > library(GGally)
> eta_squared <- etaSquared(model)
> print(eta_squared)
{r}
m2 <- clm(categorical_FOD_FODs ~ 0 + Condition_FODs:Gender_quantised_FODs - 1,
random = ~1|subject_FODs/Gender_quantised_FODs, data = indvar_FODs)
summary(m2)
fix invalid codeSat, 06 May 2023 model_full2 <- glmer(cat_SNB_new ~ 0 + developmentaldisorder_SNB * Condition_SNB -1 + (1+developmentaldisorder_SNB|subject_SNB),
data = s_complete, family = binomial)
m1 <- glm(categorical_FOD_FODs ~ Condition_FODs*Language_used_FODs, data = indvar_FODs, family = "multinomial")
Error in get(family, mode = "function", envir = parent.frame()) :
object 'multinomial' of mode 'function' was not found
> glmer(categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + Condition_FODs * Language_used_FODs + (1|subject_FODs), data = indvar_FODs, family = binomial("logit"))
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation)
['glmerMod']
Family: binomial ( logit )
Formula: categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + Condition_FODs * Language_used_FODs + (1 | subject_FODs)
Data: indvar_FODs
AIC BIC logLik deviance df.resid
17063.6 17159.1 -8523.8 17047.6 8376
Scaled residuals:
Min 1Q Median 3Q Max
-3.3123 -0.6710 0.14
fix invalid codeMon, 24 Apr 2023
fix invalid codeWed, 24 May 2023 And the result will be:
### Example 2: multiple plots in an Rmarkdown document
glmer(cat_SNB_new ~ 0 + Condition_SNB * Dotcounting_SNB - 1 + (1 | subject_SNB),
data = s_complete,
family = binomial,
control = glmerControl(optimizer = "bobyqa", tolPwrss = 1e-08))
q <- glht(m1, linfct = mcp(Condition_FODs:Language_used_FODs = "Tukey"))
q
plot(q)
fix invalid codeFri, 05 May 2023
fix invalid codeFri, 05 May 2023 ## Basic usage
### To check whether your R code is correct, use `lintr::lint()`
**If your R code is fine, this should return `NULL`:**
library(multcomp)
library(lsmeans)
posthoc <- glht(anova_result, linfct = mcp(Condition_SNB = "Tukey"))
summary(posthoc)
fix invalid codeSat, 06 May 2023 model_full2 <- glmer(cat_SNB_new ~ 0 + developmentaldisorder_SNB * Condition_SNB -1 + (1+developmentaldisorder_SNB|subject_SNB),
data = s_complete, family = binomial)
fix invalid codeWed, 24 May 2023 #### FIXME: <a id="fixme"></a>
fix invalid codeFri, 05 May 2023 - the html output should be like this:
- to cite the analysis code:
This model is a linear mixed model with maximum likelihood estimation, here we use the glmer() function from the lme4 package to fit this model. We use the default setting of the glmer() function, here the optimizer is the L-BFGS-B algorithm. In addition, we set the maximum iteration number to be 6000000 to ensure the convergence of the model.
The function add takes two arguments a and b and returns the sum of a and b.
fix invalid codeFri, 12 May 2023
fix invalid codeSat, 06 May 2023 model_full2 <- glmer(cat_SNB_new ~ 0 + developmentaldisorder_SNB * Condition_SNB -1 + (1+developmentaldisorder_SNB|subject_SNB),
data = s_complete, family = binomial)
-------------------------
**Annotated code**
-------------------------
model1 <- glmer(categorical_FOD_FODs~Condition_FODs*Gender_FODs,
random=~1|subject_FODs,
na.action=na.omit,
method = "ML",
data=indvar_FODs,
control = lmeControl(opt = "optim"))
glmer(categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + (1 | subject_FODs),
data = indvar_FODs, family = binomial, method = "bobyqa")
ggplot2:::ggplot_build.ggplot
mod_FODs_optimizer <- glmer(y ~ x1 * x2 + (1|subject_FODs), family = binomial, data = indvar_FODs,
control = glmerControl(optimizer = "bobyqa"))
f <- function(x) {
return(x[1]^2 + x[2]^2)
}
optim(c(0,0), f) #returns: $par
#[1] 7.105427e-15 7.105427e-15
#$value
#[1] 1.110223e-16
#$counts
#function gradient
#$convergence
#[1] 0
#$message
#NULL
p <- summary(fit4)$coefficients[,4]
r.squaredGLMM(fit4)
anova(fit4, test = "Chisq")
anova(fit4, test = "LRT")
# extract the p values, then put them in a data frame
df.p <- data.frame(
p = p
)
# extract the coeffients, then put them in a data frame
df.coef <- data.frame(
coef = c(summary(fit4)$coefficients[,1])
)
# extract the coeffients, then put them in a data frame
df.se <- data.frame(
se = c(summary(fit4)$coefficients[,2])
)
# extract the coeffients, then put them in a data frame
df.z <- data.frame(
z = c(summary(fit4)$coefficients[,3])
)
# extract the names and put them in a vector
df.names <- rownames(summary(fit4)$coefficients
x <- rnorm(10)
# optimized glmer
mod_FODs <- glmer(y ~ x1 * x2 + (1|subject_FODs), family = binomial, data = indvar_FODs,
control = glmerControl(optimizer = "bobyqa"))
# print model
summary(mod_FODs)
#> boundary (singular) fit: see help('isSingular')
#> Generalized linear mixed model fit by maximum likelihood
#> (Laplace Approximation) [glmerMod]
#> Family: binomial ( logit )
#> Formula: y ~ x1 * x2 + (1 | Subject_FODs)
#> Data: indvar_FODs
#>
#> AIC BIC logLik deviance df.resid
#> 48.717 56.704 -17.358 34.717 20
#>
#> Scaled residuals:
#> Min 1Q
Sample Answer: There was no significant interaction between condition and age, F(45,3139)=0.046, p=1.000. There was no main effect of condition, F(4,3139)=0.347, p=0.846, ηp²=.004, nor was there a main effect of age, F(15,3139)=0.084, p=1.000, ηp²=.004. The adjusted R² for the model was .001.
> summary(mod)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: binomial ( logit )
Formula: categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + (1 | subject_FODs)
Data: indvar_FODs
AIC BIC logLik deviance df.resid
4616.5 4671.5 -2299.3 4598.5 3312
Scaled residuals:
Min 1Q Median 3Q Max
-2.8934 -0.8258 -0.2503 0.6317 6.8485
Random effects:
Groups Name Std.Dev.
subject_FODs (Intercept) 1.403e-07
Number of obs: 3321, groups: subject_FODs, 44
Fixed effects:
Estimate Std. Error z value Pr(>|
1. The main effects for condition and age are not significant at the p < .05 level. The interaction between age and condition is not significant at the p < .05 level.
2. The effect size for condition is .004, the effect size for age is .004 and the effect size for the interaction between age and condition is .01.
# Perform chi-squared test
chisq_result <- chisq.test(cont_table)
chisq_result
Pearson's Chi-squared test
data: cont_table
X-squared = 105.87, df = 6, p-value < 2.2e-16
### Code analysis:
#### Code:
def add(a, b):
return a + b
#### Explanation:
The function add takes two arguments a and b and returns the sum of a and b.
## Findings:
* ...
* ...
* ...
## Discussion:
* ...
* ...
* ...
## Conclusion:
* ...
* ...
* ...
## References:
* ...
* ...
* ...
## Appendix:
* ...
* ...
* ...
## Data Archiving statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.
## Authors contributions:
* ...
* ...
* ...
## Acknowledgements:
* ...
* ...
* ...
## Example 2: Assign a value
fix invalid codeSat, 06 May 2023 model_full2 <- glmer(cat_SNB_new ~ 0 + developmentaldisorder_SNB * Condition_SNB -1 + (1+developmentaldisorder_SNB|subject_SNB),
data = s_complete, family = binomial)
numDF denDF F-value p-value
Condition_FODs 4 3271 0.4038075 0.8060
developmentaldisorder_FODs 1 43 0.3213196 0.5738
Condition_FODs:developmentaldisorder_FODs 3 3271 0.0147161 0.9976
> shapiro.test(resid(model1)) Shapiro-Wilk normality test data: resid(model1) W = 0.98738, p-value < 2.2e-16
mod <- glmer('categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + (1 | subject_FODs)',
data = indvar_FODs, family = binomial)
mod2 <- glmer('categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs + (1 | subject_FODs)',
data = indvar_FODs, family = binomial,
control = glmerControl(optimizer='Nelder_Mead', optCtrl=list(maxfun=1e5)))
anova(mod, mod2)
> logit_FODs<-glm(categorical_FOD_FODs~Condition_FODs + Language_used_FODs, data=indvar_FODs, family=multinomial(logit))
> summary(logit_FODs)
Call:
glm(formula = categorical_FOD_FODs ~ Condition_FODs + Language_used_FODs,
family = multinomial(logit), data = indvar_FODs)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.57216 -0.82433 -0.57809 -0.40977 3.87957
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.72358 0.05834 -12.422 < 2e-16 ***
Condition_FODsB 0.15697
fix invalid codeFri, 05 May 2023 > s_complete$developmentaldisorder_SNB <- factor(s_complete$developmentaldisorder_SNB)
> table(s_complete$developmentaldisorder_SNB)
No Yes
2712 796
> levels(s_complete$developmentaldisorder_SNB)
[1] "No" "Yes"
{r}
#print(table(indvar_FODs$Gender_FODs, indvar_FODs$categorical_FOD_FODs))
#chisq.test(table(indvar_FODs$Gender_FODs, indvar_FODs$categorical_FOD_FODs))
# Perform chi-squared test
chisq_result <- chisq.test(cont_table)
# Print output of chi-squared test
chisq_result
fix invalid codeMon, 01 May 2023 Model_DD <- glmer(cat_SNB_new ~ Condition_SNB * Gender_SNB + (1 + Condition_SNB | subject_SNB),
data = s_complete,
family = binomial,
control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 20000)))
fix invalid codeFri, 12 May 2023 {r}
indvar_SNB[is.na(indvar_SNB$Dotcounting_SNB),]
clm1 <- clm(categorical_FOD_FODs ~ 0 + Condition_FODs:Language_used_FODs, data = indvar_FODs)
A linear mixed-effects model was fitted with the FOD score as dependent variable, the condition (condition_FODs) and the developmentaldisorder (developmentaldisorder_FODs) as predictors, the subject (subject_FODs) as random effect, and the number of trials as control parameter. The residuals of the model were tested for normality using the Shapiro-Wilk test, which indicated that the model residuals were normally distributed (W = 0.98738, p < 2.2e-16). The model showed no main effects (p > 0.05).
{r}
explanation of the code:
model <- lmer(num_SNB ~ 0 + Condition_SNB*Language_used_SNB + (1|subject_SNB), data = s_complete)
anova(model)
shapiro.test(resid(model))
avfourfods$labels$fill
fix invalid codeMon, 01 May 2023 > levels(s_complete$Gender_SNB)
[1] "Female" "Male" "Non-Binary / Genderqueer"
> anova(Model_DD_1)
Analysis of Variance Table
npar Sum Sq Mean Sq F value
Condition_SNB 4 136.301 34.075 34.0754
Gender_SNB 2 1.684 0.842 0.8421
Condition_SNB:Gender_SNB 6 7.063 1.177 1.1771
#Load data
data("FOD_data")
#Preprocess data
FOD_data<-FOD_data[,-c(1,3,5)]
#Split data
set.seed(1)
indvar_FODs<-FOD_data[sample(nrow(FOD_data), size = 0.7*nrow(FOD_data), replace = FALSE),]
outvar_FODs<-FOD_data[-which(rownames(FOD_data) %in% rownames(indvar_FODs)),]
#Partial effect plot
#Generalized linear model
FODs_glm<-glm(categorical_FOD_FODs~Condition_FODs*Language_used_FODs, data = indvar_FODs, family = "binomial")
summary(FODs_glm)
#Multinomial logistic regression
FODs_multinomial<-multinom(categorical_FOD_FODs~Condition_FOD
fix invalid codeFri, 21 Apr 2023 \begin{table}[H]
\centering
\begin{tabular}{rrrrr}
\hline
& Estimate & Std. Error & z value & Pr($>$$|$z$|$) \\
\hline
Condition\_FODsA:nativelang\_FODsDanish & -0.32 & 0.46 & -0.69 & 0.49 \\
Condition\_FODsB:nativelang\_FODsDanish & 0.11 & 0.46 & 0.23 & 0.82 \\
Condition\_FODsC:nativelang\_FODsDanish & 0.32 & 0.46 & 0.69 & 0.49 \\
Condition\_FODsD:nativelang\_FODsDanish & 0.11 & 0.46 & 0.23 & 0.82 \\
Condition\_FODsA:nativelang\_FODsEnglish & -0.54 & 0.48 & -1.13 & 0.26 \\
anova_result <- aov(SNB ~ Condition, data = csv_data)
posthoc <- glht(anova_result, linfct = mcp(Condition_SNB = "Tukey"))
df_posthoc <- tidy(posthoc)
df_posthoc$`Tukey HSD` <- round(df_posthoc$`Tukey HSD`, digits = 2)
df_posthoc
fix invalid codeTue, 18 Apr 2023 </details>
### Exercise 5
The following data contains the results of a competition where participants had to perform a task. The dataframe consists of the number of times participants completed the task and the amount of time they spent performing the task.
<details>
<summary>Hint 1</summary>
<details>
<summary>Hint 2</summary>
<details>
<summary>Hint 3</summary>
<details>
<summary>Hint 4</summary>
<details>
<summary>Hint 5</summary>
<details>
<summary>Hint 6</summary>
<details>
<summary>Hint 7</summary>
<details>
<summary>Hint 8</summary>
<details>
<summary>Hint 9</summary>
<details>
<summary>Hint 10</summary>
<details>
<summary>Hint 11</summary>
<details>
<summary>Hint 12</summary>
<
fix invalid codeTue, 18 Apr 2023 > str(mydata$Language_used_FODs)
chr [1:3520] "Turkish" "Turkish" "Turkish" "Turkish" "Turkish" "Turkish" "Turkish" "Turkish" "Turkish" "Turkish" "Turkish" ...
> mydata$Language_used_FODs <- factor(mydata$Language_used_FODs)
> library(multcomp)
> glht1 <- glht(model, linfct = mcp(Language_used_FODs="Tukey"))
Error in mcp2matrix(model, linfct = linfct) :
Variable(s) ‘Language_used_FODs’ of class ‘character’ is/are not contained as a factor in ‘model’.
this did not solve the problem. please solve this problem
fix invalid codeFri, 05 May 2023 > Model_MSG <- glmer(cat_SNB_new ~ 0 + Condition_SNB * Dotcounting_SNB - 1 + (1 | subject_SNB),
+ data = s_complete,
+ family = binomial,
+ control = glmerControl(optimizer = "bobyqa", maxfun = 5000))
Warning message:
In commonArgs(par, fn, control, environment()) :
maxfun < 10 * length(par)^2 is not recommended.
fix invalid codeMon, 24 Apr 2023 av_data1_all<-s_complete %>%
group_by(Gender_SNB, num_SNB)%>%
dplyr::summarize(prop=mean(num_SNB), synch = sum(num_SNB),
n = n(), .groups='drop')
# This is formatted as code
# variables that are categorical
class(data$variable)
# variables that are continuous
class(data$variable)
# check for missing values
colSums(is.na(data))
# check for outliers
boxplot(data$variable)
# check for normality
qqnorm(data$variable)
qqline(data$variable)
# convert categorical variable to a factor variable
data$variable <- factor(data$variable)
# linear regression
model <- lm(outcome ~ predictor)
# ANOVA
model <- aov(outcome ~ predictor)
# T-test
model <- t.test(outcome ~ predictor)
# Chi-square test
model <- chisq.test(outcome ~ predictor)
# simple correlations
cor(data$variable1, data$variable2)
# repeated measures ANOVA
model <- aov(outcome ~ predictor + Error(subject/predictor), data = data)
# paired T-test
model <- t.test(variable1 ~ variable2)
#
## Screenshots
Here is a screenshot of the output of the code:
# Questions?
If you have general questions about the project or the data, please contact the researchers who conducted the study.
If you have general questions about this repository, please contact Cathy Tian (cathytian@hku.hk).
fix invalid codeFri, 05 May 2023 2. fix invalid code:
fix invalid codeMon, 24 Apr 2023 {r av_data1_all}
av_data1_all<-s_complete %>%
group_by(Gender_SNB, num_SNB)%>%
dplyr::summarize(prop=mean(num_SNB), synch = sum(num_SNB),
n = n(), .groups='drop')
d<-
av_data1_all %>% ggplot( aes(x=num_SNB, y=prop, color=Age_SNB, group=Age_SNB))+
geom_smooth(av_data1_all,method=glm, method.args= list(family = binomial(logit)),
mapping = aes(y=prop,group=Age_SNB),linetype="dashed", position ="identity",se = FALSE)+ #mapping = aes(y=prop,group=Age, col='average'),linetype="dashed", position ="identity",se = FALSE)+
geom_smooth(
fix invalid codeSat, 22 Apr 2023 model <- lme(FOD_FODs ~ Condition_FODs, random = ~1|subject_FODs, data = indvar_FODs)
av_data1_all <- s_complete %>%
group_by(Gender_SNB, num_SNB) %>%
summarise(mean_num_SNB = mean(num_SNB),
sum_synch = sum(synch),
count = n()) %>%
ungroup()
avfourfods <- ggplot(meanFODs_indvar_FODs, aes(x = Condition_FODs, y = mean, fill = Condition_FODs)) +
geom_bar(stat = "identity", position = "dodge") +
geom_errorbar(aes(ymin = mean - sd/sqrt(count), ymax = mean + sd/sqrt(count)),
width = 0.2, position = position_dodge(width = 0.9)) +
xlab("Condition") +
ylab("Mean of FODs") +
ggtitle("Average of first order differences (FODs): Four conditions") +
theme_minimal() +
theme(panel.border = element_rect(color = "black", fill = NA, size = 1)) +
labs(fill = "Conditions") +
geom_hline(yintercept = 0, color = "black", linetype = "solid", size = 2)
# Write here the conclusions of your analysis,
# you can use code to illustrate the results of your analysis
# in this section, I analyze the independant variables.
#
m1<-glm(categorical_FOD_FODs ~ indvar_FODs$Condition_FODs + indvar_FODs$Language_used_FODs,data=indvar_FODs,family=binomial)
summary(m1)
library(lme4)
model = glmer(categorical_FOD_FODs ~ Condition_FODs * Gender_FODs , family=binomial(), data = df )
summary(model)
