Я работаю над проектом анализа данных, который требует, чтобы я использовал симуляции Монте-Карло стандартной линейной регрессии. Для целей этого проекта переменные, которые мы ищем, - это размер эффекта и размер выборки моделируемых данных. Тем не менее, я новичок в выполнении симуляций Монте-Карло и могу использовать некоторую помощь в устранении неполадок моего кода.
Итак, чтобы понять проблему, я сначала использовал Vignette, который поставляется с пакетом MonteCarlo, просто чтобы понять, какпакет работает. https://cran.r -project.org / web / packages / MonteCarlo / vignettes / MonteCarlo-Vignette.html
Насколько я понимаю, пакет Monte Carlo запускает функцию, которую вы определили числомвремя, указанное вами, затем возвращает определенные пользователем результаты. Я попытался сгенерировать функцию, которая создала элементы, необходимые для функции MonteCarlo (): средство генерации данных, статистический метод, решение и результат. Модель завершается на 23%, чем вылетает. Код, который я собрал из разных источников в Интернете, представлен ниже.
####Creating a function to run regressions in a Monte Carlo simulation.
###Generating function to pull p-values out of regression models:
lmp <- function (modelobject) {
if (class(modelobject) != "lm") stop("Not an object of class 'lm' ")
f <- summary(modelobject)$fstatistic
p <- pf(f[1],f[2],f[3],lower.tail=F)
attributes(p) <- NULL
return(p)
}
###Creating function for regression analysis to be passed to the MonteCarlo() function.
Regression<-function(n, corr){
# create the initial x variable
x1 <- rnorm(n, 15, 5)
# x2, x3, and x4 in a matrix, these will be modified to meet the criteria
x2 <- scale(matrix( rnorm(n), ncol=1 ))
# put all into 1 matrix for simplicity
x12 <- cbind(scale(x1),x2)
# find the current correlation matrix
c1 <- var(x12)
# cholesky decomposition to get independence
chol1 <- solve(chol(c1))
newx <- x12 %*% chol1
# create new correlation structure (zeros can be replaced with other r vals)
newc <- matrix(
c(1 , corr,
corr, 1), ncol=2 )
# check that it is positive definite
chol2 <- chol(newc)
sample <- newx %*% chol2 * sd(x1) + mean(x1)
###Defining test statistic to be calculated on each sample.
stat<-lm(sample[,1]~sample[,2])
# get test decision:
decision<-abs(lmp(stat))<.05
# return result:
return(list("decision"=decision))
}
#define parameter grid for monte carlo simulation:
n_grid<-c(50,100,250,500)
corr_grid<-seq(0,.9,0.05)
# collect parameter grids in list:
param_list=list("n"=n_grid, "corr"=corr_grid)
####Running A Monte Carlo Simulation
library(MonteCarlo)
MC_result<-MonteCarlo(func=Regression, nrep=1000, param_list=param_list)
summary(MC_result)
Я надеялся получить объект MonteCarlo, который я мог превратить в таблицу, подобную той, что в виньетке Монте-Карло. , что-то вроде этого:
## \begin{table}[h]
## \centering
## \resizebox{ 1 \textwidth}{!}{%
## \begin{tabular}{ rrrrrrrrrrrrrrr }
## \hline\hline\\\\
## scale && \multicolumn{ 6 }{c}{ 1 } & & \multicolumn{ 6 }{c}{ 2 } \\
## n/loc & & 0 & 0.2 & 0.4 & 0.6 & 0.8 & 1 & & 0 & 0.2 & 0.4 & 0.6 & 0.8 & 1 \\
## & & & & & & & & & & & & & & \\
## 50 & & 0.06 & 0.28 & 0.81 & 0.99 & 1.00 & 1.00 & & 0.05 & 0.12 & 0.28 & 0.56 & 0.80 & 0.94 \\
## 100 & & 0.04 & 0.49 & 0.97 & 1.00 & 1.00 & 1.00 & & 0.06 & 0.19 & 0.53 & 0.83 & 0.98 & 1.00 \\
## 250 & & 0.05 & 0.88 & 1.00 & 1.00 & 1.00 & 1.00 & & 0.05 & 0.36 & 0.89 & 1.00 & 1.00 & 1.00 \\
## 500 & & 0.05 & 0.99 & 1.00 & 1.00 & 1.00 & 1.00 & & 0.05 & 0.59 & 1.00 & 1.00 & 1.00 & 1.00 \\
## \\
## \\
## \hline\hline
## \end{tabular}%
## }
## \caption{ decision }
## \end{table}
Вместо этого я получил следующую таблицу:
> ###Creating a LateX Table
> MakeTable(output=MC_result, rows="n", cols="corr", digits=3, include_meta=FALSE)
\begin{table}[h]
\centering
\resizebox{ 1 \textwidth}{!}{%
\begin{tabular}{ rrrrrrrrrrrrrrrrrrrrr }
\hline\hline\\\\
n/corr & & 0 & 0.05 & 0.1 & 0.15 & 0.2 & 0.25 & 0.3 & 0.35 & 0.4 & 0.45 & 0.5 & 0.55 & 0.6 & 0.65 & 0.7 & 0.75 & 0.8 & 0.85 & 0.9 \\
& & & & & & & & & & & & & & & & & & & & \\
50 & & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 \\
100 & & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 \\
250 & & 0.000 & 0.000 & 0.000 & 0.000 & 0.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 \\
500 & & 0.000 & 0.000 & 0.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 \\
\\
\\
\hline\hline
\end{tabular}%
}
\caption{ decision }
\end{table}
Проблема с этой таблицей в том, что она совершенно невероятна для предмета. Не имеет смысла, что линейная регрессия будет генерировать только 100% значимых p-значений или 0% значимых p-значений. Должна быть область таблицы, где некоторые комбинации корреляции и размера выборки приводят к тому, что некоторые значения p значимы, а некоторые значения p незначительны. То есть некоторые пропорции должны быть между 0 и 1, а не равны 0 и 1.
Мы будем признательны за то, что я делаю неправильно с этой функцией, чтобы привести к такому результату.