статистический анализ циклической регрессии цикла и сохранение данных R

Question

Помогите мне с петлей? Я относительно новым для R. Короткая версия данных выглядит IKE это:

sNumber blockNo running TrialNo wordTar wordTar1 Freq Len code code2 
1  1  1  5   spouse violent 5011 6 1  2 
1  1  1  5   violent  spouse 17873 7 2  1 
1  1  1  5   spouse aviator 5011 6 1  1 
1  1  1  5   aviator  wife 515 7 1  1 
1  1  1  5    wife aviator 87205 4 1  1 
1  1  1  5   aviator  spouse 515 7 1  1 
1  1  1  9  stability usually 12642 9 1  3 
1  1  1  9   usually requires 60074 7 3  4 
1  1  1  9   requires  client 25949 8 4  1 
1  1  1  9   client requires 16964 6 1  4 
2  2  1  5   grimy  cloth 757 5 2  1 
2  2  1  5   cloth  eats 8693 5 1  4 
2  2  1  5    eats whitens 3494 4 4  4 
2  2  1  5   whitens  woman  18 7 4  1 
2  2  1  5   woman penguin 162541 5 1  1 
2  2  1  9    pie customer 8909 3 1  1 
2  2  1  9   customer sometimes 13399 8 1  3 
2  2  1  9  sometimes reimburses 96341 9 3  4 
2  2  1  9  reimburses sometimes  65 10 4  3 
2  2  1  9  sometimes gangster 96341 9 3  1 

У меня есть код для порядкового регрессионного анализа для одного участника в течение одного исследования (глаз отслеживания данных - eyeData), который выглядит как это:

#------------set the path and import the library----------------- 
setwd("/AscTask-3/Data") 
library(ordinal) 

#-------------read the data---------------- 
read.delim(file.choose(), header=TRUE) -> eyeData 

#-------------extract 1 trial from one participant--------------- 
ss <- subset(eyeData, sNumber == 6 & runningTrialNo == 21) 

#-------------delete duplicates = refixations----------------- 
ss.s <- ss[!duplicated(ss$wordTar), ] 

#-------------change the raw frequencies to log freq-------------- 
ss.s$lFreq <- log(ss.s$Freq) 

#-------------add a new column with sequential numbers as a factor ------------------ 
ss.s$rankF <- as.factor(seq(nrow(ss.s))) 

#------------ estimate an ordered logistic regression model - fit ordered logit model---------- 
m <- clm(rankF~lFreq*Len, data=ss.s, link='probit') 
summary(m) 

#---------------get confidence intervals (CI)------------------ 
(ci <- confint(m)) 

#----------odd ratios (OR)-------------- 
exp(coef(m)) 

Файл eyeData представляет собой огромный массив данных, состоящий из 91832 наблюдений с 11 переменными. Всего насчитывается 41 участник с 78 испытаниями каждый. В моем коде я извлекаю данные из одного теста от каждого участника, чтобы запустить anaysis. Однако для выполнения анализа вручную для всех испытаний для всех участников требуется много времени. Не могли бы вы помочь мне создать цикл, который будет читаться во всех 78 пробах из всех 41 участника и сохранить результаты статистики (я хочу сохранить сводку (m), ci и coef (m)) в одном файл.

Заранее благодарю вас!

Answer 1

Вы можете создать уникальный идентификатор для каждого испытания каждого участника. Затем вы можете перебрать все уникальные значения этого идентификатора и соответствующим образом подмножить данные. Затем вы запускаете регрессии и сохранить результат в виде объекта R

eyeData$uniqueIdent <- paste(eyeData$sNumber, eyeData$runningTrialNo, sep = "-") 
uniqueID <- unique(eyeData$uniqueIdent) 
for (un in uniqueID) { 
    ss <- eyeData[eyeData$uniqueID == un,] 
    ss <- ss[!duplicated(ss$wordTar), ] #maybe do this outside the loop 
    ss$lFreq <- log(ss$Freq) #you could do this outside the loop too 
    #create DV 
    ss$rankF <- as.factor(seq(nrow(ss))) 
    m <- clm(rankF~lFreq*Len, data=ss, link='probit') 
    seeSumm <- summary(m) 
    ci <- confint(m) 
    oddsR <- exp(coef(m)) 
    save(seeSumm, ci, oddsR, file = paste("toSave_", un, ".Rdata", sep = "")) 
    # add -un- to the output file to be able identify where it came from 
} 

Вариацию это может включать объединение выхода каждую итерации в списке (создать пустой список в начале), а затем после выполнения оценки и команды postestimation объединить элементы в списке и рекурсивно заполнить ранее созданные список «gatherRes»:

gatherRes <- vector(mode = "list", length = length(unique(eyeData$uniqueIdent) ##before the loop 
gatherRes[[un]] <- list(seeSum, ci, oddsR) ##last line inside the loop 

Если вы обеспокоены скоростью, вы могли бы рассмотреть возможность писать функцию, которая делает все это и использовать lapply (или mclapply).

Answer 2

Решение проблемы с использованием пакета plyr (оно должно быть быстрее, чем цикл for).

Поскольку вы не приводите воспроизводимый пример, я буду использовать данные iris в качестве примера.

Сначала сделайте функцию для расчета вашей интересующей статистики и возврата их в список. Например:

# Function to return summary, confidence intervals and coefficients from lm 
lm_stats = function(x){ 
    m = lm(Sepal.Width ~ Sepal.Length, data = x) 

    return(list(summary = summary(m), confint = confint(m), coef = coef(m))) 
} 

Затем используйте функцию dlply, используя переменные, представляющие интерес в качестве группировки

data(iris) 
library(plyr) #if not installed do install.packages("plyr") 

#Using "Species" as grouping variable 
results = dlply(iris, c("Species"), lm_stats) 

Это будет возвращать список списков, содержащий вывод summary, confint и coef для каждого вида.

Для вашего конкретного случая, функция может выглядеть следующим образом (не проверено):

ordFit_stats = function(x){ 

    #Remove duplicates 
    x = x[!duplicated(x$wordTar), ] 

    # Make log frequencies 
    x$lFreq <- log(x$Freq) 

    # Make ranks 
    x$rankF <- as.factor(seq(nrow(x))) 

    # Fit model 
    m <- clm(rankF~lFreq*Len, data=x, link='probit') 

    # Return list of statistics 
    return(list(summary = summary(m), confint = confint(m), coef = coef(m))) 
} 

А потом:

results = dlply(eyeData, c("sNumber", "TrialNo"), ordFit_stats)