In ddarmon/confdist: confdist, a package for creating confidence functions
xbar <- 10
sigma <- 5
n <- 10
curve(pnorm((mu - xbar)/(sigma/sqrt(n))), from = -20, to = 20, xname = 'mu')
p-value for right sided hypothesis test
mean is a 10, sd = 5/sqrt(10)
pnorm(15 - xbar)/(sigma/sqrt(n))
curve(1-pnorm((mu - xbar)/(sigma/sqrt(n))), from = -20, to = 20, xname = 'mu')
p-value for left sided hypothesis test
1-pnorm(15 - xbar)/(sigma/sqrt(n))
6/23/20
Weight Loss Data
#install.packages("readr")
library(readr)
weight.change.data <- read.csv("Downloads/weight_change_data (1).csv")
Low Carb
hist(weight.change.data$weightchange[weight.change.data$diet=="Low Carb"])
Low Carb Confidence Distribution
xbar.lowcarb <- mean(weight.change.data$weightchange[weight.change.data$diet=="Low Carb"])
s.lowcarb <- sd(weight.change.data$weightchange[weight.change.data$diet=="Low Carb"])
n.lowcarb<- length(weight.change.data$weightchange[weight.change.data$diet=="Low Carb"])
pt.lowcarb <- pt(weight.change.data$weightchange[weight.change.data$diet=="Low Carb"], df= n.lowcarb-1)
pt.lowcarb
Low Carb Confidence Distribution Curve
curve.lowcarb <- curve(1-pt((xbar.lowcarb-mu)/(s.lowcarb/sqrt(n.lowcarb)), df = n.lowcarb-1), from = -20, to = 0, xname = 'mu')
curve.lowcarb
curve(1-pt((xbar.lowcarb-mu)/(s.lowcarb/sqrt(n.lowcarb)), df = n.lowcarb-1), from = -20, to = 0, xname = 'mu')
for (mu in seq(from=-20, to=0, by= 0.5) ) {
lowcarb.t.test <- t.test(weight.change.data$weightchange[weight.change.data$diet=="Low Carb"], alternative = "greater", mu = mu)
abline(h=lowcarb.t.test$p.value, v= mu)
}
Check
Low Carb Confidence Curve
cc.lowcarb <- function(mu) abs(2*(1-pt((xbar.lowcarb-mu)/(s.lowcarb/sqrt(n.lowcarb)), df = (n.lowcarb-1)))-1)
curve(cc.lowcarb, from = -20, to=0, xname = 'mu')
for (conf.int in seq(from = 0, to= 1, by =0.2)) {
lowcarb.t.test.2 <- t.test(weight.change.data$weightchange[weight.change.data$diet=="Low Carb"], conf.level= conf.int)
cc.lowcarb <- function(mu) abs(2*(1-pt((xbar.lowcarb-mu)/(s.lowcarb/sqrt(n.lowcarb)), df = (n.lowcarb-1)))-1)
abline(h=conf.int, v= lowcarb.t.test.2$conf.int)
}
Low Fat
hist(weight.change.data$weightchange[weight.change.data$diet=="Low Fat"])
xbar.lowfat <- mean(weight.change.data$weightchange[weight.change.data$diet=="Low Fat"])
s.lowfat <- sd(weight.change.data$weightchange[weight.change.data$diet=="Low Fat"])
n.lowfat <- length(weight.change.data$weightchange[weight.change.data$diet=="Low Fat"])
Low Fat Confidence Distribution
curve.lowfat <- curve(1-pt((xbar.lowfat-mu)/(s.lowfat/sqrt(n.lowfat)), df = n.lowfat-1), from = -20, to = 0, xname = 'mu')
curve.lowfat
Low Fat Confidence Curve
cc.lowfat <- function(mu) abs(2*(1-pt((xbar.lowfat-mu)/(s.lowfat/sqrt(n.lowfat)), df = (n.lowfat-1)))-1)
curve(cc.lowfat, from = -20, to=0, xname = 'mu')
6/24/20
Low Carb Confidence Density
low.carb.conf.dens <- function(mu) 1/(s.lowcarb/sqrt(n.lowcarb))*dt((xbar.lowcarb-mu)/(s.lowcarb/sqrt(n.lowcarb)), df = n.lowcarb -1)
curve(low.carb.conf.dens(mu), from = -20, to = 0, xname = "mu")
Check
integrate(function(mu) low.carb.conf.dens(mu), lower= -Inf, upper = -13)
curve(low.carb.conf.dens(mu), from = -20, to = 0, xname = "mu")
curve.lowcarb <- curve(1-pt((xbar.lowcarb-mu)/(s.lowcarb/sqrt(n.lowcarb)), df = n.lowcarb-1), from = -20, to = 0, xname = 'mu')
for (mu in seq(from = -20, to = 0, by = .5)) {
int.lowcarb <- integrate(function(mu) low.carb.conf.dens(mu), lower= -Inf, upper = mu)
abline(h= int.lowcarb$value, col = "blue")
abline(v= mu, col = "red")
abline(h = 0.5, v= -13.4225)
}
Confidence Distribution for Welch's Two-Sample T-Test
welch.curve.lowcarb <- curve(1-pt((xbar.lowcarb-xbar.lowfat-delta)/((s.lowcarb/sqrt(n.lowcarb))+(s.lowfat/sqrt(n.lowfat))), df = 601), from = -20, to = 0, xname = 'delta')
welch.curve.lowcarb
Check Welch Confidence Distribution
lowcarb.welch.test <- t.test(weight.change.data$weightchange[weight.change.data$diet=="Low Carb"],weight.change.data$weightchange[weight.change.data$diet=="Low Fat"], alternative = "greater", mu = delta)
lowcarb.welch.test
curve(1-pt((xbar.lowcarb-xbar.lowfat-delta)/((s.lowcarb/sqrt(n.lowcarb))+(s.lowfat/sqrt(n.lowfat))), df = 603.97), from = -10, to = 10, xname = 'delta')
for (delta in seq(from=-10, to=10, by= 1) ) {
lowcarb.welch.test <- t.test(weight.change.data$weightchange[weight.change.data$diet=="Low Carb"],weight.change.data$weightchange[weight.change.data$diet=="Low Fat"], alternative = "greater", mu = delta)
abline(h=lowcarb.welch.test$p.value, v= delta)
}
6/25/20
Confidence Curve for Two-Sample T-Test
cc.two.sample <- function (delta) abs(2*(1-pt((xbar.lowcarb-xbar.lowfat-delta)/((s.lowcarb/sqrt(n.lowcarb))+(s.lowfat/sqrt(n.lowfat))), df = 603.97))-1)
curve(cc.two.sample, from = -10, to = 10, xname = 'delta')
Density for Two-Sample T-Test
two.sample.conf.dens <- function(delta) 1/((s.lowcarb/sqrt(n.lowcarb))+(s.lowfat/sqrt(n.lowfat)))*dt((xbar.lowcarb-xbar.lowfat-delta)/((s.lowcarb/sqrt(n.lowcarb)+(s.lowfat/sqrt(n.lowfat)))), df = 603.97)
curve(two.sample.conf.dens(delta), from = -10, to = 10, xname = 'delta')
Check Density
integrate(function(delta) two.sample.conf.dens(delta), lower= -Inf, upper = -5)
qnorm(0.9995)
t.confdist
t.confdist <- function(x, y = NULL, paired = FALSE, plot = TRUE) {
if(!is.null(y)){
if(paired)
xok <- yok <- complete.cases(x,y)
else {
yok <- !is.na(y)
xok <- !is.na(x)
}
y <- y[yok]
} else {
if (paired) stop("'y' is missing for paired test")
xok <- !is.na(x)
yok <- NULL
}
x <- x[xok]
if (paired) {
x <- x-y
y <- NULL
}
xbar <- mean(x)
sx <- sd(x)
nx <- length(x)
vx <- var(x)
if(is.null(y)){
# Confidence Distribution
x1.dist.func <- function(mu) 1-pt((xbar - mu)/(sx/sqrt(nx)), df = (nx-1))
# Confidence Density
x1.dens.func <- function(mu) (1/(sx/sqrt(nx)))*(dt((xbar-mu)/((sx/sqrt(nx))), df=(nx-1)))
# Confidence Curve
x1.conf.curv.func <- function(mu) abs(2*(1-pt((xbar - mu)/(sx/sqrt(nx)), df = (nx-1)))-1)
# Confidence Quantile
x1.quantile.func <- function(p) xbar + (sx/sqrt(nx))*qt(p, df = (nx-1))
x.list <- list(pconf = x1.dist.func, dconf = x1.dens.func, cconf = x1.conf.curv.func, qconf = x1.quantile.func)
if(plot == TRUE){
t.crit <- qt(0.9995, df = nx-1)
curve(x1.dist.func(mu), from = xbar - t.crit*(sx/sqrt(nx)), to = xbar + t.crit*(sx/sqrt(nx)), xname = "mu", main = "Confidence Distribution")
curve(x1.dens.func(mu), from = xbar - t.crit*(sx/sqrt(nx)), to = xbar + t.crit*(sx/sqrt(nx)), xname = "mu", main = "Confidence Density")
curve(x1.conf.curv.func, from = xbar - t.crit*(sx/sqrt(nx)), to = xbar + t.crit*(sx/sqrt(nx)), xname = "mu", main = "Confidence Curve")
}
return(x.list)
} else if (!is.null(y)){
ybar <- mean(y)
sy <- sd(y)
ny <- length(y)
vy <- var(y)
stderrx <- sqrt(vx/nx)
stderry <- sqrt(vy/ny)
stderr <- sqrt(stderrx^2 + stderry^2)
df.welch <- stderr^4/(stderrx^4/(nx-1) + stderry^4/(ny-1))
# Confidence Distribution
welch.dist.func <- function(delta) 1-pt((xbar - ybar - delta)/(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), df = df.welch)
# Confidence Density
welch.dens.func <- function(delta) (1/sqrt((sx^2/(nx))+(sy^2/(ny)))*(dt((xbar - ybar- delta)/sqrt((sx^2/(nx))+(sy^2/(ny))), df=df.welch)))
# Confidence Curve
welch.conf.curv.func <- function(delta) abs(2*(1-pt((xbar - ybar - delta)/(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), df = df.welch))-1)
# Confidence Quantile
welch.quantile.func <- function(p) (xbar - ybar) + sqrt((sx^2/(nx))+(sy^2/(ny)))*qt(p, df = df.welch)
welch.list <- list(pconf = welch.dist.func, dconf = welch.dens.func, cconf = welch.conf.curv.func, qconf = welch.quantile.func)
if(plot == TRUE){
t.crit <- qt(0.9995, df = nx-1)
curve(welch.dist.func, from = (xbar-ybar) - t.crit*(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), to = (xbar-ybar) + t.crit*(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), xname = "delta", main = "Confidence Distribution")
curve(welch.dens.func, from = (xbar-ybar) - t.crit*(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), to = (xbar-ybar) + t.crit*(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), xname = "delta", main = "Confidence Density")
curve(welch.conf.curv.func, from = (xbar-ybar) - t.crit*(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), to = (xbar-ybar) + t.crit*(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), xname = "delta", main = "Confidence Curve", ylab = "Level of Confidence")
}
return(welch.list)
}
}
t.confdist.summary
t.confdist.summary <- function(mean,sd,n, plot = TRUE){
stopifnot(length(mean)==length(sd) && length(mean)==length(n))
if (length(mean)==1){
xbar <- mean[1]
sx <- sd[1]
nx <- n[1]
vx <- sx^2
# Confidence Distribution
x1.dist.func <- function(mu) 1-pt((xbar - mu)/(sx/sqrt(nx)), df = (nx-1))
# Confidence Density
x1.dens.func <- function(mu) (1/(sx/sqrt(nx)))*(dt((xbar-mu)/((sx/sqrt(nx))), df=(nx-1)))
# Confidence Curve
x1.conf.curv.func <- function(mu) abs(2*(1-pt((xbar - mu)/(sx/sqrt(nx)), df = (nx-1)))-1)
# Confidence Quantile
x1.quantile.func <- function(p) xbar + (sx/sqrt(nx))*qt(p, df = (nx-1))
x.list <- list(pconf = x1.dist.func, dconf = x1.dens.func, cconf = x1.conf.curv.func, qconf = x1.quantile.func)
if(plot == TRUE){
t.crit <- qt(0.9995, df = nx-1)
curve(x1.dist.func(mu), from = xbar - t.crit*(sx/sqrt(nx)), to = xbar + t.crit*(sx/sqrt(nx)), xname = "mu", main = "Confidence Distribution")
curve(x1.dens.func(mu), from = xbar - t.crit*(sx/sqrt(nx)), to = xbar + t.crit*(sx/sqrt(nx)), xname = "mu", main = "Confidence Density")
curve(x1.conf.curv.func, from = xbar - t.crit*(sx/sqrt(nx)), to = xbar + t.crit*(sx/sqrt(nx)), xname = "mu", main = "Confidence Curve")
}
return(x.list)
} else if (length(mean)==2){
xbar <- mean[1]
sx <- sd[1]
nx <- n[1]
vx <- sx^2
ybar <- mean[2]
sy <- sd[2]
ny <- n[2]
vy <- sy^2
stderrx <- sqrt(vx/nx)
stderry <- sqrt(vy/ny)
stderr <- sqrt(stderrx^2 + stderry^2)
df.welch <- stderr^4/(stderrx^4/(nx-1) + stderry^4/(ny-1))
# Confidence Distribution
welch.dist.func <- function(delta) 1-pt((xbar - ybar - delta)/(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), df = df.welch)
# Confidence Density
welch.dens.func <- function(delta) (1/sqrt((sx^2/(nx))+(sy^2/(ny)))*(dt((xbar - ybar- delta)/sqrt((sx^2/(nx))+(sy^2/(ny))), df=df.welch)))
# Confidence Curve
welch.conf.curv.func <- function(delta) abs(2*(1-pt((xbar - ybar - delta)/(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), df = df.welch))-1)
# Confidence Quantile
welch.quantile.func <- function(p) (xbar - ybar) + sqrt((sx^2/(nx))+(sy^2/(ny)))*qt(p, df = df.welch)
welch.list <- list(pconf = welch.dist.func, dconf = welch.dens.func, cconf = welch.conf.curv.func, qconf = welch.quantile.func)
if(plot == TRUE){
t.crit <- qt(0.9995, df = df.welch)
curve(welch.dist.func, from = (xbar-ybar) - t.crit*(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), to = (xbar-ybar) + t.crit*(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), xname = "delta", main = "Confidence Distribution")
curve(welch.dens.func, from = (xbar-ybar) - t.crit*(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), to = (xbar-ybar) + t.crit*(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), xname = "delta", main = "Confidence Density")
curve(welch.conf.curv.func, from = (xbar-ybar) - t.crit*(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), to = (xbar-ybar) + t.crit*(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), xname = "delta", main = "Confidence Curve", ylab = "Level of Confidence")
}
return(welch.list)
}else{
stop("vector lengths must be equal to or less than 2")
}
}
Testing with weight loss data
xbar.lowcarb.bmi <- -2.07
xbar.lowfat.bmi <- -1.75
s.lowfat.bmi <- ((xbar.lowfat.bmi)+1.97)*sqrt(n.lowfat)/(qt(0.975,df=n.lowfat-1))
s.lowcarb.bmi <- ((xbar.lowcarb.bmi)+2.3)*sqrt(n.lowcarb)/(qt(0.975, df = n.lowcarb-1))
s.lowcarb.bmi
s.lowfat.bmi
Testing with lowfat data
conf.lowfat.bmi <- t.confdist.summary(mean=xbar.lowfat.bmi, sd=s.lowfat.bmi, n=n.lowfat)
abline(h=0.95)
abline(v=-1.97)
abline(v=-1.52)
Testing with low carb data
conf.lowcarb.bmi <- t.confdist.summary(mean=xbar.lowcarb.bmi, sd=s.lowcarb.bmi, n=n.lowcarb)
abline(h=0.95)
abline(v=-2.3)
abline(v=-1.85)
Testing both
conf.lowcarb.lowfat.bmi <- t.confdist.summary(mean=c(xbar.lowcarb.bmi, xbar.lowfat.bmi), sd=c(s.lowcarb.bmi,s.lowfat.bmi), n=c(n.lowcarb,n.lowfat))
ddarmon/confdist documentation built on Aug. 8, 2020, 4:44 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
xbar <- 10 sigma <- 5 n <- 10
curve(pnorm((mu - xbar)/(sigma/sqrt(n))), from = -20, to = 20, xname = 'mu')
p-value for right sided hypothesis test
mean is a 10, sd = 5/sqrt(10)
pnorm(15 - xbar)/(sigma/sqrt(n))
curve(1-pnorm((mu - xbar)/(sigma/sqrt(n))), from = -20, to = 20, xname = 'mu')
p-value for left sided hypothesis test
1-pnorm(15 - xbar)/(sigma/sqrt(n))
6/23/20
Weight Loss Data
#install.packages("readr")
library(readr)
weight.change.data <- read.csv("Downloads/weight_change_data (1).csv")
Low Carb
hist(weight.change.data$weightchange[weight.change.data$diet=="Low Carb"])
Low Carb Confidence Distribution
xbar.lowcarb <- mean(weight.change.data$weightchange[weight.change.data$diet=="Low Carb"]) s.lowcarb <- sd(weight.change.data$weightchange[weight.change.data$diet=="Low Carb"]) n.lowcarb<- length(weight.change.data$weightchange[weight.change.data$diet=="Low Carb"])
pt.lowcarb <- pt(weight.change.data$weightchange[weight.change.data$diet=="Low Carb"], df= n.lowcarb-1) pt.lowcarb
Low Carb Confidence Distribution Curve
curve.lowcarb <- curve(1-pt((xbar.lowcarb-mu)/(s.lowcarb/sqrt(n.lowcarb)), df = n.lowcarb-1), from = -20, to = 0, xname = 'mu') curve.lowcarb
curve(1-pt((xbar.lowcarb-mu)/(s.lowcarb/sqrt(n.lowcarb)), df = n.lowcarb-1), from = -20, to = 0, xname = 'mu') for (mu in seq(from=-20, to=0, by= 0.5) ) { lowcarb.t.test <- t.test(weight.change.data$weightchange[weight.change.data$diet=="Low Carb"], alternative = "greater", mu = mu) abline(h=lowcarb.t.test$p.value, v= mu) }
Check
Low Carb Confidence Curve
cc.lowcarb <- function(mu) abs(2*(1-pt((xbar.lowcarb-mu)/(s.lowcarb/sqrt(n.lowcarb)), df = (n.lowcarb-1)))-1)
curve(cc.lowcarb, from = -20, to=0, xname = 'mu') for (conf.int in seq(from = 0, to= 1, by =0.2)) { lowcarb.t.test.2 <- t.test(weight.change.data$weightchange[weight.change.data$diet=="Low Carb"], conf.level= conf.int) cc.lowcarb <- function(mu) abs(2*(1-pt((xbar.lowcarb-mu)/(s.lowcarb/sqrt(n.lowcarb)), df = (n.lowcarb-1)))-1) abline(h=conf.int, v= lowcarb.t.test.2$conf.int) }
Low Fat
hist(weight.change.data$weightchange[weight.change.data$diet=="Low Fat"])
xbar.lowfat <- mean(weight.change.data$weightchange[weight.change.data$diet=="Low Fat"]) s.lowfat <- sd(weight.change.data$weightchange[weight.change.data$diet=="Low Fat"]) n.lowfat <- length(weight.change.data$weightchange[weight.change.data$diet=="Low Fat"])
Low Fat Confidence Distribution
curve.lowfat <- curve(1-pt((xbar.lowfat-mu)/(s.lowfat/sqrt(n.lowfat)), df = n.lowfat-1), from = -20, to = 0, xname = 'mu') curve.lowfat
Low Fat Confidence Curve
cc.lowfat <- function(mu) abs(2*(1-pt((xbar.lowfat-mu)/(s.lowfat/sqrt(n.lowfat)), df = (n.lowfat-1)))-1) curve(cc.lowfat, from = -20, to=0, xname = 'mu')
6/24/20
Low Carb Confidence Density
low.carb.conf.dens <- function(mu) 1/(s.lowcarb/sqrt(n.lowcarb))*dt((xbar.lowcarb-mu)/(s.lowcarb/sqrt(n.lowcarb)), df = n.lowcarb -1) curve(low.carb.conf.dens(mu), from = -20, to = 0, xname = "mu")
Check
integrate(function(mu) low.carb.conf.dens(mu), lower= -Inf, upper = -13)
curve(low.carb.conf.dens(mu), from = -20, to = 0, xname = "mu")
curve.lowcarb <- curve(1-pt((xbar.lowcarb-mu)/(s.lowcarb/sqrt(n.lowcarb)), df = n.lowcarb-1), from = -20, to = 0, xname = 'mu') for (mu in seq(from = -20, to = 0, by = .5)) { int.lowcarb <- integrate(function(mu) low.carb.conf.dens(mu), lower= -Inf, upper = mu) abline(h= int.lowcarb$value, col = "blue") abline(v= mu, col = "red") abline(h = 0.5, v= -13.4225) }
Confidence Distribution for Welch's Two-Sample T-Test
welch.curve.lowcarb <- curve(1-pt((xbar.lowcarb-xbar.lowfat-delta)/((s.lowcarb/sqrt(n.lowcarb))+(s.lowfat/sqrt(n.lowfat))), df = 601), from = -20, to = 0, xname = 'delta') welch.curve.lowcarb
Check Welch Confidence Distribution
lowcarb.welch.test <- t.test(weight.change.data$weightchange[weight.change.data$diet=="Low Carb"],weight.change.data$weightchange[weight.change.data$diet=="Low Fat"], alternative = "greater", mu = delta) lowcarb.welch.test
curve(1-pt((xbar.lowcarb-xbar.lowfat-delta)/((s.lowcarb/sqrt(n.lowcarb))+(s.lowfat/sqrt(n.lowfat))), df = 603.97), from = -10, to = 10, xname = 'delta') for (delta in seq(from=-10, to=10, by= 1) ) { lowcarb.welch.test <- t.test(weight.change.data$weightchange[weight.change.data$diet=="Low Carb"],weight.change.data$weightchange[weight.change.data$diet=="Low Fat"], alternative = "greater", mu = delta) abline(h=lowcarb.welch.test$p.value, v= delta) }
6/25/20
Confidence Curve for Two-Sample T-Test
cc.two.sample <- function (delta) abs(2*(1-pt((xbar.lowcarb-xbar.lowfat-delta)/((s.lowcarb/sqrt(n.lowcarb))+(s.lowfat/sqrt(n.lowfat))), df = 603.97))-1) curve(cc.two.sample, from = -10, to = 10, xname = 'delta')
Density for Two-Sample T-Test
two.sample.conf.dens <- function(delta) 1/((s.lowcarb/sqrt(n.lowcarb))+(s.lowfat/sqrt(n.lowfat)))*dt((xbar.lowcarb-xbar.lowfat-delta)/((s.lowcarb/sqrt(n.lowcarb)+(s.lowfat/sqrt(n.lowfat)))), df = 603.97) curve(two.sample.conf.dens(delta), from = -10, to = 10, xname = 'delta')
Check Density
integrate(function(delta) two.sample.conf.dens(delta), lower= -Inf, upper = -5)
qnorm(0.9995)
t.confdist
t.confdist <- function(x, y = NULL, paired = FALSE, plot = TRUE) { if(!is.null(y)){ if(paired) xok <- yok <- complete.cases(x,y) else { yok <- !is.na(y) xok <- !is.na(x) } y <- y[yok] } else { if (paired) stop("'y' is missing for paired test") xok <- !is.na(x) yok <- NULL } x <- x[xok] if (paired) { x <- x-y y <- NULL } xbar <- mean(x) sx <- sd(x) nx <- length(x) vx <- var(x) if(is.null(y)){ # Confidence Distribution x1.dist.func <- function(mu) 1-pt((xbar - mu)/(sx/sqrt(nx)), df = (nx-1)) # Confidence Density x1.dens.func <- function(mu) (1/(sx/sqrt(nx)))*(dt((xbar-mu)/((sx/sqrt(nx))), df=(nx-1))) # Confidence Curve x1.conf.curv.func <- function(mu) abs(2*(1-pt((xbar - mu)/(sx/sqrt(nx)), df = (nx-1)))-1) # Confidence Quantile x1.quantile.func <- function(p) xbar + (sx/sqrt(nx))*qt(p, df = (nx-1)) x.list <- list(pconf = x1.dist.func, dconf = x1.dens.func, cconf = x1.conf.curv.func, qconf = x1.quantile.func) if(plot == TRUE){ t.crit <- qt(0.9995, df = nx-1) curve(x1.dist.func(mu), from = xbar - t.crit*(sx/sqrt(nx)), to = xbar + t.crit*(sx/sqrt(nx)), xname = "mu", main = "Confidence Distribution") curve(x1.dens.func(mu), from = xbar - t.crit*(sx/sqrt(nx)), to = xbar + t.crit*(sx/sqrt(nx)), xname = "mu", main = "Confidence Density") curve(x1.conf.curv.func, from = xbar - t.crit*(sx/sqrt(nx)), to = xbar + t.crit*(sx/sqrt(nx)), xname = "mu", main = "Confidence Curve") } return(x.list) } else if (!is.null(y)){ ybar <- mean(y) sy <- sd(y) ny <- length(y) vy <- var(y) stderrx <- sqrt(vx/nx) stderry <- sqrt(vy/ny) stderr <- sqrt(stderrx^2 + stderry^2) df.welch <- stderr^4/(stderrx^4/(nx-1) + stderry^4/(ny-1)) # Confidence Distribution welch.dist.func <- function(delta) 1-pt((xbar - ybar - delta)/(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), df = df.welch) # Confidence Density welch.dens.func <- function(delta) (1/sqrt((sx^2/(nx))+(sy^2/(ny)))*(dt((xbar - ybar- delta)/sqrt((sx^2/(nx))+(sy^2/(ny))), df=df.welch))) # Confidence Curve welch.conf.curv.func <- function(delta) abs(2*(1-pt((xbar - ybar - delta)/(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), df = df.welch))-1) # Confidence Quantile welch.quantile.func <- function(p) (xbar - ybar) + sqrt((sx^2/(nx))+(sy^2/(ny)))*qt(p, df = df.welch) welch.list <- list(pconf = welch.dist.func, dconf = welch.dens.func, cconf = welch.conf.curv.func, qconf = welch.quantile.func) if(plot == TRUE){ t.crit <- qt(0.9995, df = nx-1) curve(welch.dist.func, from = (xbar-ybar) - t.crit*(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), to = (xbar-ybar) + t.crit*(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), xname = "delta", main = "Confidence Distribution") curve(welch.dens.func, from = (xbar-ybar) - t.crit*(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), to = (xbar-ybar) + t.crit*(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), xname = "delta", main = "Confidence Density") curve(welch.conf.curv.func, from = (xbar-ybar) - t.crit*(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), to = (xbar-ybar) + t.crit*(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), xname = "delta", main = "Confidence Curve", ylab = "Level of Confidence") } return(welch.list) } }
t.confdist.summary
t.confdist.summary <- function(mean,sd,n, plot = TRUE){ stopifnot(length(mean)==length(sd) && length(mean)==length(n)) if (length(mean)==1){ xbar <- mean[1] sx <- sd[1] nx <- n[1] vx <- sx^2 # Confidence Distribution x1.dist.func <- function(mu) 1-pt((xbar - mu)/(sx/sqrt(nx)), df = (nx-1)) # Confidence Density x1.dens.func <- function(mu) (1/(sx/sqrt(nx)))*(dt((xbar-mu)/((sx/sqrt(nx))), df=(nx-1))) # Confidence Curve x1.conf.curv.func <- function(mu) abs(2*(1-pt((xbar - mu)/(sx/sqrt(nx)), df = (nx-1)))-1) # Confidence Quantile x1.quantile.func <- function(p) xbar + (sx/sqrt(nx))*qt(p, df = (nx-1)) x.list <- list(pconf = x1.dist.func, dconf = x1.dens.func, cconf = x1.conf.curv.func, qconf = x1.quantile.func) if(plot == TRUE){ t.crit <- qt(0.9995, df = nx-1) curve(x1.dist.func(mu), from = xbar - t.crit*(sx/sqrt(nx)), to = xbar + t.crit*(sx/sqrt(nx)), xname = "mu", main = "Confidence Distribution") curve(x1.dens.func(mu), from = xbar - t.crit*(sx/sqrt(nx)), to = xbar + t.crit*(sx/sqrt(nx)), xname = "mu", main = "Confidence Density") curve(x1.conf.curv.func, from = xbar - t.crit*(sx/sqrt(nx)), to = xbar + t.crit*(sx/sqrt(nx)), xname = "mu", main = "Confidence Curve") } return(x.list) } else if (length(mean)==2){ xbar <- mean[1] sx <- sd[1] nx <- n[1] vx <- sx^2 ybar <- mean[2] sy <- sd[2] ny <- n[2] vy <- sy^2 stderrx <- sqrt(vx/nx) stderry <- sqrt(vy/ny) stderr <- sqrt(stderrx^2 + stderry^2) df.welch <- stderr^4/(stderrx^4/(nx-1) + stderry^4/(ny-1)) # Confidence Distribution welch.dist.func <- function(delta) 1-pt((xbar - ybar - delta)/(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), df = df.welch) # Confidence Density welch.dens.func <- function(delta) (1/sqrt((sx^2/(nx))+(sy^2/(ny)))*(dt((xbar - ybar- delta)/sqrt((sx^2/(nx))+(sy^2/(ny))), df=df.welch))) # Confidence Curve welch.conf.curv.func <- function(delta) abs(2*(1-pt((xbar - ybar - delta)/(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), df = df.welch))-1) # Confidence Quantile welch.quantile.func <- function(p) (xbar - ybar) + sqrt((sx^2/(nx))+(sy^2/(ny)))*qt(p, df = df.welch) welch.list <- list(pconf = welch.dist.func, dconf = welch.dens.func, cconf = welch.conf.curv.func, qconf = welch.quantile.func) if(plot == TRUE){ t.crit <- qt(0.9995, df = df.welch) curve(welch.dist.func, from = (xbar-ybar) - t.crit*(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), to = (xbar-ybar) + t.crit*(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), xname = "delta", main = "Confidence Distribution") curve(welch.dens.func, from = (xbar-ybar) - t.crit*(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), to = (xbar-ybar) + t.crit*(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), xname = "delta", main = "Confidence Density") curve(welch.conf.curv.func, from = (xbar-ybar) - t.crit*(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), to = (xbar-ybar) + t.crit*(sqrt((((sx)^2)/nx)+((sy)^2)/ny)), xname = "delta", main = "Confidence Curve", ylab = "Level of Confidence") } return(welch.list) }else{ stop("vector lengths must be equal to or less than 2") } }
Testing with weight loss data
xbar.lowcarb.bmi <- -2.07 xbar.lowfat.bmi <- -1.75 s.lowfat.bmi <- ((xbar.lowfat.bmi)+1.97)*sqrt(n.lowfat)/(qt(0.975,df=n.lowfat-1)) s.lowcarb.bmi <- ((xbar.lowcarb.bmi)+2.3)*sqrt(n.lowcarb)/(qt(0.975, df = n.lowcarb-1)) s.lowcarb.bmi s.lowfat.bmi
Testing with lowfat data
conf.lowfat.bmi <- t.confdist.summary(mean=xbar.lowfat.bmi, sd=s.lowfat.bmi, n=n.lowfat) abline(h=0.95) abline(v=-1.97) abline(v=-1.52)
Testing with low carb data
conf.lowcarb.bmi <- t.confdist.summary(mean=xbar.lowcarb.bmi, sd=s.lowcarb.bmi, n=n.lowcarb) abline(h=0.95) abline(v=-2.3) abline(v=-1.85)
Testing both
conf.lowcarb.lowfat.bmi <- t.confdist.summary(mean=c(xbar.lowcarb.bmi, xbar.lowfat.bmi), sd=c(s.lowcarb.bmi,s.lowfat.bmi), n=c(n.lowcarb,n.lowfat))
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.