library(stats)
library(nlme)
library(car)
### Simulation von heteroskesastischen Daten

x<- c(rep(0,50),rep(1,50))
x
y<- 1+ 2*x+0.5*rnorm(100)*(x==0)+ 2*rnorm(100)*(x==1)
mean(y[x==0])
sd(y[x==0])
mean(y[x==1])
sd(y[x==1])
xf<-as.factor(x)
### einfaches lineares Modell 
erg1<-(lm(y~xf))
summary(erg1)
plot(erg1)
### heteroskedatie erkennbar 

erg2<-gls(y~xf,weights=varIdent(form=~1|xf),method="REML")
summary(erg2)
plot(erg2)

### Schätzer identisch 
### jetzt mit weiterer Einflussgröße 
z<-rnorm(100)+x
y<- 1+ 2*x-0.5*z+0.5*rnorm(100)*(x==0)+ 2*rnorm(100)*(x==1)
erg3<-lm(y~xf+z)
summary(erg3)
library(effects)
plot(allEffects(erg3))
plot(erg3)
plot(erg3$residuals~xf)
### Heteroskedastie erkennbar###
erg4<-gls(y~xf+z,weights=varIdent(form=~1|xf),method="REML")
summary(erg4)
plot(erg4)
### einfache residuen ##
plot(erg4$residuals~xf)
#### Normalisierte residuen (mit der geschätzen Varianzmatrix normiert )
plot( residuals(erg4,type="normalized")~xf )
### Standardisierte Residuen
plot(residuals(erg4,type="pearson")~xf)
qqnorm(residuals(erg4,type="normalized"))
res<-residuals(erg4,type="normalized")