x <- rnorm(200,0,1.142)
z <- rnorm(200,0,2)
y <- factor(rbinom(200,1,(1/(1+exp(-1*(-2 + 1.5*x -0.5*z))))))
Pi <- matrix(data = c(0.9,0.1,0.3,0.7), nrow =2, byrow =FALSE)
dimnames(Pi) <- list(levels(y),levels(y))
ystar <- misclass(data.frame(y), list(y = Pi), k=1)[,1]
naive.model <- glm(ystar ~ x + z, family = binomial, x=TRUE, y =TRUE)
true.model <- glm(y ~ x + z, family = binomial)
simex.model <- mcsimex(naive.model, mc.matrix = Pi, SIMEXvariable = "ystar")
op <-par(mfrow = c(2,3))
invisible(lapply(simex.model$theta, boxplot, notch=TRUE, outline =FALSE, names=c(0.5,1,1.5,2)))
plot(simex.model)
par(op)
## example for a function which can be supplied to the function mcsimex()
## "xm" is the variable which is to be misclassified
my.mc <- function(datas,k){
xm <- datas$"xm"
p1 <- matrix(data = c(0.75,0.25,0.25,0.75), nrow =2, byrow = FALSE)
colnames(p1) <- levels(xm)
rownames(p1) <- levels(xm)
p0 <- matrix(data = c(0.8,0.2,0.2,0.8), nrow =2, byrow =FALSE)
colnames(p0) <- levels(xm)
rownames(p0) <- levels(xm)
xm[datas$y=="1"] <- misclass(data.frame(xm=xm[datas$y=="1"]),list(xm=p1), k=k)[,1]
xm[datas$y=="0"] <- misclass(data.frame(xm=xm[datas$y=="0"]),list(xm=p0), k=k)[,1]
xm <- factor(xm)
return(data.frame(xm))
}
Run the code above in your browser using DataLab