library(agricolae)
library(car)
## Warning: package 'car' was built under R version 4.4.3
## Cargando paquete requerido: carData
library(multcompView)
## Warning: package 'multcompView' was built under R version 4.4.3
library(multcomp)
## Cargando paquete requerido: mvtnorm
## Cargando paquete requerido: survival
## Cargando paquete requerido: TH.data
## Cargando paquete requerido: MASS
## Warning: package 'MASS' was built under R version 4.4.3
##
## Adjuntando el paquete: 'TH.data'
## The following object is masked from 'package:MASS':
##
## geyser
library(PMCMR)
## Warning: package 'PMCMR' was built under R version 4.4.3
## PMCMR is superseded by PMCMRplus and will be no longer maintained. You may wish to install PMCMRplus instead.
##
## Adjuntando el paquete: 'PMCMR'
## The following object is masked from 'package:agricolae':
##
## durbin.test
library(PMCMRplus)
## Warning: package 'PMCMRplus' was built under R version 4.4.3
library("ggpubr")
## Warning: package 'ggpubr' was built under R version 4.4.3
## Cargando paquete requerido: ggplot2
library(nortest)
library(readxl)
## Warning: package 'readxl' was built under R version 4.4.3
#----------------------------------------------
# $\bar X$ y $S^2$ son variables aleatorias
#---------------------------------------------
set.seed(8)
n=100
tsim=1000
mu=3000
sigma2=300
sigma=sqrt(300)
#----------------------------------
# a) 1000 muestras de tama?o 100
#----------------------------------
muestra=matrix(0, nrow=tsim, ncol=n)
for(i in 1:tsim)
{
muestra[i,]=rnorm(n, mu, sigma)
}
dim(muestra)
## [1] 1000 100
summary(muestra[1, ])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2948 2988 2999 2998 3010 3041
boxplot(muestra[1, ], col=2)
media_muestral=NULL
var_muestral=NULL
sd_muestral =NULL
for(i in 1:tsim)
{
media_muestral[i]= mean(muestra[i,])
var_muestral[i] = var(muestra[i,])
sd_muestral[i] = sd(muestra[i,])
}
#------------------------------------
# b) Histograma de X_1, X_50 y X_100
#------------------------------------
hist(muestra[,1], freq=FALSE, xlab=expression(x[1]), main=expression(X[1]))
li=mu-4*sigma
ls=mu+4*sigma
peso=seq(li, ls, length=2000)
densidad_peso=dnorm(peso, mu, sigma)
hist(muestra[,50], freq=FALSE, xlab=expression(x[50]), main=expression(X[50]))
li=mu-4*sigma
ls=mu+4*sigma
peso=seq(li, ls, length=2000)
densidad_peso=dnorm(peso, mu, sigma)
##lines(peso, densidad_peso, col=2, lwd=2)
hist(muestra[,100], freq=FALSE, xlab=expression(x[100]), main=expression(X[100]))
li=mu-4*sigma
ls=mu+4*sigma
peso=seq(li, ls, length=2000)
densidad_peso=dnorm(peso, mu, sigma)
##lines(peso, densidad_peso, col=2, lwd=2)
#--------------------------------------------
# c) Histograma de Xbarra y (n-1)S2/sigma^2
#--------------------------------------------
hist(media_muestral, freq=FALSE, main=expression(bar(X)),
xlab=expression(bar(x)), ylab="Funci?n de densidad")
xbarra=seq(mu-4*(sigma/sqrt(n)), mu+4*(sigma/sqrt(n)), length=2000)
summary(xbarra)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2993 2997 3000 3000 3003 3007
densidad_xbarra=dnorm(xbarra, mu, sigma/sqrt(n))
##lines(xbarra, densidad_xbarra, col=2, lwd=2)
X2=((n-1)*var_muestral)/sigma2
hist(X2, freq=FALSE, main="X2", xlab="x2")
x2_valores=seq(60, 160, length=2000)
##lines(x2_valores, dchisq(x2_valores, df=(n-1)), col=2, lwd=2)
#----------
# d) Z y T
#----------
Z=(media_muestral-mu)/(sigma/sqrt(n))
T=(media_muestral-mu)/(sd_muestral/sqrt(n))
hist(Z, freq=FALSE)
z=seq(-3,3, length=2000)
##lines(z, dnorm(z,0,1), col=2, lwd=2)
hist(T, freq=FALSE)
t=seq(-3, 3, length=2000)
##lines(t, dt(t, df=(n-1)), col=2, lwd=2)
?dt
## starting httpd help server ...
## done
#-----------------------------------
# e) Intervalos de confianza del 95%
#-----------------------------------
li_ic=NULL
ls_ic=NULL
for(i in 1:tsim)
{
li_ic[i]=media_muestral[i]-1.96*sigma/sqrt(n)
ls_ic[i]=media_muestral[i]+1.96*sigma/sqrt(n)
}
ic=cbind(li_ic, ls_ic)
conteo=NULL
for (i in 1:tsim)
{
conteo[i]=ifelse(li_ic[i]<mu & ls_ic[i]>mu,1,0)
}
conteo
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [38] 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [75] 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1
## [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [149] 0 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0
## [186] 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [223] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1
## [260] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [297] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 0 1 1 1 1
## [334] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [371] 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [408] 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [445] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [482] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [519] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [556] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [593] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1
## [630] 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0 1 1
## [667] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [704] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [741] 1 1 1 1 1 1 1 1 1 1 0 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [778] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0 1 1 1 1 0
## [815] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [852] 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1
## [889] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1
## [926] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 1 1 1 1
## [963] 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 0 1 1
## [1000] 1
sum(conteo)/tsim
## [1] 0.96
#-------------------------------
# Prueba T-student una muestra
#-------------------------------
x=muestra[1, 1:20] # 20 datos de la normal con media mu=3000 y sigma=300
t.test (x, mu=3000, alternative="two.sided")
##
## One Sample t-test
##
## data: x
## t = -0.88378, df = 19, p-value = 0.3879
## alternative hypothesis: true mean is not equal to 3000
## 95 percent confidence interval:
## 2988.281 3004.761
## sample estimates:
## mean of x
## 2996.521
#----------------------------------------
# Prueba T-student de muestras pareadas
#----------------------------------------
#----------------------------------------------
# Prueba T-student dos muestras independientes
#----------------------------------------------
x=muestra[1, 1:20] # 20 datos de la normal con media mu=3000 y sigma=300
y=muestra[2, 25:40]
t.test (x, y, mu=0, alternative="two.sided", var.equal = TRUE)
##
## Two Sample t-test
##
## data: x and y
## t = -0.054671, df = 34, p-value = 0.9567
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -14.14263 13.40165
## sample estimates:
## mean of x mean of y
## 2996.521 2996.891
#------------------------------
# An?lisis de varianza (ANOVA)
#------------------------------
psist=c(180, 173, 175, 182, 181,
172, 158, 167, 160, 175,
163, 170, 158, 162, 170,
158, 146, 160, 171, 155,
147, 152, 143, 155, 160)
trat=c(rep("1", 5), rep("2", 5), rep("3", 5), rep("4", 5), rep("5", 5))
datos=as.data.frame(cbind(trat, psist))
plot(psist~trat, xlab="Tratamiento", ylab="Presi?n sist?lica", col=2, pch=16)
anova=aov(psist~trat)
summary(anova)
## Df Sum Sq Mean Sq F value Pr(>F)
## trat 4 2010.6 502.7 11.24 6.06e-05 ***
## Residuals 20 894.4 44.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
muestra=rf(1000, 4,20)
hist(muestra,freq=FALSE)
x=seq(0,15, length=1000)
plot(x, df(x, 4, 20), type="l", col=2)
#abline(v=11.24, lty=2, col=4)
alpha=0.05
vc=qf(1-alpha, 4, 20)
abline(v=vc, lty=2, col=5)
boxplot(anova$residuals, horizontal=TRUE, xlab="Residuales del ANOVA", col=2)
hist(anova$residuals, xlab="Residuales del ANOVA", col=4, ylab="Frecuencia", main="")
shapiro.test(anova$residuals)
##
## Shapiro-Wilk normality test
##
## data: anova$residuals
## W = 0.98911, p-value = 0.9927
# Normalidad
# Se prueba con base en los residuales del modelo
lillie.test(anova$residuals)
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: anova$residuals
## D = 0.064912, p-value = 0.9969
shapiro.test(anova$residuals)
##
## Shapiro-Wilk normality test
##
## data: anova$residuals
## W = 0.98911, p-value = 0.9927
# Homogeneidad de varianzas
require(graphics)
plot(psist ~ trat, data = datos)
#hartleyTest(psist~trat, data=datos)
bartlett.test(psist~trat, data=datos)
##
## Bartlett test of homogeneity of variances
##
## data: psist by trat
## Bartlett's K-squared = 2.6811, df = 4, p-value = 0.6125
leveneTest(psist ~ trat)
## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 4 0.483 0.748
## 20
# Pruebas de Comparaci?n M?ltiple
LSD.test(anova, "trat",console=TRUE)
##
## Study: anova ~ "trat"
##
## LSD t Test for psist
##
## Mean Square Error: 44.72
##
## trat, means and individual ( 95 %) CI
##
## psist std r se LCL UCL Min Max Q25 Q50 Q75
## 1 178.2 3.962323 5 2.990652 171.9616 184.4384 173 182 175 180 181
## 2 166.4 7.368853 5 2.990652 160.1616 172.6384 158 175 160 167 172
## 3 164.6 5.272571 5 2.990652 158.3616 170.8384 158 170 162 163 170
## 4 158.0 9.027735 5 2.990652 151.7616 164.2384 146 171 155 158 160
## 5 151.4 6.655825 5 2.990652 145.1616 157.6384 143 160 147 152 155
##
## Alpha: 0.05 ; DF Error: 20
## Critical Value of t: 2.085963
##
## least Significant Difference: 8.822417
##
## Treatments with the same letter are not significantly different.
##
## psist groups
## 1 178.2 a
## 2 166.4 b
## 3 164.6 b
## 4 158.0 bc
## 5 151.4 c
Tukey=TukeyHSD(anova, "trat", ordered = TRUE)
plot(Tukey, col=2)
