Nothing
var.strata <- function(strata, y = NULL, rh = strata$args$rh, rh.postcorr = FALSE,
model=c("none", "loglinear", "linear", "random"), model.control = list())
{
### Fonction externe : voir fiche d'aide pour la documentation
# Validation des arguments et initialisation de variables
# (Ls, takenone, certain, nmodel, beta, sig2, ph, gamma, epsilon, rhL) :
call.ext <- match.call()
out <- valid_args(obj_fct = as.list(environment()), call.ext)
# Variables que valid_args tire de strata
Ls <- out$Ls; takenone <- out$takenone; certain <- out$certain; L <- out$L; rhL <- out$rhL;
# Variables relatives au model
nmodel <- out$nmodel; beta <- out$beta; sig2 <- out$sig2; ph <- out$ph; pcertain <- out$pcertain;
gamma <- out$gamma; epsilon <- out$epsilon;
# Variable pour la sortie : liste des arguments
args <- out$args;
# Initialisation de variables tirees de strata :
##### ETAPE TRES IMPORTANTE
x <- if (!is.null(call.ext$y)) y else strata$args$x
# Si y a ete donne en entree, on doit faire les calculs sur y plutot que x, alors on modifie x
#####
N <- length(x) ## Nombre total d'observations
Nc <- length(certain)
Nnoc <- N - Nc ## nombre d'observations sans la strate certain (dans le vecteur xnoc)
xnoc <- if (is.null(certain)) x else x[-certain] ## observations sans la strate certain
strata.rhL <- if (takenone > 0) c(rep(1,length(takenone), strata$args$rh)) else strata$args$rh
stratumIDnoc <- if (is.null(certain)) strata$stratumID else as.numeric(as.character(strata$stratumID[-certain]))
Nh <- strata$Nh
bias.penalty <- if(is.null(strata$args$bias.penalty)) 1 else strata$args$bias.penalty
# Initialisation de quelques simples stat calculees sur les donnees (EX, EX2, EYc)
out <- init_stat(obj_fct = as.list(environment()))
EX <- out$EX; EX2 <- out$EX2; EYc <- out$EYc;
# Calcul des moments anticipees
moments <- .C(get_momentY_C, as.double(xnoc), as.integer(stratumIDnoc), as.integer(Nnoc),
as.integer(Nh), as.integer(L), as.integer(Nc), as.double(EYc), as.integer(takenone),
as.integer(nmodel), as.double(beta), as.double(sig2), as.double(ph), as.double(gamma),
as.double(epsilon), as.double(EX), as.double(EX2),
EYh = as.double(rep(0,L)), VYh = as.double(rep(0,L)), phih = as.double(rep(0,L)),
psih = as.double(rep(0,L)), TY = as.double(0), TAY = as.double(0), NAOK = TRUE)
EYh <- moments$EYh
VYh <- moments$VYh
TY <- moments$TY
TAY <- moments$TAY
# Ajustement a posteriori pour la non-reponse, si demande et si CV cible
if (rh.postcorr && !is.null(strata$args$n)){
warning("no posterior correction for non-response is available for stratified design with a target n")
dopostcorr <- FALSE
## Si le plan d'echantillonnage comporte un n cible, celui-ci ne serait plus atteint si on appliquait
## notre correction a posteriori pour non-reponse qui gonfle les nh. On ne le fait donc pas.
} else {
dopostcorr <- rh.postcorr
}
nhnonint <- if (dopostcorr) strata$nhnonint*strata.rhL/rhL else strata$nhnonint
nh <- if (dopostcorr) pmin(ceiling(nhnonint), Nh) else strata$nh
## Remarque : Ici, on n'arrondit pas nh de la meme facon que dans le reste du package.
# Pour le calcul d'autres valeurs pour la sortie
n <- get_n(nhcalcul = nh, obj_fct = as.list(environment()))
RRMSE <- NA ## car on veut que get_stat_out fasse le calcul du RRMSE
out_stat <- get_stat_out(obj_fct = as.list(environment()))
# Sortie des resultats
out <- list(nh=nh, n=n, nhnonint=nhnonint, certain.info = c(Nc = Nc, meanc = EYc), meanh = EYh, varh = VYh, mean=TY/N)
out <- c(out, out_stat)
out <- c(out, list(call = call.ext, date = date(), args = args))
class(out) <- "var.strata"
return(out)
}
print.var.strata <- function(x,...)
{
# Section des arguments fournis
cat("Given arguments:\n")
cat("strata = "); print(x$call$strata)
if (!is.null(x$args$y)) {
cat("y = "); print(x$call$y)
}
cat("rh.postcorr =",x$args$rh.postcorr)
if (is.null(x$args$y)) {
cat("\nmodel =",x$args$model)
nparam <- length(x$args$model.control)
if (nparam>0) {
cat(" : ")
for (i in 1:nparam) {
cat(names(x$args$model.control)[i],"=",x$args$model.control[[i]])
if (i<nparam) cat(" , ")
}
}
}
# Section du tableau de stratification
takenone <- if (is.null(x$args$strata$args$takenone)) 0 else x$args$strata$args$takenone
L<-x$args$strata$args$L+takenone
rh <- if(takenone==0) x$args$rh else c(NA,x$args$rh)
ph <- c(x$args$model.control$ptakenone,x$args$model.control$ph,x$args$model.control$pcertain)
type <- ifelse(x$nh == 0, "take-none", ifelse(x$nh == x$args$strata$Nh, "take-all", "take-some"))
tableau<-data.frame(rep("|",L),type,rh,x$args$strata$Nh,x$nh,x$nh/x$args$strata$Nh,rep("|",L),x$meanh,x$varh, stringsAsFactors=FALSE)
colnames(tableau) <- c("|","type","rh","Nh","nh","fh","|","E(Y)","Var(Y)")
if(!is.null(x$args$strata$args$certain)) {
tableau <- rbind(tableau,
list("|","certain",1,x$certain.info["Nc"],x$certain.info["Nc"],1,"|",x$certain.info["meanc"],NA))
}
tableau<-rbind(tableau,c(NA,NA,NA,sum(tableau$Nh),x$n,x$n/sum(tableau$Nh),NA,NA,NA))
rownames(tableau) <- if(!is.null(x$args$strata$args$certain)) c(paste("stratum",1:L),"","Total") else
c(paste("stratum",1:L),"Total")
if (identical(x$args$model,"loglinear")) tableau<-cbind("ph"=c(ph,NA),tableau)
tableau[,unlist(lapply(tableau,is.numeric))]<-round(tableau[,unlist(lapply(tableau,is.numeric))],2)
### Correction pour afficher correctement les NA
tableauc<-format(tableau)
for (i in 1:(dim(tableauc)[1]-1)) tableauc[i,] <- ifelse(substr(tableauc[i,],nchar(tableauc[i,])-1,nchar(tableauc[i,]))=="NA","-",tableauc[i,])
tableauc[dim(tableauc)[1],] <- ifelse(substr(tableauc[dim(tableauc)[1],],nchar(tableauc[dim(tableauc)[1],])-1,
nchar(tableauc[dim(tableauc)[1],]))=="NA","",tableauc[dim(tableauc)[1],])
### Fin de la correction
cat("\n\nStrata information:\n")
print(tableauc,na.print="")
cat("\nTotal sample size:",x$n,"\n")
cat("Anticipated population mean:",x$mean,"\n")
# Section sur les moments anticipes
sortie <- if (is.null(x$args$strata$args$takenone)) 1 else { if(0==x$args$strata$args$takenone) 2 else 3 }
if (sortie%in%c(1,2)) {
cat("Anticipated CV:",x$RRMSE,"\n")
if (2==sortie) cat("Note: CV=RRMSE (Relative Root Mean Squared Error) because takenone=0.\n")
} else {
est<-cbind(x$relativebias,x$propbiasMSE,x$RRMSE,x$args$strata$args$CV)
dimnames(est) <- if(length(est)==4) list(" ",c("relative bias","prop bias MSE","RRMSE","target RRMSE"))
else list(" ",c("relative bias","prop bias MSE","RRMSE"))
cat("\nAnticipated moments of the estimator:\n")
print.default(est, print.gap = 2, quote = FALSE, right=TRUE)
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.