2015-01-23 3 views
0

У меня есть следующий код с помощью R:вычисление коэффициента GLM матрица, диагональные и среднеквадратическая ошибка

x0 <- matrix(rnorm(100,1)) 
x <- as.matrix(cbind("Intercept"=1, x0)) 
n <- dim(x0)[[1]] 
z <- cbind(rep(1,n),x0) 
p <- dim(x0)[[2]]+1 

for(i in 1:n) { 
    gstart <- glm(y~x0,family=binomial)$coef 
} 

Я хочу, чтобы вычислить оценки предыдущей обобщенной линейной модели в n образцах и создать матрицу оценок для в n случаи, то вычисления bias и mean square error, где матрица параметров задаются следующим кодом:

n=100 #is the number of samples 
parameter.mat<-cbind(rep(2,n),rep(2,n)) 
+2

Если 'y' быть определен где-то, и' x' и 'p' не используются и 'n' определяется дважды? Это немного запутывает – user20650

ответ

1

Я думаю,, что вы хотите изучить разницу между коэффициентами, возвращаемыми glm, и средними непараметрическими коэффициентами bootstrap. В приведенном ниже примере, первый дает способ, используя boot package, а затем с помощью цикла (похожий на ваш вопрос)

# some example data - set seed for reproducibility 
set.seed(1) 
dat <- data.frame(y = rbinom(100, 1, 0.5), x = rnorm(100)) 

# samples 
n <- 1000 

# glm estimates 
mod <- glm(y ~ x, family="binomial", data=dat)$coef 


# alternative method using boot package ----------------------------------- 
library(boot) 

# function to extract model coefficients 
f <- function(dat, i) glm(y ~ x, family="binomial", data=dat[i, ])$coef 

# run bootstrap 
set.seed(1) 
boot(dat, f, R=n) 

# manual bootstrap - sample with replacement ----------------------------- 
out <- vector("list", length=n) 

for(i in 1:n) { 
    newdat <- dat[sample(1:nrow(dat), , T), ] 
    out[[i]] <- glm(y ~ x, family="binomial", data=newdat)$coef 
    } 

# matrix of bootstrap coefficients 
bc <- do.call("rbind", out) 

# bootstrap means 
bc.mn <- colMeans(bc) 
bias <- mod - bc.mn