# 包括随机部分的微观模拟 GLM(microsimulation GLM including stochastic part)

``````set.seed(1)
library(MASS)

d <- mvrnorm(n=3000, mu=c(30,12,60), Sigma=matrix(data=c(45, 5, 40, 5, 15, 13, 40, 13, 300), nrow=3))
d[,1] <- d[,1]^2
``````

``````m <- glm(formula=d[,1]~d[,2] + d[,3], family=gaussian(link="sqrt"))
``````

``````p_linear <- m\$coefficients[1] + m\$coefficients[2]*d[,2] + m\$coefficients[3]*d[,3]
``````

``````sum((predict(object=m, type="response")-p_linear^2)^2)==0
``````

``````sd_residuals <- sd(p_linear^2 - d[,1])
p <- p_linear^2 + rnorm(n=1000, mean=0, sd=sd_residuals)
p_simulate <- simulate(m)\$sim_1
``````

``````par(mfrow=c(1,2), mar=c(4,2,2,1))
plot(d[,1], p, col="blue", pch=20, cex=0.5, xlim=c(-500, 3000), ylim=c(-500, 3000))
abline(lm(p~d[,1]), col="blue")
points(d[,1], p_simulate, col="red", cex=0.5)
abline(lm(p_simulate~d[,1]), col="red")
abline(a=0, b=1)
``````

...并比较两个预测：

``````plot(p, p_simulate, col="green", cex=0.5, xlim=c(-500, 3000), ylim=c(-500, 3000))
abline(lm(p_simulate~p), col="green")
abline(a=0, b=1)
``````