====== RStudio ferramenta didática ======

Essa página é destinada a coleção de funções R para auxílio no ensino de estatística com o editor RStudio. O objetivo principal é reunir aqui implementações empregando a função ''manipulate::manipulate()'' para fazer gráficos interativos. Em segundo plano ficam dicas de como otimizar o uso desse editor.

Lista de implementações a fazer:
  * Gráfico de densidade de probabilidade com controle nos parâmetros;
  * Histograma com controle no número de classes, intervalo de classe e tipo de frequência;
  * ✔ Gráfico de densidade controlando o bandwidth e tipo de função kernel (Walmes);
  * ✔ Boxplot com controle no critério de representação dos extremos (Walmes);
  * Gráfico para ilustrar poder do teste controlando a diferença entre as médias;
  * ✔ Gráfico da densidade normal padrão com destaque para área acumulada até o quantil (Walmes);
  * Gráfico da reta ajustada e pontos ilustrando alavancagem;
  * Gráfico para ilustrar obtenção de valores iniciais para usar no ajuste de modelos de regressão não linear;
  * Gráfico qqplot alterando o valor de lambda da tranformação boxcox;
  * ✔ Gráfico tridimensional variando o ângulo de observação (Walmes);
  * Aproximação da binomial pela normal controlando o valor de p;
  * ✔ Outras Aproximações pela normal;
  * Convergência da média de realizações binomial, Poisson, beta, etc, para uma distribuição normal controlando tamanho da amostra;
  * ✔ Gráfico para estudo de medidas de influência em modelos de regressão linear (Walmes);
  * ✔ Teste de normalidade aplicado aos dados e aos resíduos (Walmes);
  * ✔ Teste de correlação para dados e resíduos de experimentos (Walmes);
  * ✔ Taxa de erro tipo I e tipo II (Walmes).

Dicas sobre o editor:
  * o atalho <fc #008000>''alt+-''</fc> faz o sinal de atribuição <fc #008000>''<''''-''</fc>;

==== Gráfico tridimensional variando o ângulo de observação ====

<code R>
# por Walmes ------------------------------------------------

require(manipulate) 
require(lattice)
da <- expand.grid(x=seq(0,10,l=30), z=seq(0,10,l=30))
da$y <- with(da, x+z+0.2*x*z) # gera dados

manipulate(
  ## faz o gráfico tridimensional
  wireframe(y~x+z, da,
            screen=list(z=z.angle, x=-60)),
  ## controla o valor do z.angle          
  z.angle=slider(0, 360, 10))

#------------------------------------------------------------
</code>

==== Gráfico da densidade normal padrão com destaque para área acumulada até o quantil ====

<code R>
# por Walmes ------------------------------------------------

require(manipulate)

manipulate(
  {
  curve(dnorm(x, 0, 1), -5, 5, ylab="f(x)")
  x <- seq(-5, q, by=0.05)
  fx <- dnorm(x, 0, 1)
  polygon(c(x, rev(x)),
          c(fx, rep(0, length(fx))),
          col="gray90")
  Pr <- round(pnorm(q, 0, 1), digits=3)
  legend("topleft", bty="n",
         legend=substitute(P(X<q)==Pr, list(q=q, Pr=Pr)))
  },
  q=slider(-5, 5, step=0.1, initial=0)
  )

#------------------------------------------------------------
</code>

==== Gráfico de densidade controlando o bandwidth e tipo de função kernel ====

<code R>
# por Walmes ------------------------------------------------

require(manipulate)

x <- rgamma(300, 3, 7)
manipulate(
  {
  plot(density(x, bw=bw, kernel=kernel))
  if(show.rug==TRUE) rug(x)
  },
  kernel=picker("gaussian", "epanechnikov", "rectangular",
                "triangular", "biweight","cosine",
                "optcosine"),
  bw=slider(0.01, 0.15, step=0.003, initial=0.05),
  show.rug=checkbox(TRUE, "show rug")
  )

#------------------------------------------------------------
</code>

==== Aproximações pela normal ====
<code R>
require(manipulate)
##
## Poisson e normal
##
manipulate(
  {
    XM <- floor(lam+4*sqrt(lam))
    curve(dpois(x, lambda=lam), 0, XM, n=XM+1,
          ylab="P[X=x] / f(x)", type="h")
    curve(dnorm(x, m=lam, sd=sqrt(lam)), 0, XM,
          add=T, col=2)
    legend("topright", c("Poisson", "Normal"), lty=1, col=1:2)
    title(substitute(lambda == l, list(l=lam)))
  },
  lam=slider(0.5, 30, step=0.5)
  )

##
## t e normal
##
manipulate(
    {
    curve(dt(x, df=df), -4, 4, ylim=c(0, 0.4), ylab="densidade f(x)")
    curve(dnorm(x), -4, 4, add=T, col=2)
    legend("topright", 
      c(substitute(t[nu == a], list(a=df)), expression(N(0,1))), 
      lty=1, col=1:2)
  },
    df = slider(1, 40)
  )

##
## Chi^2 e normal
##
manipulate(
  {
    curve(dchisq(x, df=df), 0, df+4*sqrt(2*df), 
          ylab="densidade f(x)")
    curve(dnorm(x, m=df, sd=sqrt(2*df)), 0, df+4*sqrt(2*df), 
          add=T, col=2)
    legend("topright", 
    leg=eval(substitute(c(expression(chi[nu==df]^2),
             expression(N(mu==df,sigma^2==df2))), 
             list(df2=2*df, df=df)),
),
          lty=1, col=1:2)
  },
  df=slider(1,50)
)
</code>

==== Gráfico de duas densidades normais com destaque para áreas acumuladas até o quantil ====
(Veja o enunciado deste exemplo em [[http://www.leg.ufpr.br/~silvia/CE055/node49.html]])

<code R>
# por Silvia ------------------------------------------------

manipulate(
  {
  curve(dnorm(x, 1.25, 0.12), 0.8, 2, ylab="f(x)",xlim=c(0.8,2))
  x <- seq(0.8, q, by=0.01)
  fx <- dnorm(x, 1.25, 0.12)
  polygon(c(x, rev(x)),
          c(fx, rep(0, length(fx))),
          col="gray90")
  Pr1 <- round(pnorm(q, 1.25, 0.12), digits=3)
  legend("topleft", bty="n",
         legend=substitute(P(X1<q)==Pr1, list(q=q, Pr1=Pr1)))
  curve(dnorm(x, 1.55, 0.13), 1.1, 2.0, ylab="f(x)",add=TRUE)
  x <- seq(0.8, q, by=0.01)
  fx <- dnorm(x, 1.55, 0.13)
  polygon(c(x, rev(x)),
          c(fx, rep(0, length(fx))),
          col="gray70")
  Pr2 <- round(pnorm(q, 1.55, 0.13), digits=3)
  legend("topright", bty="n",
         legend=substitute(P(X2<q)==Pr2, list(q=q, Pr2=Pr2)))
},
  q=slider(0.8, 2.0, step=0.01, initial=0.8)
  )
 
#------------------------------------------------------------
</code>

==== Boxplot com controle no critério de representação dos extremos ====

<code R>
# por Walmes ------------------------------------------------

require(manipulate)

manipulate(
  {
    x <- rep(1:10, 2)
    x[20] <- extreme
    gr <- gl(2, 10)
    bp <- boxplot(x~gr, outline=outline, range=range,
                  notch=notch, plot=FALSE)
    inf <- bp$stats[4,2]
    sup <- inf+range*diff(bp$stats[c(2,4),2])
    ylim <- extendrange(r=c(min(x), max(c(x,sup))), f=0.05)
    boxplot(x~gr, outline=outline, range=range,
            notch=notch, ylim=ylim)
    arrows(1.5, inf, 1.5, sup, angle=90, code=3, length=0.1)
  },
  extreme=slider(10, 30, step=0.5, initial=10),
  range=slider(1, 4, step=0.1, initial=1.5),
  outline=checkbox(TRUE, "show.outlier"),
  notch=checkbox(FALSE, "show.interval")
  )
  
#------------------------------------------------------------
</code>

==== Gráfico para estudo de medidas de influência em modelos de regressão linear ====

<code R>
# por Walmes ------------------------------------------------

require(manipulate)
data(anscombe)
ans0 <- anscombe[,c("x1","y1")]
ans0 <- ans0[order(ans0$x1),]
rownames(ans0) <- NULL
ans1 <- ans0
rx <- 2*diff(range(ans1$x1))
ry <- 2*diff(range(ans1$y1))
cols <- rep(1,nrow(ans0))

layout(matrix(c(1,2,1,3,4,5),2,3))

manipulate({
  ans1[po,] <- ans0[po,]+c(dx,dy)
  plot(ans1)
  points(ans1[po,], col="red", pch=19)
  abline(a=3, b=0.5, col="gray50", lty=2)
  m1 <- lm(y1~x1, data=ans1)
  abline(m1, col=2)
  h <- hatvalues(m1)[po]
  r <- residuals(m1)[po]
  legend("topleft", bty="n",
         legend=c(substitute(h==ha, list(ha=round(h,3))),
                  substitute(r==re, list(re=round(r,3)))))
  legend("bottomright", legend=po, bty="n")
  plot(m1)
  },
  po=slider(1,nrow(ans0), step=1, initial=5),
  dx=slider(-rx+0.001, rx, initial=0),
  dy=slider(-ry, ry, initial=0))

#------------------------------------------------------------
</code>

==== Teste de normalidade aplicado aos dados e aos resíduos ====

<code R>
# por Walmes ------------------------------------------------

par(mfrow=c(2,1))

manipulate({
  m <- rep(seq(0,by=h1,length.out=nlev), nrep)
  x <- rnorm(m, m, sd)
  xp <- qqnorm(x); qqline(x)
  rug(xp$x); rug(xp$y, side=2)
  legend("topleft", legend=shapiro.test(x)$p, bty="n")
  m0 <- lm(x~factor(m))
  xp <- qqnorm(residuals(m0)); qqline(residuals(m0))
  rug(xp$x);rug(xp$y, side=2)
  legend("topleft", bty="n",
         legend=shapiro.test(residuals(m0))$p)
  },
  h1=slider(0.001, 10, initial=1),
  nlev=slider(2, 15, initial=5),
  nrep=slider(2, 25, initial=5),
  sd=slider(0.01, 10, initial=1))
  
#------------------------------------------------------------
</code>

==== Teste de correlação para dados e resíduos de experimentos ====

<code R>
# por Walmes ------------------------------------------------

require(MASS)
par(mfrow=c(2,1))

manipulate({
  m <- rep(seq(0,by=h1,length.out=nlev), nrep)
  x <- mvrnorm(length(m), mu=c(0,0),
               Sigma=matrix(c(1,cor,cor,1),2,2))
  x[,1] <- x[,1]+m; x[,2] <- x[,2]+m
  plot(x)
  legend("topleft", legend=cor.test(x[,1], x[,2])$est, bty="n")
  m0 <- aov(x~factor(m))
  r <- residuals(m0)
  plot(r)
  legend("topleft", legend=cor.test(r[,1], r[,2])$est, bty="n")
  },
  h1=slider(0,19.99,initial=0.01),
  nrep=slider(10,300,initial=20),
  nlev=slider(2,15,initial=5),
  cor=slider(-0.99,0.99,initial=0))

#------------------------------------------------------------
</code>

==== Taxa de erro tipo I e tipo II  ====

<code R>
# por Walmes ------------------------------------------------

require(manipulate)

n <- 70
xx <- seq(qnorm(0.9), 4, l=n)
yy <- dnorm(xx, 0, 1)

area <- function(x,y){
  da <- rbind(cbind(x,y), c(x[1],y[1]))
  DA <- sapply(1:length(x),
               function(o){
                 dir <- da[o,1]*da[o+1,2]
                 esp <- da[o,2]*da[o+1,1]
                 c(dir, -esp)
               }
               )
  abs(sum(apply(DA, 1, sum)/2))
}

manipulate({
  curve(dnorm(x, 0, 1), -4, 9, col=1,
        ylim=c(-0.15,0.5), yaxt="n", ylab="f(x)")
  axis(2, at=seq(0,0.4,0.1))
  curve(dnorm(x, i, 1), col=2, add=TRUE)
  xx2 <- seq(-4+i, qnorm(0.9, i, 1)+2, l=n)
  yy2 <- dnorm(xx2, i, 1)
  xx2 <- c(min(xx2), xx2, max(xx2))
  yy2 <- c(0,yy2,0)
  yy2 <- pmin(yy2, dnorm(xx2, 0, 1))
  polygon(xx2, yy2, col="gray90")
  polygon(c(xx, rev(xx)),
          c(yy, rep(0, length(yy))), den=10)
  segments(i, 0, i, 0.4, col=2)
  segments(0, 0, 0, 0.4, col=1)
  tipo2 <- round(area(xx2,yy2), 3)
  text(2, -0.075, pos=1, label=expression(alpha==0.10))
  arrows(2, -0.075, 2, -0.01, length=0.1)
  text(i/2, 0.45, pos=3, label=substitute(beta==b, list(b=tipo2)))
  arrows(i/2, 0.45, i/2, max(yy2)+0.01, length=0.1)
  },
  i=slider(0, 6, step=0.01, initial=0)
  )

#------------------------------------------------------------
</code>