Does random effect always produce the same result as fixed effect? I have two models - one is including a categorical covariate as a fixed effect, the other includes it as a random effect:
require(nlme)
set.seed(123)
n <- 100
k <- 5
cat <- as.factor(rep(1:k, n))
cat_i <- 1:k # intercept per kategorie
x <- rep(1:n, each = k)
sigma <- 0.2
alpha <- 0.001
y <- cat_i[cat] + alpha * x + rnorm(n*k, 0, sigma)
plot(x, y)

m2 <- lm(y ~ cat + x)
summary(m2)

m3 <- lme(y ~ x, random = ~ 1|cat, na.action = na.omit)
summary(m3)

As you can see, both models m2 and m3 produce exactly the same coefficient estimate for x (including SE). Also the residual standard error is the same. The same result is produced when I simulate some missing data:
# simulate missing data
y[c(1:(n/2), (n*k-n/2):(n*k))] <- NA

m2 <- lm(y ~ cat + x)
summary(m2)

m3 <- lme(y ~ x, random = ~ 1|cat, na.action = na.omit)
summary(m3)

So can we say in general that adding effect as random will have the same impact on the other coefficients and the overall inference as adding it as fixed? If not, can you please provide a simple example (or change the provided one) when this fails?
 A: The ideas behind the answer are the same as in this answer. Again you are a victim of your design. Try the following code:
 require(nlme)
 set.seed(123)
 n <- 100
 k <- 5; 
 cat <- as.factor(sample(1:5, n*k, replace=T) ) #This should be a bit unbalanced.
 #cat <- as.factor(rep(1:k, n))
 cat_i <- 1:k # intercept per kategorie
 x <- rep(1:n, each = k)
 sigma <- 0.2
 alpha <- 0.001
 y <- cat_i[cat] + alpha * x + rnorm(n*k, 0, sigma)

 m2 <- lm(y ~ cat + x)
 m3 <- lme(y ~ x, random = ~ 1|cat, na.action = na.omit)
 round(digits=7,fixef(m3)[2]) == round(digits=7, coef(m2)[6]) 
 #Not in this case but results are really close because of your design matrix
 #If you did the same comparison with your original structure you would get TRUE.
 #    x 
 #FALSE 

 #y[c(1:(n/2), (n*k-n/2):(n*k))] <- NA
 # simulate missing data in an obviously unbalanced way
 ym <- y;  ym[1:66] <- NA
 m2m <- lm(y ~ cat + x) 
 m3m <- lme(ym ~ x, random = ~ 1|cat, na.action = na.omit)

 round(digits=7,fixef(m3m)[2]) == round(digits=7, coef(m2m)[6]) 
 #The results are further apart because your design changed more significantly.

As I said to you here, you are seem to ignore shrinkage in general and is one of the main reasons of you getting different estimates. Also, it needs to be emphasized that the main difference between a fixed and a random effect is that the first one is an unknown constant to estimate while a random effect is a random variable; it does not make sense to "estimate a value" for it. 
