title: Assessing the accuracy of a deterministic mathematical model

text: How can assess the accuracy of the output of a deterministic mathematical model? For example, a climate model can predict the mean annual temperature (MAT) for a specific location. I can use the model to predict thirty years of MAT in New York City, $T_\text{model}$. Now let's say I have the observed MAT in New York City for the same 30 years, $T_\text{obs}$. Is there a statistical test for the hypothesis $T_\text{obs}=T_\text{model}$? Can I assess the accuracy of the mathematical model?


title: Model-selection for linear mixed models over alternative sets of parameters (nlme function in R)

text: My models look like:               lme1 = lme(y~X+Y+V, random=~1|Subject, data=mydata, method ="ML")     lme2 = lme(y~X+Y+V2+V3, random=~1|Subject, data=mydata, method ="ML")     lme3 = lme(y~X+Y+V4, random=~1|Subject, data=mydata, method ="ML")      where X and Y are factors, but V, V2, V3,and V4 are continuous variables (modeled as covariates). I am using `Method ="ML"` in the hope that I could compare the likelihood values across the models. My research question has to do with whether V4 (in lme3) was a better predictor than V2 and V3 together, V2+V3 was better than V, etc. What goodness of fits measure is valid here? Can I use AIC values to compare models of different sets of parameters? I've also found some references on computing $R^2$ for mixed models. In particular, I am interested in the likelihood ratio test $R^2$ (Magee, 1990) which computes a $R^2$ by comparing each of these models to the null model. Using this method, I'd be comparing all three of my models to the same null model with just `y~1`. Is it then a valid approach to compare the $R^2$s generated? I am not a statistician but I would like to use a valid (at least justifiable) measure for my analysis. Any feedback would be greatly appreciated.


title: What is the time complexity of Lasso regression

text: What is the asymptotic time complexity of Lasso regression as either the number of rows or columns grows?


title: Means and imputation of log-normal variables

text: The geometric mean is the appropriate measure of central tendency for log- normally distributed variables. However, the arithmetic mean still has some use in relation to log-normal variables - in inferring totals from survey data, for instance. It is my understanding that the arithmetic mean of the sample of a log-normal distribtion is the same as the arithmetic mean of the population, with normal errors providing the sample is large enough for CLT to be invoked. If you operate a simple missing data imputation method of replacing with a mean, you would inflate the geometric mean if you used the arithmetic mean (and give high leverage to imputed data), and deflate the arithmetic mean (and therefore the total) if you used the geometric mean. My question is what is most appropriate to use, the geometric mean or the arithmetic mean? Or is such a simple imputation method completely inappropriate for log-normal data?


title: How to specify a multilevel regression using lmer function in R?

text: I am trying to analyse the effects of clipping (33% and 66%) and season (1, 2 and 3) on species richness. I have a control (no clipping, in no season) that the values are nested in. I am having difficulty building the code for this.   This is what I have managed:               lmm1 <- lmer(Species_Richness ~ Seasonf + Clipf + (Seasonf|Controlf) +          (Clipf|Controlf), data = Data)      How can I specify my model?


title: For chi-square on any 2 by X contingency table, should no more than 20% of the cells  be less than 5?

text: In their 1999 book, Yates, Moore & McCabe specify that no more than 20% of expected counts should be less than 5 and all individual expected counts should be greater than 1. Does this mean that, for any 2 by X contingency table that no more than 20% of the cells should be less than 5? So, for a 2 by 5 table no more than 2 (2 cells out of 10 = 20%) should be less than 5?


title: Best multivariate polynomial regression

text: I have a dataset (x,y) where x is a n-dimensional vector and y is an m-dimensional vector. (m=3, n>2) My goal is to find the best polynomial in x fitting the (x,y) dataset. The dimension of x is pretty big (right now it is 25), and I don't want to enter manually all the possibilities (ie x1*x3*x5, x1*x4*x6, ...). I can use Matlab, Mathematica and R. How can I do this? Also, I would be interested in hearing your suggestions about the following problem: how can I choose from the result the most relevant coefficients? (maybe x1*x2 is more relevant than x2*x3) Thank you


title: Is it possible to have only one strong factor loading to one factor in Factor Analysis?

text: I did a study which produced 6 factors. One of the factors has only 1 variable at 0.75, with the next closest on 0.54. My class was told to only include >0.6 for interpretation and omit the rest. Is it possible to have a factor representing one variable?


title: How to find cutoff values in multinomial regression?

text: **Updated** I have a data set with five continuous variables measuring an index in a population. For finding underlying classifications in the population based on the five observed variables, I did a latent profiles analysis (LPA) for two and three categories model using all five variables.The task is to find cutoff values for each variable in order to classify new observations without doing the whole procedure. Then I used produced categorization from LPA to find relative weights of each of the five variables using logistic (for two-category solution) and multinomial (for three-category solution) regression.That was for finding the relative contribution of each variable to the classification.Next step was to find cutoff values for each variable for future use. For finding cutoff values for each variable, I could find cutoff for each variable in two-category solution by using logistic regression and setting the regression equation equal to 0.5 and then finding corresponding value of the variable.In other words,I used predicted categorization by LPA as a criterion (DV) for each variable in logistic regression.I didn't use ROC because it didn't give me reasonable results. My question is, how to do the same thing for three-category solution? I tried multinomial regression but I didn't came up with a solution. with three categories I need to have two cutoff values in order to define category boundaries for each variable. Or simply, how to determine cutoff values for a variable based on a predefined categorizes to match with the categories?


title: How to calculate the probability of a random sequence of bits being a specific sequence?

text: If, for example, I have `0011` as a set of known bits $x$, how do I determine the probability that a sequence of randomly generated bits is equal to $x$? Thanks for the help! I'm sure this is a dumb question, but probability has always been my weakness (which my algorithms class is highlighting!). _EDIT:_ This is what I should've asked: How do you find the probability of randomly getting `0011` on the left-side of an 8-bit sequence?


Is there an analytic form for the Hellinger distance between von Mises distributions?


title: Repeated measures using R

text: How do I carry out a repeated measures analysis using R?


title: Two-Sample Problem Questions

text: I'm just preparing for my exam, but I have no idea how to solve this problem:   1. Durable Press cotton fabrics are treated to improve their recovery from wrinkles after washing. Unfortunately, the treatment also reduces the strength of the fabric. A study compared the breaking strength of fabrics treated by two commercial durable press processes. Five swatches of the same fabric were assigned at random to each process. Here are the data, in pounds of pull needed to tear the fabric. Permafresh: 29.9, 30.7, 30.0, 29.5, 27.6 Hylite: 28.8, 23.9, 27.0, 22.1, 24.2 c) Test the hypothesis null hypothesis U1 = U2 against the two-sided alternative twice: once using all the data and again without the outlier in the Permafresh sample. Do the two tests lead to similar conclusions? Can we safely conclude that one treatment has significantly higher mean breaking strength than the other? Help would really be appreciated. Thank you so much.


title: How can I efficiently approximate the sum of Bernoulli random variables for any number of summands in partial sum?

text: I have almost the same question as: How can I efficiently model the sum of Bernoulli random variables? But: (1) The number of random variables for summation is ~ N=20 (case 1) or N=90 (case 2). (2) $p_i$ ~ 0.13 (case 1) (3) The precision of the model based on Poisson law is not enough. (4) We need that our approx would be the good enough to model partial sums like these as well: $\sum_{i=k,N}{X_i}$, ( $k=1,N$ ) (5) We have empirical data for every $X_i$. The diagram shows that there is almost linear dependence for $Pr(X_i=1)$ for i=1,6 and then almost constant function or some small linear dependence (for i=7,20 ). Really I'm not sure in (3) since the Le Cam's inequality looks like very general... What class of model could we try?


title: McNemar's test implementation in Java

text: Does any body know of a Java implementation of McNemar's Test?


title: Constructing multilevel regression design matrix

text: A common formulation of multilevel/hierarchical regression models is $y = Xb + Zc + e$, where $X$ is an $n \times p$ matrix of $p$ individual level predictors, $Z$ is an $n \times q$ matrix of $q$ group level predictors, $y$ is an $n \times 1$ vector of observations at the individual level and $e$ is the error term. Suppose there are $J$ groups and the $q$ group level predictors are continuous. I would like to understand how to construct a full design matrix from this formulation, for example to apply a gradient descent optimization algorithm. Is this formulation equivalent to a design matrix constructed as: $p$ columns for the individual level predictors, $J$ columns of indicator variables for the $J$ groups (varying intercepts), $pJ$ columns for the interaction between the group indicators and the individual level predictors (varying slopes), and $q$ columns for the group level predictors (i.e., columns with repeated values for within-group observations). Are there alternative ways to construct the design matrix when the number of groups is large ($J>10000$) to reduce the number of parameters?


title: Regression model with multimodal outcome

text: All, I am trying to create a regression model where the (continuous) outcome is multimodal: ![enter image description here](http://i.stack.imgur.com/FLsNy.jpg) The outcome is the retail price of a certain product, and prices tend to fall around distinct amounts (750, 1000, 1250, 1500, etc). There are, however, a few prices in between so the prices are not distinct. I have run a linear model with satisfying results, though the extra prices between the modes give me pause. I also tried binning the prices down to a few groups representing the modes and it works somewhat well. Is there a better or worse way to model this? is there some sort of better or worse methodology for binning the outcome? Thank you


title: Calculating R-square for Elastic Net

text: I am trying to do "variable selection" using Elastic Net method (Matlab Lasso function with alpha of 0.5). I have 75 predictors in total (some are correlated with each other, hence using Elastic Net instead of Lasso), and I would like to get a subset of them, which are good predictors for my outcome. So my question is: How can I calculate something like R-Square that show how much of my outcome is explained by these selected variables? 1) If I use the selected variables in a multiple linear regression model, is the R-Square gonna be valid, since my variables are correlated? 2) Can I calculate cross-validated R-Square (using leave-one-out) to get a more accurate R-Square? 3) Is there any other way than calculating R-Square that I show my variable selection method predicts well?


title: What are the typical scales of multivariate regression?

text: I am about to receive a dataset for multivariate regression where we are trying to find $X\in\Re^{m\times n}$ such that it best fits the loss function $$\|AX-B\|_F^2$$ where $A\in \Re^{l\times m},B\in\Re^{l\times n}$. For some other purposes, what I was wondering was, what are the typical proportions for $A,X,B$ in most datasets? I am looking for ballpark estimates; I would assume $A$ is usually very tall but how many columns would it have? How many columns would $B$ have? If there are any freely available datasets for me to play with, if you can, please give link. I was not able to find any purely numerical datasets on the UCI ML repository.


title: Can AR models be written as MA models?

text: I heard that $AR(1)$ can be written as $MA(\infty)$. Can $AR(p), p \in \mathbb N, p \ge 2$ be written as $MA(q)$ for some $q$? Thanks!


title: Log-linear regression vs. logistic regression

text: Can anyone provide a clear list of differences between log-linear regression and logistic regression? I understand the former is a simple linear regression model but I am not clear on when each should be used.


title: Ratio estimation model in 2-stage cluster sampling

text: I've been reading about stratified sampling, 2-stage SRS sampling, and ratio estimation in finite populations and I have a question. When the ratio estimator is introduced, it seems that in order for it to perform well it is necessary that the population follows the following model (which I'll denote $\xi$): $$y_{ij} = \beta x_{ij} + \epsilon_{ij}$$ with the additional requirement that $E(\epsilon_{i,j}) = 0$ and $V(\epsilon_{i,j}) = \sigma^2 x_{i,j}$ where $x$ refers to the auxiliary information. I'm using subscripts $i, j$ since I'm thinking of this in the context of clusters. Here's my understanding of this issue as a whole: from the design-based perspective (perhaps I should call this the model-assisted design-based perspective?), inference about the ratio point estimate is based on all possible samples under the design (2-stage SRS) and is not at all based upon this model, although if this model is met by the population then we'll see much better performance. If instead we _assume_ this model (i.e. take the model-based route) then inference is based upon it and we actually need it to be true. As I understand it, this will lead to the same point estimate but potentially different variance estimates because we are either estimating the variance over all possible samples (design-based) or over all possible populations realizable under this model with these model parameters (i.e. $\beta$ and $\epsilon_{i,j}$). We'll also have that the ratio estimate is design-biased and model-unbiased. Please correct me if any of this is wrong. Here's my actual question. When we move on to cluster sampling, this model does not seem to be mentioned at all. As I understand it, the auxiliary information is the size of the clusters (denoted $M_i$). That the relationship passes through the origin, i.e. $y_{i,j} = 0$ iff $M_i = 0$, is clearly built in, but what about the variance relationship? I don’t see where variance proportional to cluster size comes in. Yet the ratio estimator seems to be viewed as a very reasonable (simple) choice so it seems that the conditions for it to perform well must be met. Any clarification would be tremendously appreciated.


title: Causal Analysis

text: I want to do "Nutritional Causal Analysis". Basically my objective is to identify possible factors associated with acute malnutrition in children under five years of age in developing countries. What types of factors are related with malnutrition? What will be conceptual framework of casual path of malnutrition?


title: ANOVA 'factor level comparison' methods

text: We had an exercise in statistics class where we did a small ANOVA analysis. We were comparing how the two levels of a factor influenced the response under different conditions (combinations of other factors). My analysis ended up as a tree where a leaf corresponded to each condition and I could see were/how my factor of interest affected the response. It was quite tedious making sense of the R outputs even in such a small case (and the 'tree' was in my mind). What techniques are usually used to make a 'factor-level comparison' like this in a more general case with lots of factors and levels? (Intuitive graphs and conclusions) Can it be generalized to ANCOVA?


title: Understanding formula for the standardized selection ratio (SSR)

text: Manyl et al defined an index called Standardized Selection Ratio (SSR). Here is the source. I don't quite understand how this index is calculated and how the p.values are calculated. Can you help me. The interesting part starts at page 40 (and ends few pages after) of the book. Below is an example of an article where this index was used. They described the SSR index but it doesn't make much sense to me either. * * * On this article they used the Standardized Selection Ratio (SSR) in order to know what is the preferred host (anemone) of anemonefish. Here is a quotation coming from their methods: > The “preferred host” of anemonefish was assessed calculating the > “Standardized Selection Ratio (SSR)” (values between 0 and 1) (Manly et al. > 1993). Manly’s standardized selection ratio represents the probability that > an individual will use a particular habitat type, taking into account the > different resource availability. For each anemonefish species ($i$) > inhabiting an anemone species ($j$), SSR was calculated as: $$SSR = \frac{w_i}{\sum{w_j}} $$ where $w_i = \frac{o_i}{p_j}$ > $o_i$ is the relative frequency of the anemonefish species $i$ and $p_j$ the > relative frequency of the anemone species $j$. Higher values of SSR indicate > a strong preference for the selected resource. The Log- Likelihood statistic > ($\chi^2$L) (Manly et al. 1993) was used to check the significance of the > observed distribution under a null hypothesis of a random host choice.


title: Would this be a problem that Survival Analysis can tackle? The problem deals with localized durations and uncensored data

text: I am trying to fit a model in R that doesn't have censoring. In other words, I know that each event will eventually happen. There are several experimental groups with people in each group receiving surgery. However, each group receives the surgery on different dates. I only have data for who died. Hence, my data set so happens to look like:               > data.frame     Group     Person   D = Day of Death since GROUP received surgery     Age     1         7        5                                                   45     1         8        12                                                  62     1         9        85                                                  55     1         10       102                                                 34     2         11       3                                                   44     2         12       7                                                   72     2         13       25                                                  81     2         14       56                                                  63     2         15       78                                                  57     3         16       2                                                   66     3         17       15                                                  28     3         18       88                                                  77     3         19       150                                                 54      The above is different than usual survival analysis data in that the days till death are dependent on the group a person belongs to. Meaning, each person within a group received the surgery on the same day, with each group's surgery date possibly different than another group's surgery date. Also, what is different here is that I only have data for people who died. Hence, there isn't any censoring (I believe). The Age variable is a covariates. If I were to try to model the hazard function as exponential, would it make sense to do this in R? Meaning, if I thought the hazard function looked like: $h_i(t) = \beta_0 + \beta_1Age_i + \beta_2log(Age_i)$ where the log is to account for any potential logarithmic effects. And if my survial function were: $S(t) = e^{-h_i(t)*t}$, would the following R model do just that? `survreg(Surv(D) ~ Age+log(Age), data = data.frame, dist ="exponential")` Would I have any troubles? I ask this because this isn't really your standard survival analysis problem. Thank you!


title: Understanding bootstrapping for validation and model selection

text: I think I understand how the fundamentals of bootstrapping work, but I'm not sure I understand how I can use bootstrapping for model selection or to avoid overfitting. For model selection, for example, would you just choose the model that yields the lowest error (maybe variance?) across its bootstrap samples? Are there any texts that discuss how to use bootstrapping for model selection or validation? **EDIT:** See this thread, and the answer by @mark999 for more context behind this question.


title: Explanation of picking an orthogonal array in the Taguchi method

text: I'm working on my thesis and I am trying to lower the amount of runs I need to do in an experiment so I thought I'd use the Taguchi method, however I don't understand how to use the selector table. I've been over and over how this is supposed to work and I still don't think I'm doing it correctly. I have a 3^k method where k = 5. My question stems from the fact that a teacher told me I could (should) use an L27 to have the interactions of 4 out of the 5 factors and I can't figure out how he thought of that given that, as I understand it, I should just be able to match the number of factors to my levels and that would give me my new array (L18) but apparently that's not the best for me to do.


title: Fractional integration and cointegration with R

text: I am currently reading about fractional integration in the context of error correction models. Does anyone know of a package in r that has a function for Sowells maximum likelihood estimator of d?


title: Incidence density sampling in R

text: I was wondering if there is an R package to perform incidence density sampling for a case-control design. That is, for matching in a case-control design; accounting for multiple matching covariates.


title: Mean Reciprocal Rank with GBM in R

text: Let's say I'm optimizing MRR with a GBM in R:               library(gbm)     generate.data <- function(N) {            # create query groups, with an average size of 25 items each       num.queries <- floor(N/25)       query <- sample(1:num.queries, N, replace=TRUE)            # X1 is a variable determined by query group only       query.level <- runif(num.queries)       X1 <- query.level[query]            # X2 varies with each item       X2 <- runif(N)            # X3 is uncorrelated with target       X3 <- runif(N)            # The target       Y <- X1 + X2            # Add some random noise to X2 that is correlated with       # queries, but uncorrelated with items            X2 <- X2 + scale(runif(num.queries))[query]            # Add some random noise to target       SNR <- 5 # signal-to-noise ratio       sigma <- sqrt(var(Y)/SNR)       Y <- Y + runif(N, 0, sigma)       Y <- ifelse(Y>median(Y), 1, 0)            data.frame(Y, query=query, X1, X2, X3)     }          set.seed(10)     data.train <- generate.data(1000)     gbm.mrr <- gbm(Y~X1+X2+X3,          # formula                     data=data.train,     # dataset                     distribution=list(   # loss function:                       name='pairwise',   # pairwise                       metric="mrr",     # ranking metric:                       max.rank=1,                       group='query'),    # column indicating query groups                     n.trees=2000,        # number of trees                     shrinkage=0.005,     # learning rate                     interaction.depth=3, # number per splits per tree                     bag.fraction = 0.5,  # subsampling fraction                     train.fraction = 1,  # fraction of data for training                     n.minobsinnode = 10, # minimum number of obs for split                     keep.data=TRUE,      # store copy of input data in model                     cv.folds=5,          # number of cross validation folds                     verbose = FALSE,     # don't print progress                     n.cores = 1)         # use a single core          best.iter.mrr <- gbm.perf(gbm.mrr, method='cv')     title('Training of pairwise model with mrr metric')      (Code adapted from here) This produces the following plot, where the black line is in-sample MRR, and the green line is out-of-sample MRR. Clearly, GBM is mimizing MRR: it chooses the model with the smallest out-of-sample MRR (~0.3): ![RPlot](http://i.stack.imgur.com/x07om.png) I'm confused as to why gbm MINIMIZES mrr. It's my understanding that a MRR of 1 is perfect (best document is the first result, 1/1=1) and a MRR of close to zero is very bad (e.g. best document is the 1000th result, 1/1000 is close to zero). Is my understanding of mmr correct? If so, then why does the gbm package minimize it?


title: How to assess overdispersion in Poisson GLMM, lmer( )

text: I have a GLMM with Poisson distribution and random spatial block. My experimental design is 2x2 factorial, with 4 blocks, resulting in 16 total data points. Here is the specification of the model in R using the lme4 package.               lmer(rich ~ morph*caged + (1|block), family=poisson, data=bexData)      When I call summary on this object, I am returned                  AIC   BIC logLik deviance      18.58 22.44 -4.288    8.576     Random effects:      Groups Name        Variance Std.Dev.      block  (Intercept)  0        0           Number of obs: 16, groups: block, 4      I have left out the fixed effect parameter tests and correlations for brevity. Here are my primary questions:   1. Can you use this output to calculate overdispersion?      * I have read that overdispersion can be calculated as the residual deviance divided by the residual degrees of freedom. Is that 8.576 / (16 - 4)? (Zuur et al., Mixed Effects Models)   2. If this calculation is correct, the estimator phi = 0.715. This indicates that there is not overdispersion in my data.      * Does this indicate that there is underdispersion?      * Is this a problem?      * Can anybody offer advice as to thresholds for over/underdispersion at which corrections to the models should be made? Zuur has said in one book that 5 is a common cutoff. Do people agree with that?     * How can such corrections be made?   3. I've also noticed here that the variance for the random effect is 0.      * Does this mean that there are precisely _no_ error correlations between data points within my blocking factor?     * If this is so, why would a generalised linear model of the form shown at bottom have an AIC substantially higher, around 55?     * is AIC a reasonable method for choosing GLMM over GLM (as suggested by Zuur)? .               glm(rich ~ morph*caged, data=bexData, family=poisson)


title: Why is my bootstrap function for paired samples t test in R not returning the same result as SPSS?

text: I want to verify using a bootstrap approach the hypothesis about a lack of the intervention effect within a group of six patients: pre-(X) post(Y) intervention measurement. My data               ID      X      Y     1  9.856  8.992     2 19.512  4.573     3  1.936  1.572     4 14.575  1.529     5  8.476 12.000     6  1.862  1.417      Within R (2.15.1) Ive written the following code, using the t-test for paired data, which relies on resampled pairs:               boot.p.value <- function(data, S)     {         boot.t.stat <- as.numeric()         t.stat <- t.test(x=data[,1], y=data[,2], paired=TRUE)$statistic             for(s in 1:S)             {                 boot.data <- data[sample(1:nrow(data), replace=TRUE),] ## resample pairs                 boot.t.stat[s] <- t.test(x=boot.data[,1],y=boot.data[,2], paired=TRUE)$statistic         }         p.value <- sum(1*(boot.t.stat >= t.stat))/S         return(p.value)     }      Where: > boot.p.value(data, S=1000) [1] 0.518 When repeated the resulting p-values values stay between .4 and .6. For the same data set the SPSS ver. 19 for the paired samples t-test provides bootstrap-based p = 0.182, for 1000 resamples. Why this difference?


title: Include lower level (2-way and 3-way) interactions in 3-level multilevel model?

text: I have a three level multilevel model (therapists, patients, repeated measures) in which I have included 4-way interactions, for instance: therapist self-efficacy * treatment condition * patient on track (yes/no) * time This tests the hypothesis that therapist who are in condition A and are higher no self-efficacy have a steeper slope with patients that are not on track Most of the lower order interactions are non-significant and also not that relevant for interpretation. Should I still include all of them? I feel like I would put a lot of parameters in the model, whereas it would be difficult to interpret the lower order interactions that are not significant. Most fora state that it would be adviceable to include lower order interaction, but there are also some places where it is stated that there is no statistical reason to include them. Any advice or references on this situation?


title: Factor analysis or just stick with compare means?

text: I would be grateful if somebody could assist me with the questions below This is for my dissertation (I am relatively new to SPSS and indeed to statistics). I have a small set of cases (N=41) of companies who took my questionnaire. There are 2 5-point Likert scale sets of questions: - one set is intended to determine whether the company uses performance measurement for a cost advantage strategy or a differentiation strategy (8 questions) - one set is intended to measure if the company finds performance measurement useful (6 questions) The companies have already been classified in the high (Hi) performers and low (Lo) performers. I have computed the mean scores for each group and was able to compare the means using the t-test for the two groups (Hi/Lo). I also did a factor analysis, decided on two factors for the 1st set (KMO=0.74) and one for the second and ended up with factor scores as variables. These range from say -2.6 to +1.3. Now I do not know what I can do with these results. I would really like to be able to determine that if a company scored high on the second set of questions (Usefulness) and high say on the cost advantage scale that this company used performance measurement for cost advantage strategy. Sort of a usefulness scale "enhancing" the other two (if this makes sense). I would then be able to find out if the high performer group follows a particular strategy or not (same for the low performers). Could regression or correlation do that? Thanks in advance, Nic


title: EM for truncated/censored data, how do those CDF's cancel?

text: I am practicing basic EM problems for qualifying exams. I am having trouble with truncated/censored data problems. There is something I don't understand about the complete-data likelihood. Consider the basic formulation: > Observed data X1...Xr and Z1 ... Zq have pmf/pdf f(~|p) and cdf F(~|p) > Observed data are all less than or equal to c and missing data are all > greater than c Then, the conditional densities of an X and a Z given their constraints w.r.t. c are then f(x|p)/F(c|p) I(X <= c) and f(z|p)/(1-F(c|p)) I(Z > c) respectively. From what I've seen in online examples, in the complete data likelihood those cdfs in the denominators somehow cancel. At least I think they do, because I tried working it with the cdfs, and at least in the case of the Poisson distribution, the M-step gets messy, since I have to take the derivative over p of logF(c|p) and log(1-F(c|p)) (log of a finite sum of the Poisson pmf), which I don't think have closed forms (could be wrong). In the _Grouping, Censoring, and Truncation_ section of Dempster's EM paper, they use a general case where the complete likelihood looks like a multinomial. Aside from what I have above (sets X and Z of size q and r) they allow you have an additional set of truncated missing values of unknown size, then include the unknown size as a missing data point. In that general case they also get a complete likelihood that doesn't have any of those pesky CDF- denominators, though they don't explain how they cancel. Can someone explain this? Thanks.


title: Proper Sampling - can I collect a two-group sample this way without issues?

text: I need to collect a two group sample for a comparison analysis (perhaps using logistic regression). The population that I need to extract a sample from is all firms from country A with activities in country B. The firms are classified into two categories: having a subsidiary in country B (S), or not having a subsidiary in country B (NS). I expect the share of S firms to be small relative to NS firms (but I have no way of knowing for sure). I already hold the entire population of S firms (because this data was available to me). However data on NS firms is not readily available and I have to collect that, and I will probably not get access to identify and collect all NB firms. So my situation is I have the entire population of S firms, and need to collect enough NS firms for subsequent analysis to be significant. Most likely my final sample will consist of all S firms and some share of the population of NS firms. Without much experience in doing these kinds of studies, I can't help to think that there is some kind is bias/reliability issue when sampling this way (one group: entire group population, other group: some part of group population). I have learned that _if_ it so happens that the population of NS firms is indeed much larger than S firms (again there is no way to know without data for the entire population of firms), and I e.g. end up with similar-sized samples of each group, there will be a case of oversampling the minority group. However I cannot find any remarks anywhere that consider this a problem for a comparison study, as a correct sample representation of the entire population is less important in this manner. Is my concern justified? Or is it fine to do it that way for e.g. logistic regression? If not, how can I get around the issue?


title: How to decide the optimal AR-model order?

text: I'm trying to create AR-model on wheather data and I wondered is there a method or algorithm which can find the optimal order for an AR-model? I'm using Matlab for my data-analysis, is there a function which can do this? Thank you


title: Randomness and Probability

text: Suppose the probability of a "random" event is very small (call this probability $p$). In real life, true randomness seems impossible. So would the actual true probability of the event be greater than $p$?


title: Modified ECDF in R

text: The ecdf (empirical cumulative distribution function in CDF) in R, instead of giving $P(X \le x)$ for a random variable $X$, it gives the proportion of observations in the data that are $\le X$. I tried to look up modified ECDF but couldn't find any. Is there any standard function to do the mathematical ECDF? Or any workaround would be appreciated!


title: Reference level and Interpretation

text: Suppose that we have a categorical variable $X$ that can take three values: $0$, $1$, or $2$. We use $X=0$ as the reference level. The following dummy variables are created: $$X_1 = 1 \ \text{if} \ X=1 \ \text{or} \ X_1=0 \ \text{if} \ X \neq 1$$ $$X_2 = 1 \ \text{if} \ X=2 \ \text{or} \ X_2=0 \ \text{if} \ X \neq 2$$ In a linear regression model, suppose we have $$Y = \beta_1 X_1 + \beta_2 X_2 + \epsilon$$ We know that $\beta_1$ compares $X=1$ vs $X=0$ and $\beta_2$ compares $X=2$ vs $X=0$. How would we get the coefficients that compares $X=2$ vs $X=1$? Would it be $\beta_1+\beta_2$?


title: Assumption for an M/M/1 queue

text: When a queueing system is modeled as an M/M/1 queue, it is assumed that the arrival time of jobs has Poisson distribution and the service rate has exponential distribution. I am wondering what features a system should have in order to model the arrival rate as Poisson? I known that Poisson is the only distribution that its inter-arrival time of jobs is exponentially distributed which is memoryless. Are there any better and more intuitive features for it? Also even in a more complex modeling (M/G/1), only the service rate is changed to general which means that the Poisson arrival rate is good enough where G/M/1 or G/G/1 is not as appealing as the previous ones.


title: Random block design ANOVA in R

text: Quick question: I have found this tutorial which recommends a two-factorial design for the following setup: three menu items (fixed factor) are tested in six restaurants (random factor). The guide recommends using               aov( response ~ item + restaurant )      I think that the correct way of doing that would be               aov( response ~ item + Error( restaurant / item ) )      I am confused. **My question: am I right, or am I wrong and the author of the tutorial is right?**


title: Higher Moments of an Unknown Density Function

text: Given a R.V $X$ and that $\ E(X) = 0 $ and $\ E(X^2) = \sigma^2 $. Is there anyway to compute $\ E(X^3) $ without knowing the density function of $X$?


title: Linear regression and arithmetic mean

text: I do experiments with a certain parameter x. The result is y. I assume y is linearly related to x. Suppose I can do 1000 experiments, which method will give me a better estimation of the linear relation?   * Select 1000 different values of x, get a single y for each x, and do linear regression?   * Select 100 different values of x, run 10 experiments for each x, average the y values for each x, and then do linear regression on the 100 averages?   * Select 100 different values of x, run 10 experiments for each x, and do linear regression without averaging first? What if I am not sure that the relation is linear?


title: Reference for generalized linear mixed models using Laplace approximation

text: I'm fitting a generalized linear mixed model in R using the Laplace approximation. I'm looking for a reference for the Laplace approximation used for that, or a reference regarding the comparison between the PQL and Laplace methods for fitting generalized linear mixed models.


title: Data preparation for logistic regression

text: This figure: ![odds ratios](http://i.stack.imgur.com/Q5Lrp.jpg) shows my odds ratios. I believe they have to be linear for logistic regression to work and am wondering how to pre-process these data. Thanks.


title: What could be the reason for using square root transformation on data?

text: Is there any reason of what I can think of, to transform the data with a square root? I mean what I always observe is that the R^2 increases. But this is probably just because of centering the data! Any thought is appreciated!


title: How do I find the best model with a saturated dataset?

text: I have hundreds of explanatory variables and under 100 observations (saturated data set). I'd like to create a linear model in which I have two or so composite variables made up of a dozen of the explanatory variables each. How do I find the best variables to use for the composites without going through every combination? Currently my $R^2 = .64$, I'd like to improve that. The model looks something like this: $$Y = B_1(v_1 + v_2 + \dots + v_{12}) + B_2(v_{13} + v_{14} + \dots +v_{24})$$ where $B_1$ and $B_2$ are the coefficients and the entirety of ($v_1 + v_2 + \dots + v_{12}$) is acting as one variable.


title: Left skewed vs. symmetric distribution observed

text: This is pretty hard for me to describe, but I'll try to make my problem understandable. So first you have to know that I've done a very simple linear regression so far. Before I estimated the coefficient, I watched the distribution of my $y$. It is heavy left skewed. After I estimated the model, I was quite sure to observe a left-skewed residual in a QQ-Plot as wel, but I absolutely did not. What might be the reason for this solution? Where is the mistake? Or has the distribution $y$ nothing to do with the distribution of the error term?


title: How to calculate variance of three or more attributes with known frequencies

text: Question was originally posted on Stack Overflow. I want to calculate variance of data at a given time where I know the frequency of the conditional and decision attributes. At time $t$ every conditional attribute can have 3 possible values while decsion attribute consists of two classes. See the picture below for more details: ![enter image description here](http://i.stack.imgur.com/XZnue.png)   * How can I calculate the variance of this data snapshot?   * Sample code or an algorithm would be most appreciated.


title: Extended binomial distribution

text: I'm trying to implement the extended binomial density function with support on c( 0 : (floor(N) + 1)), but I'm running into (I think) precision issues, as running:               ########################     #---DENSITY FUNCTION---#     ########################     debinom <- function(k, n, p, sum) {         if (k <=  n) {       return( choose(n, k) * p^k * (1-p)^(n-k) )       } else {         return (1.0 - sum)         }     }#END: pebinom          ##########################################     #---CUMULATIVE DISTRIBUTION FUNCTION 2---#     ##########################################     pebinom <- function(x, n, p) {            # point mass at 0       totalDensity = cumProb = debinom(0.0, n, p, 0.0)            k = 0       while (k <= (x)) {         density2 = debinom(k, n, p, totalDensity)         totalDensity = totalDensity + density2         cumProb = cumProb + density2         k = k + 1       }            k = k + 1       density = debinom(k, n, p, totalDensity)       cumProb = cumProb + density * (x - k)            return (cumProb)      }#END: debinom          ############     #---TEST---#     ############     for (i in 0:10) {     x = i + runif(1)     cat(x, " ", pebinom(x, 100, 0.1), "\n")     }      gives a negative probabilities for tail values. # EDIT I have changed, and mostly simplified the routines along the comments and answers I've received:               #########################################     #---PROBABILITY DISTRIBUTION FUNCTION---#     #########################################          debinom <- function(k, n, p) {          if (k <=floor(n)) {            return( choose(n, k) * p^k * (1-p)^(n-k) )            } else if(k == (floor(n)+1)) {              cumProb = 0.0         for(i in 0 : floor(n)) {           cumProb = cumProb + debinom(i, n, p)           }              return (1.0 - cumProb)              } else {            return(0.0)         }          }#END: pebinom          ########################################     #---CUMULATIVE DISTRIBUTION FUNCTION---#     ########################################     pebinom <- function(x, N, P) {          cumProb = 0     for(i in 0 : (floor(x)) ) {      cumProb = cumProb + debinom(i, N, P)     }          return(cumProb)     }


title: Probability of webpage accessed in a particular hour

text: I want to programatically calculate at which particular hours of a day the probability of a webpage hits(accessed) are high. Which statistic formula should I use to calculate the peak hours of the web page, if I have already the below data about a page . e.g. for page xyz, at left I have hours and at right I have hits , For different pages the hits are different. Page xyz hits count data:               hr=hits     1=0     2=0     3=0     4=0     5=14     6=0     7=0     8=5     9=5     10=8     11=10     12=10     13=12     14=7     15=5     16=5     17=3     18=0     19=0     20=0     21=0     22=0     23=0     24=0


title: How to test for independence in within-subject allele lengths?

text: I have data on a repeat length polymorphism. For each subject, I have two values, giving the (discrete) number of repeats on each allele. I want to test whether the length of the repeat on one allele and that on the other are independent of each other. My problem is that there is no natural order in the repeat lengths of one subject. If I wanted to do say a chi-squared test, there would be no justification for assigning one repeat length to the row variable and the other to the column variable of the contingency table I apply the test to, or vice versa. How do you test for independence of two-dimensional categorical data when, for each sampling unit, the two values are interchangeable?


title: Self-study question

text: I'm currently working on a self study worksheet. I understand most parts of the solution for part III, but I can't seem to make out how this comes about: ![enter image description here](http://i.stack.imgur.com/emSUF.jpg) **QUESTION:** ![enter image description here](http://i.stack.imgur.com/6g1QI.jpg) **ANSWER:** ![enter image description here](http://i.stack.imgur.com/ESOp2.jpg)


title: Fisher Exact Test Correlation Value

text: I have run the Fisher Exact test for several larger than 2x2 matrices. I got a Value number and a P value instead of just a P value. Does this mean the test ran out of memory and is invalid? or...Why did I get a Value?


title: How good is Monte Carlo Simulation when the variable distribution is unknown?

text: I am reading the book "how to measure everything", there is a chapter when the author encourages the usage of Monte Carlo simulation in simulating the future events in order to get a better understanding of the risk/ uncertainty involved. The theory sounds great in paper, and I am tempted to apply Monte Carlo techniques to startup analysis, to study the performance of employees, sales campaign, markets... everything! But in practice, I doubt whether Monte Carlo simulation is useful when the distribution of the variables involved in an equation is unknown. Think of the stock price return distribution-- it doesn't fit into a Gaussian Distribution, but what distribution we can use? And when the variable distribution is unknown, how useful is Monte Carlo simulation? My guess is not very useful, but I am right? A relevant, but probably narrower question is, is there any research done along the line of sensitivity analysis? By "sensitivity analysis" I mean we assume known distributions for the variables, but we perturb them a little by slightly changing the distribution shape or slightly changing some of the constants, will the outcome of the simulation results remain roughly the same, or completely different?


title: Single-pass algorithm for kurtosis

text: Here is a simple test I've run on MATLAB to check the validity of a single pass (online) algorithm for computing $3$rd moment and $4$th moment.               randn('state',0);          num2 = 0;     num1 = 0;     delta = 0;       M1 = 0;     M2 = 0;     M3 = 0;     M4 = 0;          Xvec = zeros(1, 100000);          for j = 1:length(Xvec)              X = 3*randn(1);              % Single pass algorithm            num2  = num1;         num1  = num1 + 1;         delta = X - M1;         delta_n  = delta / num1;         delta_n2 = delta_n * delta_n;         term1 = delta *delta_n * num2;         M1 = M1 + delta_n;         M4 = M4 + term1 * delta_n2 * (num1 * num1 - 3.0*num1 + 3.0) + ...         6.0 * delta_n2 * M2 - 4.0 * delta_n * M3;         M3 = M3 + term1 * delta_n * (num1 - 2.0) - 3.0 * delta_n * M2;         M2 = M2 + term1;              Xvec(j) = X;     end     % Quantities obtained from Single Pass     avg = M1;     variance = M2 / (num1 - 1.0);     kurtosis = (num1 * M4)/(M2 * M2) - 3.0;     skewness = sqrt(num1) * M3 / (M2^1.5);          % Reference Quantities      avg1     = sum /length(Xvec);     moment_2 = moment(Xvec, 2);     moment_3 = moment(Xvec, 3);     moment_4 = moment(Xvec, 4);      The online algorithm provides a correct mean (avg = avg1) and variance (variance = moment_2). However, the values for kurtosis and skewness obtained from the online algorithm are way off from the actual $3$rd and $4$th moments. What might be going wrong? Source


title: What happens in linear regression when observations are not independent in time

text: Let's consider the example of whether a person's weight is correlated with their height (simple linear regression). What is the difference between running this correlation on the following two data sets? **Measurements independent of time** Randomly sample 100 people and record their weight and height. **Measurements with time dependence** Randomly sample 10 people, but record their weight and heights at 10 different times.


title: How to visualize an interaction from a GLMM: use whole model to obtain predicted values or just factors in the interaction?

text: I need to plot an interaction between two continuous factors for interpretation. I am using the wireframe function in the lattice package and plotting the predicted values from the model with the best fit. The model looks somewhat like this:               out<-lmer(response~apples+pears+strawberry+mangoes+(mangoes)^2+strawberry:mangoes+       (1|id),family="poisson")      However, if I use the predicted values from the whole model the wireframe's surface is different than if I only use the predicted values from a model that only contains the variables that interact (omitting also the quadratic term):               out<-lmer(response~strawberry+mangoes+strawberry:mangoes+(1|id),family="poisson")      The difference is big enough to cause concern, and I would like to ask which method gives the best interpretation of the data, and why are the surfaces so different?


title: Cross level modeling in R, different levels of predictor, moderator and outcome variable

text: In the research that I am carrying out I have a moderation in which predictor, moderator and outcome variable are on the different levels: predictor and moderator on within level (the data were collected longitudinally over 6 months, every participant has 10 results) and outcome variable on between level (the data were collected just one time in the end of the study). I work with R. Is it possible to use Growth Modeling in this case? If not, is it possible to analyze this data using any different method?


title: Question about squares of the coefficients of variation

text: Suppose that $X$ is a random variable with any distribution that takes only positive values. Can the following inequality hold for constants $c_2 >c_1 >0$? $$\text{Var}(Z^{c_2})/\text{E}[Z^{c_2}]^2 > \text{Var}(Z^{c_1})/\text{E}[Z^{c_1}]^2$$ In other words: $$\text{E}[Z^{2c_2}]/\text{E}[Z^{c_2}]^2 > \text{E}[Z^{2c_1}]/\text{E}[Z^{c_1}]^2$$ where $Z=X+1$. With many thanks for your answer.


title: Information gain in random trees

text: When splitting attributes while constructing a random tree, I use information gain in order to determine the best value to split the tree on. I add nodes to the tree until a stopping criterion is met. What is the minimum value of information gain, to be used as a stopping criterion?


title: How to use reservoir states for readout and training?

text: I’m trying to make a Liquid State Machine, I have a spiking neural network as the liquid, and a feedforward neural network that should learn to map the reservoir’s states to the output. I’ve read many articles, but I didn’t find much about how I give the reservoir’s states to the feedforward network’s inputs. One of my sources said that I should sum all the reservoir’s values and give it to all input units, but that would result in the inputs having all the same value. Other articles said that the input units should be sensitive towards different features, but my project is recognizing handwriting and I have no idea what features I could extract from that. It only consists of values between 0 and 1, according to whether the pixel is filled or not. I tried to map the liquid’s different layers to different input units of the feedforward network, but it also didn’t work. I didn’t really think it was a good idea, but I tried it just in case. My question is: how do I give the reservoir's states to the feedforward network, so it could learn anything at all? I also might have problems with the liquid itself too, since it doesn’t do anything, only fires its data randomly. Should it also change its structure according to the input data? (The articles usually say that it doesn’t have to react to the input, but I think it’s strange.)


title: Time series: one I(1) and one I(0) variable, should I use VAR/VEC, test for cointegration?

text: Like the title says, I've got two time series, one is stationary to begin with and thus has no unit root, the other time serie is stationary after one-time differencing. I want to create a model out of this and I know that when unit roots are present, I should test for cointegration. But I've read in Engle & Granger (1987) that cointegration tests are only to be done when you have two or more I(1) variables, is that correct? So I cannot find in literature if I should now use a VAR model on differences or test for cointegration and perhaps do a Vector Error Correction model. Can anyone help me? I would be very thankful!


title: Test error correlation with label frequency

text: For a university project I wrote a SVM classifier and during the benchmark phase I'm getting some results I'd like to better understand from the theoretical side. My original dataset contains a multi-labeled examples but I'm buiding a binary classifier so I pick a feature `A` and I replace all the labels with +1 if `A` is present or -1 if it isn't. For instance:               label1 label2 label3 label4 = feature1 feature2 feature3 feature4     label2 label4  = feature1 feature2 feature3 feature4     label3 label4  = feature1 feature2 feature3 feature4      If choose `label3` the dataset becomes:                1 = feature1 feature2 feature3 feature4     -1  = feature1 feature2 feature3 feature4      1  = feature1 feature2 feature3 feature4      If choose `label2` the dataset becomes:                1 = feature1 feature2 feature3 feature4      1  = feature1 feature2 feature3 feature4     -1  = feature1 feature2 feature3 feature4      Now, my original dataset contains 20k rows and:   1. If I choose a label which is in 9k rows I get a average test error of 13%   2. If I choose a label which is in 5k rows I get a average test error of 7%  I was wondering why this is happening. My guesses are:   * In the second case the separating hyperplane has a larger margin    * In this the first case the is overfitting


title: What diagnostic plots exists for quantile regression?

text: Following on my question for OLS, I wonder: what diagnostic plots exists for quantile regression? (and are there R implementation of them?) A quick google search already came up with the worm plot (which I have never heard about before), and I'd be happy to know of more methods you might know about. (are any of them one from OLS ported for quantile-regression?) Thanks


title: Validation of data partitioning

text: I've some numerical sequence of data monitored over time $v(i)$, which belongs to one of two classes, 0 or 1. I monitor the variable several times, obtaining a sequence associated to 0 or 1. I use part of the cases in which $v(i)$ belongs to class 0 to have a profile/baseline (training data/cases). I then use the cases in which $v(i)$ belongs to 1 to see if there is some difference from the baseline (test data/cases). I end then in computing a metric, based on the information showed by $v(i)$, namely if it showed some "difference" when the case was a 1, or a 0 (error), and the viceversa considering the class 0. So from the whole cases $N$ I get just one value $M(v(i),N)$. But this value is just a statistic. My question is: How to guess the distribution of this metric? I thought to partition the "test data/cases" in $k$ groups (keeping the proportions), and then estimate the distribution of $M(v(i),N/k)$ through a $T$-test, or something similar. But does it make really sense? Is it a valid way to infer the distribution of $M$? Thanks in advance for the tips and discussion I will (in case) receive :)


title: Comparing components within a single time-series (physiologic data)

text: I have physiological data (electodermal activity) from several subjects recorded continuously during a task with multiple components. I'm interested in comparing whether the response within each component is different (for example: is the response in part D greater/different than in part A). Visually, there is a trend toward increased EDA as we get deeper into the task. I'm teaching myself time-series analysis as I go (reading Brockwell and Davis). My advisor suggested running a repeated-measures ANOVA on the mean EDA in each component with participant as the repeat so that we get a sense of main effect of component. My question is this: Would it be prudent to do anything in particular to the data (i.e. detrend it, etc) prior to calculating means within components? Additionally, (and I suppose more overarchingly) does my advisor's suggestion for comparing the components make functional sense (RM ANOVA). If this were anything but time-series I would think so but I'm new to TS so any advice would be helpful. I'm doing this primarily in R but have access to SAS as well. Thanks for any help anyone can be.


title: Why the most of world correlations are positive?

text: When analysing many real datasets, I have noticed that the most of correlations is positive. For example, here is a visualization of Pearson's correlation matrix for my current data. Is there any reason for this? ![enter image description here](http://i.stack.imgur.com/Yec8K.jpg) **EDIT:** Matrix is computed from 500*140 dataset, where the most of variables are Likert items. Here is the distribution of correlations. ![enter image description here](http://i.stack.imgur.com/UOvod.png) And here is example from another area. Vitamin and mineral contents were measured in different food samples (of the same weight). ![enter image description here](http://i.stack.imgur.com/cEd4i.png) ![enter image description here](http://i.stack.imgur.com/oZX8t.png)


title: Boosting algorithms and overfitting

text: In the book of Witten and Frank on _Data Mining algorithms_ , I read: _"If boosting succeeds in reducing the error on fresh test data, it often does so in a spectacular way. One very surprising finding is that performing more boosting iterations can reduce the error on new data long after the error of the combined classifier on the training data has dropped to zero. Researchers were puzzled by this result because it seems to contradict Occam’s razor, which declares that of two hypotheses that explain the empirical evidence equally well the simpler one is to be preferred. Performing more boosting iterations without reducing training error does not explain the training data any better, and it certainly adds complexity to the combined classifier. Fortunately, the contradiction can be resolved by considering the classifier’s confidence in its predictions."_ If I understand correctly, boosting can reduce the test error after the training error has dropped to zero. Does it mean that the classifier does overfit? Can you clear me this point?


title: Compare maxima of two Gaussian samples

text: Suppose $X$ and $Y$ are both normally distributed, with $X \sim \mathcal{N}(0,1)$ and $Y \sim \mathcal{N}(c,1),$ where $c > 0$. Consider $n$ independent draws of both $X$ and $Y$. As $n \rightarrow \infty,$ what is the probability that the sample maximum of the draws of $Y$ is greater than the maximum for $X$?


title: Looking for help with ITSM software (or other comparable software)

text: I'm looking for someone who is familiar with the ITSM software. I have some data that needs to be fit with an ARIMA/SARIMA model and then forecast using Holt-Winters/Seasonal method. I then need to compare them using RMSE and MAPE. I'm completely unfamiliar with ITSM and just need some help on the commands required to perform each of the above functions. The ITSM manual that comes with the download was not very helpful. If ITSM is not your forte, is there another easily available program that would perform all of these functions for me?


title: Specifying Dirichlet prior in rjags

text: I am running the below mixing model using the `rjags` package in R, but I keep getting the error message `RUNTIME ERROR:Cannot insert node into p[1,1]. Dimension mismatch`. I have clearly specified the Dirichlet prior incorrectly but I am unsure what needs changing.               model     {     # Likelihood     for(i in 1:N) {     Y[i,1:J] ~ dmnorm(mu[i,],Sigma.inv)     mu[i,1:J] <- p[i,1:K]%*%s[1:K,1:J,i]     }          # Priors on s     for(i in 1:N) {     for(k in 1:K) {       s[k,1:J,i] ~ dmnorm(mu.s[,k],Sigma.s.inv[,,k])       }     }          # Prior on p     for(i in 1:N) {      for(k in 1:K) {      p[k,i] ~ ddirch(alpha)     }     }          }     # alpha is specified as alpha<-c(1,1) in R        Any help will be much appreciated. As requested, here is the R code to run the model:               #Data for rjags mixing model          # Y values     Y<-matrix(c(5.8565, 18.5979, 5.7579,  1.1401,   0.4945, 0.2122, 0.3984,             14.0696,    10.6925,    8.5126, 2.1779, 1.1406, 0.4051, 0.5971,             13.7417,    10.0652,    8.2977, 2.2172, 1.1204, 0.3929, 0.6143,             11.7867,    15.1036,    8.3222, 1.9592, 0.9697, 0.3454, 0.5695), byrow=TRUE,ncol=7)     colnames(Y)<-c("Al","Ca","Fe","K","Mg","Na","Ti")          # mu.s values     mu.s<-matrix(c(6.9746333,  6.6574444,   10.5605000, 13.9441083,                 35.4727333, 13.2216667, 6.7945000,  4.1493958,                 5.0452333,  7.1281111,  6.1497083,  6.9320042,                 1.1962000,  1.1212222,  2.0860833,  2.4469958,                 0.6191000,  0.4853333,  1.0108333,  0.8803833,                 0.1965333,  0.1992222,  0.4911250,  0.4069042,                 0.4584000,  0.3072222,  0.6167083,  0.6584167),               byrow=TRUE, ncol=4)     rownames(mu.s)<-c("Al","Ca","Fe","K","Mg","Na","Ti")     colnames(mu.s)<-c("CB","FD","RV","TS")          # Sigma.inv values - a precision matrix     Sigma.inv<-structure(c(1291.8709,   6061.281,   698.1364,   -5510.728,   -690.3269, -3010.418,  641.0724,     6061.2809,  171684.070, 2179.5785,  -36806.221, 10430.4946, -13140.215, -93994.2303,     698.1364,   2179.578,   1705.0527,  -3500.690,  559.1736,   -3054.477,  -11664.9575,     -5510.7282, -36806.221, -3500.6902, 30963.025,  -5429.0766, 18493.588,  -10856.3964,     -690.3269,  10430.495,  559.1736,   -5429.077,  41562.2004, -38633.413, -18496.1137,     -3010.4185, -13140.215, -3054.4766, 18493.588,  -38633.4135,    116805.134, 21412.5527,     641.0724,   -93994.230, -11664.9575,    -10856.396, -18496.1137,    21412.553,  291117.7431), dim=c(7,7))     Sigma.inv<-as.positive.definite(Sigma.inv)          # Sigma.s.inv values - a precision matrix     Sigma.s.inv<-structure(c(17.1162359691711, 4.27333154820829, 0.799319413122213,              67.8990020405825, -181.19674061241, 240.147200630862, -235.997112580624,              4.27333154820829, 1.30433955546833, 0.346955743724074, 23.1902684301427,              -49.6382239122873, 59.017476611683, -63.7180758402929, 0.799319413122213,              0.346955743724074, 1.31697237634114, 5.28599480618941, -13.9928477614376,              -1.89552256841994, -9.54585305537158, 67.8990020405825, 23.1902684301427,              5.28599480618941, 577.794716546282, -1005.75939032625, 1125.51717030183,              -1523.09284200498, -181.19674061241, -49.6382239122873, -13.9928477614376,              -1005.75939032625, 2530.24588028147, -3264.56488630485, 2931.08069507173,              240.147200630862, 59.017476611683, -1.89552256841994, 1125.51717030183,              -3264.56488630485, 5517.74069570005, -4010.71091461381, -235.997112580624,              -63.7180758402929, -9.54585305537158, -1523.09284200498, 2931.08069507173,              -4010.71091461381, 5348.51616374144, 856.390529266247, -23.7499072549466,              30.8155799569968, -3121.79812866305, -6697.97335697594, 201.209069779657,              1119.41488802919, -23.7499072549466, 1.16945366195659, -0.680279236072279,              92.6596468824124, 185.589336069471, 8.55556473205667, -50.7878188653177,              30.8155799569968, -0.680279236072279, 1.33309308251978, -109.66739253505,              -247.377133710695, 20.475196391213, 35.4240283157601, -3121.79812866305,              92.6596468824124, -109.66739253505, 11599.9856016394, 24300.5608779793,              -610.557137045431, -4613.78940593557, -6697.97335697594, 185.589336069471,              -247.377133710695, 24300.5608779793, 52721.3237164215, -1907.99103757817,              -8706.61402757383, 201.209069779657, 8.55556473205667, 20.475196391213,              -610.557137045431, -1907.99103757817, 1924.48807974801, -150.765001020582,              1119.41488802919, -50.7878188653177, 35.4240283157601, -4613.78940593557,              -8706.61402757383, -150.765001020582, 3107.68558872587, 49.0622143602674,              8.44130661570893, -68.6376288719083, -256.909898229718, 45.2409921131372,              -122.516603322888, 505.474492036306, 8.44130661570893, 5.04250385971496,              -17.8396095728346, 0.998493378138782, -23.6324766148321, -21.9783007018327,              11.6566702138881, -68.6376288719083, -17.8396095728346, 122.63962551017,              281.428831916314, -9.59515390867994, 284.182865766167, -811.53358800718,              -256.909898229718, 0.998493378138782, 281.428831916314, 2315.44570687734,              -688.357300654226, 693.150782011221, -4275.03644110308, 45.2409921131372,              -23.6324766148321, -9.59515390867994, -688.357300654226, 1005.15772657742,              -887.017891953989, 937.328239057861, -122.516603322888, -21.9783007018327,              284.182865766167, 693.150782011221, -887.017891953989, 2701.25571899862,              -3519.5533300646, 505.474492036306, 11.6566702138881, -811.53358800718,              -4275.03644110308, 937.328239057861, -3519.5533300646, 13930.3108880326,              7.66040736195613, 1.14566743309344, 2.21528817899981, 12.2202312499515,              -75.5047038963214, 15.8162629586147, -144.642630413016, 1.14566743309344,              0.80708707189904, -0.219335995078617, 0.00356901053686257, -5.58097800246653,              10.9900105678523, 5.78098900706769, 2.21528817899981, -0.219335995078617,              4.25455776136105, 1.87986945691396, -11.0086131453236, 12.770373667221,              -71.2240549969279, 12.2202312499515, 0.00356901053686257, 1.87986945691396,              73.8970279130133, -192.535629882285, 38.158382504518, -564.560539731276,              -75.5047038963214, -5.58097800246653, -11.0086131453236, -192.535629882285,              1332.63815286144, 369.902304528778, 1013.5952607372, 15.8162629586147,              10.9900105678523, 12.770373667221, 38.158382504518, 369.902304528778,              1566.66721031925, -1476.55006583692, -144.642630413016, 5.78098900706769,              -71.2240549969279, -564.560539731276, 1013.5952607372, -1476.55006583692,              8509.10242328047), .Dim = c(7L, 7L, 4L))          #N     N<-nrow(Y)     #K     K<-4     #J     J<-7     #alpha     alpha<-c(1,1)          myjagsfile<-"JAGSmodel2.txt" # the above jags model     mydata <-list(N=N, J=J, K=K, Y=Y, mu.s=mu.s, Sigma.s.inv=Sigma.s.inv,   Sigma.inv=Sigma.inv, alpha=alpha)     mypars <-c("p", "s")     myiterations<-10000      myburnin<-5000     mythin<-125     set.seed(1000)     run.mcmc <- jags(data=mydata, inits=NULL, parameters.to.save=mypars, model.file=myjagsfile, n.chains=2,                   n.iter=myiterations, n.burnin=myburnin, n.thin=mythin, DIC=TRUE)


title: Confidence interval around binomial estimate of 0 or 1

text: What is the best technique to calculate a confidence interval of a binomial experiment, if your estimate is that p=0 (or similarly p=1) and sample size is relatively small, for example n=25?


title: How to test the significance of interquartile ranges(IQR)

text: There is a question about testing whether the interquartile ranges(IQR) from the two non-normal distributions are significantly different. Now I am looking for a test of significance testing interquartile ranges. I found some approaches for comparison of IQR, ie. Westenberg-Mood test permutation test boostrap + z test from http://r.789695.n4.nabble.com/significance-test- interquartile-ranges-td4636456.html Which test is appropriate for significance test interquartile ranges? (Data information:ratio data, the sample size is unequal in five groups,the least is less than 10, the most is about 300)


title: Can I use Principal Curves Analysis to fit a Vector Cloud instead of a Point Cloud?

text: I have recently discovered Principal Curves while trying to solve the problem I will describe below. The principle of Principal Curves is to fit a cloud point to find the "path" running along that points. My goal is to use something like this to fit a set of GPS tracklogs running along a road, so that the result would be a statistical aproximation of the "real" road (or, more precisely, the region of the road where people - bikers in this case - ACTUALLY ride, which most probably would be side lanes and road shoulders). The problem is: while classical Principal Curves Analysis consider the point set to have independent _points_ , I believe it would be more appropriate to consider each _track_ to be composed of _dependent vectors_. When using points, for example, a track segment with high sampling rate might cause a bias when compared to another with lower sampling rate, although each of them describe only ONE trajectory each, and I am more interested in weighting the trajectories with one another instead of the individual points. (by the way, does that rationale of mine make sense?) I posted a question on the same topic, but with different formulation, in the Gis.StackExchange site, with some representative pictures: http://gis.stackexchange.com/questions/70623/how-can-i-statistically- calculate-the-real-road-from-a-set-of-gps-tracks Thanks for reading!


title: What does "PCA (Principal Component Analysis) spheres the data" mean?

text: I was reading some notes and it says that PCA can "sphere the data". What they define to me as "sphering the data" is dividing each dimension by the square root of the corresponding eigenvalue. I am assuming that by "dimension" they mean each basis vector into which we are projecting (i.e. the eigenvectors we are projecting to). Thus I guess they are doing: $$ u^{'}_i= \frac{u_i}{\sqrt{eigenValue(u_i)}}$$ where $u_i$ is one of the eigenvectors (i.e. one of the principal components). Then with that new vector, I am assuming they are projecting the raw data we have, say $x^{(i)}$ to $z^{(i)}$. So the projected points would now be: $$ z'^{(i)} = u^{'}_i \cdot x^{(i)}$$ They claim that doing this ensures that all the features have the same variance. However, I am not even sure if my interpretation of what they mean by sphering is correct and wanted to check if it was. Also, even if it was correct, what is the point of doing something like this? I know they claim it makes sure everyone has the same variance but, why would we want to do this and how does it achieve this?


title: Choice of path weights in SEM conceptual model using openMx

text: I am reviewing the R package OpenMx for a genetic epidemiology analysis in order to learn how to specify and fit SEM models. I am new to this so bear with me. I am following the example on page 59 of the OpenMx User Guide. Here they draw the following conceptual model ![enter image description here](http://i.stack.imgur.com/3jOmk.png) And in specifying the paths, they set the weight of the latent "one" node to the manifested bmi nodes "T1" and "T2" to be 0.6 because: > The main paths of interest are those from each of the latent variables to > the respective observed variable. These are also estimated (thus all are set > free), get a start value of 0.6 and appropriate labels.               # path coefficients for twin 1     mxPath(       from=c("A1","C1","E1"),       to="bmi1",       arrows=1,       free=TRUE,       values=0.6,       label=c("a","c","e")     ),          # path coefficients for twin 2     mxPath(       from=c("A2","C2","E2"),       to="bmi2",       arrows=1,       free=TRUE,       values=0.6,       label=c("a","c","e")     ),      The value of 0.6 comes from the estimated covariance of bmi1 and bmi2 (of strictly *mono*zygotic twin pairs). I have two questions: 1) When they say that the path is given a "starting" value of 0.6 is this like setting a numerical integration routine with initial values, like in estimation of GLMs? 2) Why is this value estimated strictly from the monozygotic twins?


title: Can first-order Markov chain be considered a special case of a hidden Markov model?

text: I am trying to apply R depmixS4 package in order to cluster time series with model based clustering. The model consists of K components, each being a first order Markov models. The Expectation-Maximization algorithm is then used to estimate model parameters. My time series are multivariate and of arbitrary length ( i.e. can be 400, can be 1). Now, there are 2 problems:   1. depmixS4 is oriented towards Hidden Markov Models, not the basic first order ones.    2. I do not completely understand the E and M-steps of model based clustering when applied to first-order Markov model components. Overwhelmingly, scientific literature talks about Hidden Markov Models.  However, it seems to me that simple first-order Markov models can be seen as a particular case of hidden ones (where the response probability distribution is just identity, response equals state - thus the states are actually visible). So the clustering process does not have to decode the hidden states. However, the E-M process itself is still unclear, and I am not sure if I can apply the depmixS4 package methods by adjusting them or should I develop my own algorithm in R. The mixture Markov model used in my research, is the following: $$p(v | \Theta) = \sum_{t=1\ldots K}p(c_k)p_k(v | \Theta_k),$$ where $v=v_1,\ldots,v_L$ - vector of arbitrarily length, and $$p_k(v |\Theta_k)=p(v_1| \Theta_{k_i})\prod_{i=2\ldots L} p(v_i |v_{i-1}, \Theta_{k_T}).$$ Given are also vector of initial state prior distribution and transition probability matrix $T_p$.


title: Weights in quantile regression for complex survey in R

text: I want to include sample weights to my quantile regression model, but I'm not sure how to do this. I've already define my weight, which are replicated weights already given in survey dataset (computed in survey package):               w<-svrepdesign(variables=data[,1:10],repweights=data[,11:30],type="BRR",        combined.weights=TRUE, weights=r.weights, rho=0.5,dbname="")      and my rq model is:               rq(y~x,tau=c(.1,.2,.3,.4,.5,.6,.7,.8,.9),data=my.data))      I tried to use `withReplicates` function, but with no success. Any suggestions?


title: Machine Learning - What to Use as the "Correct Answer" for training purposes

text: I apologize if this is too simple a question for this forum but... If I'm trying to build a neural network or other machine learning based prediction of who will win a football game and by how many points given such and such inputs do I use the actual score of past games as the correct answer for my neural network for training purposes. For example if I'm using winning percentage of each team to predict winners. I might have the following:               [Home Team Win Pct] [VisitingvTeam Win Pct]   [Home Team Win Margin (points)]           50%                           45%                 +2      50%                           75%                 -6      ....                         .....               ....      ...and so on from past games - I might have 10 years of game data let's say. Can I use the actual outcome point differential ([Home Team Win Margin]) as the "correct answer" for training purposes? Since I have no way of knowing if these two teams played 100 times under the same exact conditions if the average winning margin would be +2 or some other number since for a given game result I only have a sample size of one. Or is there some other measure of "correct answer" that I should use instead? I certainly can't think of one. Thanks for your help.


title: Fixed effect insignicant results

text: Currently iam working on my master thesis which is about risk adjusted returns in the Asia Pacific REIT market. The goal of the paper is to determine/find variables that conceive explanatory power over the IV. To determine this I also performed a FE regression based on countries. The result are more or less in line with FE for the whole sample (with firmID as dummy) and OLS regression. The suprising thing actually is that for the whole sample period (so still with country dummies) some of the variables seem to be insignificant, while on the contrary, if the sample is split up in 2 sub- periods, they show to have a significant effect in both the periods. I find this rather strange that a variable is significant for both the sub- periods, but fail to do so for the whole sample period. What could be the reason the reason for this? Is this a statistical problem or more should I be looking more in the finance direction? For those with an finance background, the variables which show this pattern are market capitalization, price-to-book ratio and dividend yield. I've uploaded the FE regression result for the Australian market. Hopefully this works clarifying. Help is highly appreciated! PS. Lagrange and Hausman test favored FE model. Whole sample regression; ![](http://www.freeimagehosting.net/newuploads/pyjot.png) Sub period 1; ![](http://www.freeimagehosting.net/newuploads/d3ull.png) Sub perdiod 2; ![](http://www.freeimagehosting.net/newuploads/betwt.png)


title: Significance of Chi-square test for large sample

text: I am working on sample size of 70 . My observed and expected values are different (Observed 586, 648, 526, 662, 658, 502, ..... so on for 70 samples and expected 570 634, 513, 647, 644, 490, .....so on for 70 samples). I used Chi-sq test to check the null that there is no difference in observed and expected value. And I got a p-value of 0.99 which means I can't reject null. Can I use any other test for this kind of problem since I can see that for all the 70 samples observed is higher than expected but I cant show this result statistically.


title: Combination of Data from Two Normal Samples not Normal?

text: I have two sets of data that hypothesis tests have shown to be normal and from the same distribution. I'm using MATLAB and for the way they give p-values, higher p-values suggest a better goodness-of-fit, as I understand it. I want to combine the two data sets to get the parameters for the combined data set. However, when I run a hypothesis test on the combined data set, hypothesis tests either reject the null assumption of Normality or come very close as seen by the test statistics. I've been looking for an explanation but have no idea where to start. Could anyone shed some light into this?


title: comparing results of lmer using geom_area plots

text: I was wondering if somebody could help with trying to visualise some data. I have 11 groups of organisms and have related them to levels of a specific bacteria in water samples. I have a strong positive relationship between the proportion of one of the groups and some non-significant relationships in others, some positive, some negative. My idea was to show on one plot how the proportion of this positively related group increases as the level of bacteria increases and others stay similar or decline. I was trying to do this using the geom_area plot on ggplot http://docs.ggplot2.org/current/geom_area.html. Using the raw proportional data is a bad idea as the relationship is only really seen when including random effects in a mixed effects model. I considered using the fitted model values, however I am not sure if this is the best idea if I want to compare model results? Many thanks for any help with this.


title: The Conjugate Beta Prior proof

text: ![The Conjugate Beta Prior](http://i.stack.imgur.com/4dw9Z.jpg) Hello. I'm having a problem with trying to figure out this proof that shows the beta distribution is conjugate to the binomial distribution (picture attached). I understand it until the third row, but I got confused with this step from the third to the fourth row. I would appreciate a simple explanation. Thank you in advance.


title: Sample vs. Population

text: I'm trying to analyse the performance of my server by collecting performance data using standard tools such as Windows Performance Monitor, and then putting this data on Excel for analysis. When I set up the data collector (as it's called), I can chose different intervals for collection, such as every 1 second, every 10 seconds and so on. My first thought was that data collected every 1 second was going to be a "population", and greater intervals turned this data into "sample", but I soon realized that I could break it further into milliseconds and whatnot. Is the data I collect at any given interval going to be considered "population" or "sample"?


title: How to add standard error or confidence interval to a plot of a predicted vs observed values?

text: Is it possible to add standard error or confidence interval to a plot of a predicted vs observed values derived from a multiple regression model? I believe that I have seen such plots as an output in Statistica, but am unsure how to create them in R. I believe I have a solution (below), but am unsure that I have done this correctly. Basically, I have created a new `dataframe` with predictor variable in the range of their possible values. My worry with such an approach is that the prediction is based on the rows of data, and does not really address situations where the variables are randomly selected. Many thanks for your help. **Example:**               set.seed(1)     n <- 200     x1 <- rnorm(n, mean=10, sd=3)     x2 <- rnorm(n, mean=20, sd=5)     e <- rnorm(n, mean=10, sd=3)          y <- 5 + 2*x1 + 0.5*x2 + e          fit <- lm(y ~ x1 + x2)     summary(fit)          #plot of predicted vs observed     pred1 <- predict(fit, se.fit=TRUE)     plot(pred1$fit ~ y)     abline(0,1, col=8, lwd=2)          #new dataframe sequence of each predictor variable in their range     df.new <- data.frame(x1=seq(min(x1), max(x1),,100), x2=seq(min(x2), max(x2),,100))     pred2 <- predict(fit, df.new, se.fit=TRUE)          #plot of predicted vs observed w/ standard error interval?     png("pred_vs_obs.png", width=6, height=6, units="in", res=200)     plot(pred1$fit ~ y)         abline(0,1, col=8, lwd=2)         lines(pred2$fit+1.96*pred2$se.fit ~ pred2$fit, col=2, lty=2, lwd=2)     lines(pred2$fit-1.96*pred2$se.fit ~ pred2$fit, col=2, lty=2, lwd=2)     dev.off()      ![enter image description here](http://i.stack.imgur.com/AgUkO.png) **Edit:** The following code elaborates on my hesitation with the method that I used. The relationship between Standard Error (SE) and y is not a precise; i.e. various values of y that are relatively close together, have widely differing SE (black symbols in figure below, `pred1`), while the above method predicts a single SE for each predicted y (red symbols, `pred2`). Furthermore, using several different combinations of x1 and x2 that always result in the same y-value, I get a single (but different!) SE (green symbol, `pred3`). What is going on here? Is there a more correct way of doing this with some sort of permutation method?               #? Do different solutions to a given predicted value always give the same standard error?     y.tmp <- rep(40,20)     x1.tmp <- seq(0,10, length(y.tmp))     x2.tmp <- (y.tmp - fit$coeff[1] - fit$coeff[2]*x1.tmp) / fit$coeff[3]          df3 <- data.frame(x1=x1.tmp, x2=x2.tmp)     pred3 <- predict(fit, df3, se.fit=TRUE)          YLIM <- range(pred1$se.fit, pred2$se.fit, pred3$se.fit)         png("fit.se_vs_fit.png", width=6, height=6, units="in", res=200)         plot(pred1$se.fit ~ pred1$fit, ylim=YLIM, lwd=2)         points(pred2$se.fit ~ pred2$fit, col=2, lwd=2)         points(pred3$se.fit ~ pred3$fit, col=3, lwd=2)     legend("topright", legend=c("orig. data", "range of x1 & x2", "various comb. of x1 & x2 \nto acheive y=40"), col=1:3, pch=1, lwd=2, lty=0)     dev.off()      ![enter image description here](http://i.stack.imgur.com/4TcVs.png)


title: Comparing the value of two factors

text: I have made a factor analysis of one construct (i.e. difference between two countries) and got two distinct factors (1=cultural differences, 2=economical environment). I extracted the scores as I want to check if the score of cultural differences is bigger than economical. How do I do that? I have tried making a simple T test but it does not return any good result, as the factor scores are standardized and the factor means are 0 (S.Dev=1). What approach do I have to do?


title: What method(s) should I use to achive statistical ranking?

text: I'm attempting a project where I need to statistically rank available cars based on several variables such as cost, mpg, seating, milage, etc.. I wish to rank these cars in order decide which car would be the best choice (highest "worth") to buy (or best several cars if I was informing multiple people of the best cars to get). As the list of available cars changes from day to day, I will also need to re-run the code every day to allow the rankings to give me the best decision for this new day. What statistical methods should I use to go about this ranking system? I plan on determining which factors I find most important so the variables used will be subjective in choice. I thought about trying MDS or clustering but I didn't know if that would be relevant since I'm already subjectively determining what variables are to be used. I don't see how regression can be used since I can't get a handle on the "worth" of previous cars as that is what I'm trying to rank by. Also, I will be attempting this in R so any helpful packages/functions would be great to know as well. Any help with how to go about this ranking scheme would be helpful as I'm at a loss. Thanks so much


title: How to find out the relationship between variables?

text: I have two variables _Distance_ and _Slope_. These go into a model resulting in a _cost_ variable. I like to figure out the influence of the variables Distance and Slope on the Cost variable. I have several hundred test numbers for the variables. What programs or steps could I use to see the impact of the two variables on the cost? In the end I would like to know for example _Distance_ influences the _Cost_ variable too 70 %. _Slope_ influence _Cost_ too 30 %.


title: How do I compare means of the same group on two continous variables?

text: I would like to compare the means of the same group on two continuous scales. I want to see if the scores are significantly different on the scales. Which statistical test do I use for this?


title: Assigning a sample set to one of several populations based on composition

text: I have a problem that seems tractable, which I will attempt to describe with classic balls and buckets probability nomenclature. Imagine someone brings you a bag of colored balls. This sample was selected at random from one of many large buckets of balls. Each bucket contains 1,000,000 balls, and each ball is one of 1,000 colors. For any bucket, the proportion of balls for each color is fixed, and we know what it is. For example, Bucket A always has 0.02 Red, 0.013 Green, 0.009 Blue, etc. So--how would I go about estimating which bucket the sample of balls came from? I know that if I have a bag (sample) of 1,000,000 balls, I can say unambiguously which bucket that came from, as the proportions are defined. However, what do I do if I have only 2000 balls in the sample? It seems to me I should be able to calculate quite specifically the likelihood that the sample came from a particular bucket. Any help would be greatly appreciated! I'll be doing this in code (likely Python). Best regards, Dan


title: Understanding a split plot

text: Can someone explain the intuition behind the split plot. From what I understand, its is essentially restricted randomization. But I still quite understand it. Is there a resource or example anyone can give me to make it clearer?


title: How to interpret cumulative periodogram?

text: I use `cpgram` function in R to produce a cumulative periodogram of a monthly time series. Its horizental axis is a number between 0 and 6 labeled as frequency: what does these numbers mean, since the frequency in cumulative periodogram should be between 0 and 0.5 like the ones shown below (upper ones)? ![enter image description here](http://i.stack.imgur.com/Gs4O2.jpg) Also, do you know how can I create cumulative periodograms like these plots?


title: How does cross validation work in R's gbm package?

text: Can someone provide a work flow about this? For instance, suppose I am doing binary classification, For each iteration of the algorithm: Randomly sample k*N rows, where k is the bag.fraction, and N is nrow(dataset). Build a classifier using this training sample, suppose we use classification tree here. Predict the terminal node class label. This how boosting is done without a CV. If I want to do a 3-fold CV, where exactly do I apply it? Between step 1 and 2 or after step 3? Thanks!


title: Generalizing van Elteren to more than two treatments

text: I know about the following tests:   * Mantel-Haenszel for categorial stratified data   * Kruskal-Wallis for nonparametric ANOVA for more than two treatments   * Friedman for a repeated measures version of Kruskal-Wallis   * van Elteren for nonparametric comparisons of stratified data between two treatments But is there something that combines the features of these? In other words, I'm looking for a test that takes data for $g$ experimental treatments, collected across $q$ strata, where the response cannot be assumed to follow any parametric distribution, and I want to test for a difference in treatment effect while controlling for the confounder(s) that I stratified across. The sample size is small. If there isn't such a test (I could not find any references in the literature), then can we resort to some kind of nonparametric bootstrap, or a Monte Carlo method? How could we design such a procedure?


title: Using similarity matrix to measure diversity of a group

text: I wish to measure the "diversity" of a group of objects. Right now I'm using Euclidean distance to compute the similarity matrix between all the objects in the group. I'm searching for a measure of this diversity to decide if replacing an element of this matrix with a new one will make the group more or less homogeneous. I wish to end up the most heterogeneous group. Of course, I'm open to another approaches if I'm looking at this from the wrong angle. Bonus points if it's computationally efficient.


title: Tobit-like model for dependent variable with upper limits

text: Following the wiki page and the form of the likelihood function for a generalized Tobit model presented here, I am thinking about applying a model of this type to a problem I'm facing but need some help in the interpretation. I have a variable $Y$ which I assume to linearly depend on $X$ via a parameter $\beta$. My model is \begin{equation} Y = X\beta + \epsilon_i \end{equation} The problem is that some observed values I have for $Y$, the $y_i$'s, are upper limits. Can I extend the concept of "censored data" to this case and use the Tobit model to perform the linear regression? Do you have any thoughts on the possibility of further inclunding a categorical moderator variable in the regression? Is this possible? Thanks.


title: Orthogonal series density estimation

text: I am going through this paper Orthogonal series density estimation. I have a doubt in following Assume that the random variable X is supported on [0, 1], that is, P(X ∈ [0, 1]) = 1, and that the probability density f of X is square integrable. Then the density may be approximated with any wished accuracy by a partial sum $f_J(x) :=\sum_{j=0}^J \theta_j \varphi_j(x), 0\le x\le1$ Here $\{\varphi_j\} $ is the cosine basis $ \varphi_0(x)=1, \varphi_j(x)=\sqrt 2 cos(\pi j x), j=\{1,2,\cdots\}$ $\theta_j$ is the $j_{th}$ coefficient and $J$ is the cut-off Note that $\theta_0=\int_0^1f(x) dx=1$ Also $\theta_j=\int_0^1 \varphi_j(x) f(x) dx= E\{\varphi_j(x)\}$ This implies that $\theta_j$ can be estimated by sample mean estimator. $\hat{\theta}=n^{-1}\sum_{l=1}^n \varphi(X_l)$ $Var(\hat{\theta_j})=E(\hat{\theta}-\theta)^2$ $=n^{-1}[1+2^{-1/2}\theta_{2j}-\theta_j^2]=n^{-1}d$ Using $cos^2(\alpha)=(1+cos(2\alpha))/2$ Here d is defined as coefficient of difficulty. My doubt is in last step $E(\hat{\theta}-\theta)^2=n^{-1}[1+2^{-1/2}\theta_{2j}-\theta_j^2]$ I think $E({\hat{\theta_j}^2}-2\hat{\theta_j}{\theta}+{\theta}^2)=E({\hat{\theta_j}^2}-{\theta}^2)$ If we assume the points are drawn iid $E(\varphi(x_l),\varphi(x_m))=E(\varphi(x_l))E(\varphi(x_m))$ So $E((\sum_{l=1}^n\varphi(x_l))^2)=E(\sum_{l=1}^n\varphi(x_l)^2+2*\sum_{l,m=0}^n\varphi(x_l)\varphi(x_m))$ $=E(\sum_{l=1}^n\varphi(x_l)^2)+2*E(\sum_{l,m=0}^n\varphi(x_l)\varphi(x_m))$ $E({\hat{\theta_j}^2})=n^{-2}[n(1+2^{-1/2}\theta_{2j})+\big(n(n-1)/2\big)2*{\theta_j}^2]$ $E({\hat{\theta_j}^2})=n^{-1}[1+2^{-1/2}\theta_{2j}+(n-1){\theta_j}^2]$ Please let me know if my understanding is correct


title: Building the connection between cosine similarity and correlation in R

text: According to some articles (e.g. here) correlation is just a centered version of cosine similarity. I use the following code to calculate the cosine similarity matrix of the column vectors of a matrix `X` (slightly modified from here):               cos.sim <- function(ix)      {       A = X[,ix[1]]       B = X[,ix[2]]       return(t(A)%*%B/sqrt(sum(A^2)*sum(B^2)))     }        n <- ncol(X)      cmb <- expand.grid(i=1:n, j=1:n)      C <- matrix(apply(cmb,1,cos.sim),n,n)      **My question**   Which modifications of the code above are needed to get the correlation matrix `cor(X)` instead of the cosine similarity matrix. I guess the changes are minimal but I can't see them at the moment.


title: Metric optimization on discrete learning sample

text: There are a set of ("artifical") not Minkovski (triangle inequality is not guaranteed) metrics defined on set of objects. There are one etalon ("natural") metric, which estimation is known only for some objects (a finite learning subset). How can I get some combination of "artifical" metrics (new metric essencually) that would be the best (for proposed class of combinations) approximation of the "natural" metric? I mean the type of appropriate methods or something specifying the area of further googling. May be it's just two dimensional function approximation where both dimensions present the same discrete space as a distance matrix and function is an etalon metric. Specific realization of distances combination is not the point (though lower computational complexity is better), so it may be simple linear combination. And kind of objects is not important too, but in my case it's some data sets for which I define several similarity measures that must be combined into one corresponding to the expert estimation.


title: Forecasting daily electricity price

text: I'm trying to make a model to forecast the electricity price, a time series model with R and i have some questions Our data are daily price of the past 3 years from north pool countries, and we are trying to do a model to forecast the daily prices too. We don't know how to descompose the seasonality of our data, because all the test and studies need to be made at the data without seasonality. In the other hand, we have some problems to understand the ARIMA or ARMA model, because is the one we want to apply, and we don't know how to do it in R. Thanks and sorry for my bad english.


title: Monte Carlo for DOE Power Calculation When Performance is Unknown

text: I have been wanting to get away from using power calculation tools (JMP, Lenth) in favor of Monte Carlo techniques to determine the a-priori power to observe factor level effects in a DOE. I am really new to this and have a few questions. Lets start with what I have:   1. A binomial response which I plan to perform a logistic regression to   2. List of factors that I think are relevant   3. A model (Main Effects + interaction + quadratic terms)   4. An experimental design (Taken from JMP or the DoE.wrapper package)   5. A threshold requirement (0.6)   6. An effect size (ES = 0.1) I would like to determine the power to observe an effect in the response for each one of the model terms. I have read the excellent response by user gung which I wanted to use to help tailor my own approach. After reading, I had a few specific questions.   1. The example uses prior knowledge of the system performance (the rates) that I do not have. How do I continue without this knowledge?    2. The solution does not factor in effect size in any way (at least as far as I can tell). How do I factor in effect size to this process?   3. The power for each of the terms is coming from the significance of the model term to be non-zero. Is this the same as the power to observe an effect size in the response? For question 2, I would assume that you would just need to run the model assuming the null performance (p = 0.6) and then re-run the fit assuming the performance of Null + Effect Size (p = 0.6 + ES). I know this last part might be a bit jumbled up. Any help would be greatly appreciated.


title: Self Study Question on Regression

text: Good evening everyone, I am working on the following self-test question for regressions, but my answer 26.901 seemed really different from the model answer 27.2155. **_Question_** > A hotel derives its gross income from its hotel and restaurant operations. > The owners are interested in the relationship between the number of rooms > occupied on a nightly basis and the revenue per day in the restaurant . > Below are a sample of 25 days (Monday through Thursday) from last year > showing the restaurant income Y in thousand dollars and number of rooms > occupied X > > ![enter image description here](http://i.stack.imgur.com/NvVcQ.jpg) **_My Attempt_** > ![enter image description here](http://i.stack.imgur.com/pjueE.jpg) Appreciate any advice please.