title: Calculate relative price levels in 4 markets?

text: I am working in Excel and want to calculate the relative price level in 4 markets. I know a good method I can use when I am only looking at 2 markets                            USA    UK     Bubble gum    4     6     Lollypop      3      Chocolate bar 7     8      In the above example with 2 markets I would conclude that the price level is (6/4+7/8)/2 = 1,1875 higher in the UK than it is in the USA. We can express this as USA = 1 UK = 1,1875 But how would I calculate this when we have four markets and don't have prices for all items in all markets?                            USA    UK    France    Germany     Bubble gum    4     6       8     Lollypop      3             5          4     Chocolate bar 7     8      10      I need your help calculating the price levels in the example above. I'd like to do the calculations in Excel if possible.


title: Two-stage least squares approach to Deming regression

text: I am interested in statistical inference for the Deming regression model: $$ x_i=x^*_i + \epsilon_i$$ $$ y_i = (\alpha+\beta x^*_i) + \epsilon'_i$$ where the $x^*_i$'s are nonrandom fixed numbers, $\epsilon_i \sim {\cal N}(0,\sigma_x^2)$, $\epsilon'_i \sim {\cal N}(0,\sigma_y^2)$, and all variables $\epsilon_i, \epsilon'_i$ are mutually independent. I particularly focus on the case of the assumption $\sigma_x=\sigma_y$. Graham Dunn's slides show how to treat the Deming regression with a two-stage least squares approach when the $x^*_i$'s are random (this model is sometimes termed as the ultrastructural model in the literature). I would like to know how to similarly treat the case when the $x^*_i$'s are nonrandom and whether it is possible to perform this approach with John Fox's "sem" R package.


title: Is it worth to take another sample?

text: There is a random normal variable, with unknown mean and std. I want to estimate the mean by sampling. After several samples, I have an estimate, which I can make more precise by taking more samples, however each additional sample has a cost. How can I estimate, from the samples I currently have, the effect on an additional sample on the accuracy of my estimate of the mean?


title: ML estimation of a multinomial proportion with constraints

text: **The Question** A researcher surveys $n$ college students and counts how many support, how many oppose and how many are undecided about a recently introduced federal policy. Letting $X_1, X_2, X_3$ denote these counts, she model $X = (X_1, X_2, X_3)$ as $X \sim \mathsf{Multinomial}(n,p)$. I want to find the restricted MLE of $\mathbf p = (p_1, p_2, p_3)$ for a null hypothesis that the actual proportions of supporters and opposers in the entire college are equal. My mathematical formulation of this null was that this MLE would be the vector $\bf p$ that maximizes the likelihood function over $\Delta_3^0 \subsetneq \Delta_3$ where $\Delta_3$ is the 3-dimensional simplex of all probability vectors i.e. vectors that sum componentwise to 1. Then $\Delta_3^0$ would contain all $p$ of the form $\mathbf p = (a, a, 1 - 2a), a\in(0,1)$. _Note_ I realize that this a simple question for such mathematical abstractions, I have simply tried to follow my teacher's notation. **My Work** When I went through and tried to write out the likelihood function, I got $$L_{x, H_0} = \text{const. }\times a^{x_1}\cdot a^{x_2}\cdot (1 - 2a)^{x_3}$$ which should maximize at $\hat a = {x_1 + x_2 \over x_1 + x_2 + x_3}$. However, when I actually calculated for data of counts $(X_1, X_2, X_3) = (140, 165, 195)$ I got $a = 0.61$, which gives me a negative probability in my probability vector which doesn't make sense. Thanks for any help.


title: Terminology of two different tabular representations of aggregated data

text: think of a group of people aged from 21 to 35 years. Then we can show the number of persons regarding gender and age groups 21-25, 26-30 and 31-35 in (at least) two different ways: **Type A** ![enter image description here](http://i.stack.imgur.com/gLjAX.jpg) **Type B** ![enter image description here](http://i.stack.imgur.com/nWum5.jpg) What I am curious about is **what names are used to refer to those different types of table structures**? The practical relevance is that Excel can use Type A as a source for Pivot Tables but not Type B. Often though huge tables of type B are provided for statistics and I would like to google for 'How to turn Type B tables into Type A tables?'. * * * **Terminology (summarized from answers)** :   1. Type A:     * **R** : 'long data format'     * '(melted) cross tab     * **Excel** : 'flat(tened) data table'   2. Type B:     * **R** : 'wide data format'     * (traditional) cross tab The search phrase for Excel/Pivot table would be 'how to flatten a cross table'. * * *


title: Generating a population of individuals with known Gaussian and uniform random properties

text: Suppose that I have a number of variables. Each is known to be Gaussian or uniformly random with known parameters and occur with a known probability. I also have a table of correlations (or covariances if that's easier) for each one. It would be great if these parameters could come from any kind of distribution, but if not, uniform and Gaussian are fine for my interests. So, a typical problem case would be A is Gaussian with mean a1 and variance a2; B is uniform with maximum b1 and minimum b2. C is Gaussian with mean c1 and variance c2. Then I have a table of correlations A with B, A with C and B with C. How would I write an algorithm to generate a sample of random vectors (A,B,C) which satisfy the above properties. Again, it would be really great if the variables could come from any kind of distribution, but I would settle for at least Gaussian only, uniform only or a mixture of Gaussian and uniform variables. It would be nice to have an approach that would scale up easily. I would like to implement this myself, rather than use some predefined package. So, I'm interested in algorithmic details.


title: How to analyse data collected in different ways

text: Let's say I collect some quantitative case-control data. First thing that pops into my mind is using a t-test to see if these two groups differ. Let's say some more data is collected, but it's just a bit different (a different experiment collected it, or machine that collects data has some variability or whatever). What are my options on analyzing this data? First thing that pops into my mind is regression or ANOVA (using indicator variables) as a comment points out below : $E[Y] = \beta_0 + \beta_1*CASE + \beta_2*BATCH$ where $BATCH$ represents in which 'batch' the data was collected (or which experimenter, or machine, etc). Are there better ways to do this analysis? I know this is quite an open-ended question but I just wanted to see what else is available.


title: Use a baseline including all observations in lm()

text: When using factors in a linear model, I would like to retrieve the 'average' effect as baseline (intercept), rather than level 1 of all factors. Is this possible?               tab <- data.frame(gl(2,5),rnorm(10)+as.numeric(gl(2,5)))     colnames(tab) <- c('pred','dat')     mylm <- lm(dat~pred,data=tab)     coefficients(mylm)      returns something like               (Intercept)       pred2        0.5879712   1.8542153       coefficients are on average 1 for both intercept and pred2. I would like to rather have intercept = 1.5, pred1 = -0.5 and pred2 = +0.5               model.matrix(mylm)      returns                        (Intercept) pred2         1            1     0         2            1     0         3            1     0         4            1     0         5            1     0         6            1     1         7            1     1         8            1     1         9            1     1         10           1     1         attr(,"assign")         [1] 0 1         attr(,"contrasts")         attr(,"contrasts")$pred         [1] "contr.treatment"      I would like a matrix like this                         (Intercept)  pred1   pred2         1            1    1     0         2            1    1     0         3            1    1     0         4            1    1     0         5            1    1     0         6            1    0     1         7            1    0     1         8            1    0     1         9            1    0     1         10           1    0     1      Does that make sense? if not, why? Is there some command to generate such a design matrix? I've been looking for a way to do this without success...


title: Recursive feature elimination with only two classes

text: Recursive feature elimination (RFE) is a feature-selection strategy. It performs in two nested levels of cross-validation. First it tries to divide the training set into N folds. RFE puts one fold aside for testing the generalization and then trains itself with the remaining data. In my case I only have two classes and would like to do RFE. The problem here is that, since there are only two classes, the algorithm can not put one aside and then train itself with the remaining folds (only one fold remains and training on single class is not sensible.) I would be very much grateful if you could let me know what I could do here.


title: Intuition behind standard deviation

text: I'm trying to gain a better intuitive understanding of standard deviation. From what I understand it is representative of the average of the differences of a set of observations in a data set from the mean of that data set. However it is NOT actually equal to the averages of the differences as it gives more weight to observations further from the mean. Say I have the following population of values - $\\{1, 3, 5, 7, 9\\}$ The mean is $5$. If I take a measure of spread based on absolute value I get $$\frac{\sum_{i = 1}^5|x_i - \mu|}{5} = 2.4$$ If I take a measure of spread based using standard deviation I get $$\sqrt{\frac{\sum_{i = 1}^5(x_i - \mu)^2}{5}} = 2.83$$ The result using standard deviation is larger, as expected, because of the extra weight it gives to values further from the mean. But if I was just told that I was dealing with a population with a mean of $5$ and a standard deviation of $2.83$ how would I infer that the population was comprised of values something like the $\\{1, 3, 5, 7, 9\\}$? It just seems that the figure of $2.83$ is very arbitrary...I don't see how you are supposed to interpret it. Does $2.83$ mean the values are spread very wide or are they all tightly clustered around the mean... When you are presented with a statement that you are dealing with a population with a mean of $5$ and a standard deviation of $2.83$ what does that tell you about the population?


title: Cross fitting with same params but differents models

text: What is the best fitting way with 2 variables to explain ($Z_1$ and $Z_2$) depending on the same variables ($X$ and $Y$) and parameters $\theta$ but with differents models ($f$ and $g$)? For instance: $Z_1 = f(X,Y,\theta) $ and $Z_2 = g(X,Y,\theta)$. I have an estimation of $\hat{X}$, $\hat{Y}$ and $\hat{\theta}$, and two mechanistic models $f(.)$ and $g(.)$. I would like two regressions in the same time. I used a GLM, GLMM and `nls()` function in `R` which work to estimate $\hat{\theta}$ for $Z_1$ and $Z_2$ one at a time. Are Bayesian model the only way? Where may I look for?


title: How to impute data without missing at random?

text: Recently I got a global longitudinal data from several countries, and each county has one outcome variable and two predictors from 1995 to 2008. I found one of the predictors is always missing in each country in 1995, 1997, 1999, and 2001 because that variable was collected every two years before 2001. This situation seriously affects another complete predictor when using a complete case analysis because too many useful information is lost in that four years. Also, it obviously affects the fitting of a time smoothing function. I am pretty wondering whether this situation is appropriate to use the multiple imputation method to generate data. As I know, this case is not a missing at random mechanism, so the multiple imputation may not be a perfect way to deal with my case. I am looking for any advice here to find out a solution to deal with those missing data. Any suggestion is appreciated.


title: Finding the Expected Average Distance from the Maximum given a distribution

text: For a given sample set $S$ with $N$ individual samples $x_i$, I can easily find the average distance from the maximum by doing something like this: $\sigma_{max_N}:=\sqrt{\frac{1}{N}\sum\limits_{i=1}^N {\left(x_i - \max{\left(S\right)}\right)}^2}$ Trying to do so for larger and larger sample sizes, I (unsurprisingly) end up with different limit values given different distributions. For instance, for a uniform distribution on a unit interval, I get a value around 0.577 with digits after that varying wildly. A normal distribution with the same mean $\left(\frac{1}{2}\right)$ and variance $\left(\frac{1}{12}\right)$ gives me something around 1.5, though the exact value shifts around a lot more. For an exponential distribution with parameter $\lambda=2\sqrt{3}$, I end up at around 4.5. Is there a way to get an expected value of $\sigma_{max_N} \to \sigma_{max}$ for a given distribution as the sample size $N \to \infty$? It doesn't necessarily have to be analytic.


title: How to decide and calculate the skewness of an asymmetric confidence interval around a mean or a mean difference?

text: I wonder how statistical programs decide to skew the confidence interval? (I mean the confidence intervals for the means, or for the mean differences). How do they decide on the extent of this (without any bootsrapping)? I know they do it based on the skewness of the sample distribution (in two-tailed tests), but would appreciate any details on that.


title: Is output of Deamer deconvolution not a density?

text: I have a Model Y= X+e and need the density of X. The deamer package deconvolves the density for X, but if I use the simpsons rule to integrate this density, I get values which are above 1. The following example gives me a density which integral is 1.173454:               library(deamer) # deconvolution     library(Bolstad) # simpson rule          # The Y's I have are inv-Weibull distributed and the error's are inv-normal distributed.     # As the deconvolution of those would take a long time, i used runif in my example     # to simplify the problem. The following uncommented lines are what I like to deconvolve:          #library(actuar) # for rinvweibull     #y <- rinvweibull(30000, shape=5.53861156, scale=488)/1000     #y <- y[y<1.5]     #e <- 1/rnorm(30000, mean=0.0023853421, sd=0.0004784688)/1000     #e <- e[e<1.5]     #decon <- deamerSE(y, error=e, from=-0.1, to=0.3)          y <- runif(1000, min = 0.8, max=1.2)     e <- runif(1000, min = 0.1, max=0.5)     decon <- deamerSE(y, error=e, from=0.4, to=1)     plot(decon)          # following line gives me integral of density (with simpsons rule)     sintegral(decon$supp, decon$f)$value      I am not sure if this is just an estimation error or if I should consider downscaling the density so that the integrated density is 1:               # Downscaling     yValsScaled <- decon$f/area         plot(decon$supp, yValsScaled, type="l")     (areaScaled <- sintegral(decon$supp, yValsScaled)$value)      What do you think? btw: If you use a larger intervall with the _from_ and _to_ argments in _deamerSE_ function, the integrated density will be even greater (because of the periodicity of the density). Usually I thought that with deamerSE I would get a density which integral (from -inf to inf) is approximately 1. Therefore I thougt that integrating the density using a smaller intervall (e.g. with _from=0.4_ and _to=1_ in the deamerSE function) should give me a density which integral is less than 1. But as you can see it doesn't. So I am rather confused.


title: Confidence-intervals for conditions tested with a mixed-effects model

text: Sorry for a possibly ignorant question. I have fit a mixed-effects model using the lmer function from the lme4 package, and the main fixed effect (a factor with three levels) in the model was significant according to a run with pvals.fnc (from the languageR package). To illustrate the effect in an appealing way, I would like to plot three bars with the "baseline condition" (the intercept including the effect level 1), the "experimental condition 1" (the intercept + the effect level 2), the "experimental condition 2" (intercept + effect level 3), and the corresponding confidence intervals for these conditions. However, how do I do that? The function plotLMER.fnc used to do that, but has stopped working for me (lme4_0.999375-32; languageR_1.0).   Model: mylmer <\- lmer(outcome ~ (1|subject) + (1|item) + Factor, data) So I have MCMC output in the form of:               $fixed                   Estimate MCMCmean HPD95lower HPD95upper  pMCMC Pr(>|t|)     (Intercept)     0.4728   0.4718     0.2250     0.7368 0.0010   0.0008     Factor2        -0.0420  -0.0420    -0.1732     0.0931 0.5414   0.5710     Factor3        -0.1643  -0.1631    -0.3153    -0.0112 0.0328   0.0508      My idea was to construct the conditions from this data by simply using the intercept as the baseline condition, and construct the experimental conditions by adding the effect of each factor level to the intercept. The new (intercept + effect combined) CI would be constructed by turning the HPD interval into a standard deviation and then use the square root of the sum of squared standard deviations (like here) Is this approach appropriate? If I do it, the CIs become become fairly large and the effect no longer looks significant (conflicting with the pvals.fnc output). Is there a better way (code example would be great)? If not, is there a "least worst" solution that could satisfy readers who strongly request the standard bar plot with CIs? Problems in the back of my head, which I am too ignorant to properly formulate:   CI not equal to HPD.   Transforming CI to SD is inappropriate, as CI may be asymmetric.   Problems with determining CI at all in mixed models (like here).


title: Compare rank orders of population members across different variables

text: I have some data which describe residential units for people with learning disability, variables like how nice the furnishings are, the level of psychiatric symptomology on the unit, how happy the staff are, stuff like that. I want to check to see if we are measuring the right things- e.g. do units with happier staff have nicer atmospheres, do units with nice furnishings have happy staff, that kind of thing. The problem is that I only have data for 8 units (averages within each, e.g. 10 staff give their responses on how happy they are, this is the average for the unit) so I can't really use linear regression to see whether the things that we have measured affect each other. I've drawn some scatterplots for all of the data and on the whole I would say there does appear to be a linear relationship in the way I would expect. But as I say with 8 units it's not really going to generate any statistics. I had the bright idea of ordering each unit across all variables, and then comparing the rank orders somehow. If we are measuring the right things then the ranks should be similar, like this:               Unit 1 (ranks across all variables): 1,1,1,1,1,1,1,1,1     Unit 2: 2,2,2,2,2,2,2,2,2      etc. whereas if I'm wrong, and the variables aren't important to each other, I will get this:               Unit 1: 1,2,3,4,5,6,7,8     Unit 2: 8,7,6,5,4,3,2,1     Unit 3: 4,5,6,7,8,1,2,3      etc. Here's what I get:               Unit1    7    5  5.0    5    3    4    5    3     Unit2    6    2  4.0    6    5    3    2    5     Unit3    3    7  7.5    1    4    1    1    1     Unit4    4    4  3.0    7    6    7    7    8     Unit5    5    3  1.0    4    2    5    6    7     Unit6    2    6  6.0    8    8    8    8    6     Unit7    1    8  7.5    3    7    6    4    4     Unit8    8    1  2.0    2    1    2    3    2      It looks pretty good to me, except for the first column. Any thoughts on this? Is is a proper statistic that I haven't heard of? Or is there something reasonably robust I can do with these results? Sorry for the lengthy question, many thanks in advance!


title: Distribution of norms of vectors whose components have a known distribution?

text: I know very little of statistics, so forgive me if the question is badly posed. Suppose I have $n$ samples $x_1,\dotsc,x_n$ drawn from some distribution $f$ (let's say $n\approx 100$). Fix $1\leq k \leq n$ (in practice $k \ll n$, say $k\approx 5$). Consider the set $$ X = \\{\lVert (x_{i_1}, x_{i_2}, \dotsc, x_{i_k}) \rVert \quad \vert \quad 1\leq i_j \leq n \text{ (distinct)}\\}, $$ where $\lVert \cdot\rVert$ is some $L^p$-norm. Can I say anything of the distribution of the members of $X$ (in terms of $n$, $k$, $f$ and $p$, of course)?


title: Creating ROC curve for multi-level logistic regression model in R

text: I used the functions from this link for creating ROC curve for logistic regression model. Since the object produced by `glmer` in `lme4` package is a S4 object (as far as I know) and the function from the link cannot handle it. I wonder if there are similar functions for creating ROC curve for multi-level logistic regression model in R.


title: Identification of parameters problem

text: I always struggle to get the true essence of identification in econometrics. I know that we state that a parameter (say $\hat{\theta}$) can be identified if by simply looking at its (joint) distribution we can infer the value of the parameter. In a simple case of $y=b_1X+u$, where $E[u]=0,E[u|x]=0$ we can state that $b_1$ is identified if we know that its variance $Var(\hat{b})>0$. But what if $E[u|X]=a$ where $a$ is an unknown parameter? Can $a$ and $b_1$ be identified? If I expand the model to $Y=b_0+b_1X+b_2XD=u$ where $D\in\\{0,1\\}$ and $E[u|X,D]=0$, to show that $b_1,b_2,b_3$are identified, do I simply have to restate that the variance for all three parameters is greater than zero? I appreciate all the help on clearing my mind concerning identification.


title: Stacked bar plot

text: I have a data-frame whose first column is the name of an item and the second column is the frequency of that item in the dataset.                names            freq     1 tomato           7     2 potato           4     3 cabbage          5     4 sukuma-wiki      8     5 terere           20      I would like to have a stacked bar column that depicts the proportion of each entry on the chart. How do you handle coloring of the stacked bar when presented with over sixty entries? what is the easiest way to do this?


title: Maximizing probability of winning on loaded coin

text: If we have a loaded coin that plays out 75% heads, 25% tails, what would be the best way to bet on the outcome of each of $n$ trials? How could we maximize our probability of winning? Is it possible to generalize for a coin that's loaded $n:(100-n)$?


title: Hit and run MCMC

text: I'm trying to implement the hit and run MCMC algorithm, but I'm having a bit of trouble understanding how to go about it. The general idea, is as follows: > To generate a proposal jump in MH, we: > >   * Generate a direction $d$ from a distribution on the surface of the unit > sphere $\mathcal{O}$ >   * Generate a signed distance $\lambda$ along the constrained space. > However, I have no idea of how I should go about implementing this in R (or any other language). Does anyone have a snippet of code that would point me in the right direction? BTW, I not that interested in a library that does this method, I want to try and code it up myself. Many thanks.


title: How do you estimate correlations between ordinal and binary data in SPSS?

text: I'm struggling with my lab report. I have a questionnaire which contains Likert responses (1-7) as well as questions that require only yes or no answers. I want to see if there is a correlation between higher scores on the likert scale and responses to the binary answers to the other questions. Specifically - the likert response is about anxiety levels. The binary (yes/no) answer is about negative consequences due to alcohol use. I want to just run the one ordinal likert item with all the negative consequence responses to see if there is a relationship. Can anyone help me? * * * Thanks so much for the fast comments. It definitely assures me that Spearman's Rho can be used but my main problem is the fact that i am mixing binary and ordinal variables together in the one correlation matrix. I have one ordinal variable (anxiety - 1-7 on likert scale) crossed with 13 dichotomous (0 or 1 representing yes or no) to answers ("have you ever been hospitalised due to your alcohol use?). I get statistically significant answers (<.05 and < .01 but the Correlation coefficients are tiny (i.e o.2) at best. the hospitalisation is the DV an the anxiety levels the IV. I just think i'm doing something wrong here. should i collapse the 13 binary variables in to one some how? i can't see how SPSS will recognise the values between 1-7 when correlating with binary answers. i keep thinking it will just read that there is something in that box (ordinal) but not recognise the value (i.e 1-7) . but people keep saying that SPSS does in fact recognise this difference and account for it.


title: How to calculate the degree of freedom for Williams' Test

text: Williams' test is a multiple comparison test. My question is, when compare the Williams' test statistic with the tabulated critical values, which degree of freedom should I look for. In the SAS example, a substance has been tested at seven levels($k=6$) in a randomized block design of eight blocks and the mean square is with (7 - 1)(8 - 1) = 42 degrees of freedom. From my understanding of degree of freedom, shouldn't it be $7\times8-7$ since the $s$ used in the test statistic as an estimate of $\sigma$ is the pooled variance? More formally, The null hypothesis is: $$ H_0: \mu_0=\mu_1=...=\mu_k $$ The alternative is: $$ H_1: \mu_0\leq\mu_1\leq...\leq\mu_k $$ with at least one strict inequality. The MLE of the treatment means $\mu_i$ are obtained by the sample means $\bar{x}_i$ and the number of observations in the samples $n_i$, by the formular $$ \mu_i^*=\max_{l\leq u\leq i}\min_{i\leq v \leq k}\sum_{j=u}^{v}n_j \bar{x}_j/\sum_{j=u}^{v}n_j $$ The test statistic is: $$ T_i=\frac{\mu_i^*-\bar{x}_0}{\sqrt{s^2/n_i+s^2/n_0}}$$ where, $s$ is an unbiased estimated of $\sigma$, the within group standard deviation that is independent from $\bar X_i$. The null hypothesis is rejected and the fact that the $i$-th dose level is the minimum effective dose is concluded if $$T_j> t_{j,\alpha}, \mbox{for all } j \geq i$$ where, $t_{j,\alpha}$ is the upper $\alpha$th percentile of the distribution $T_j$, with $\nu$ degree of freedom. My question is how to calculate $\nu$? And why? I thought $\nu$ should be $\sum n_i -(k+1)$, but my simulation results agree more with the SAS way of using $n_0\times k$, where all $n_i$'s are equal. Another point, for all $j$'s , I should use the same degree of freedom since I calculated $s$ based on all samples not after removing the already effective high dose level samples. Am I right?


title: Averaging across independent strata: confidence interval

text: I've $k$ independent strata, with populations $n_1,\dots,n_k$, and have estimated the mean in each stratum $\bar{x_1}, \dots \bar{x_k}$ and obtained the standard errors $se_1,\dots,se_k$ (using STATA mean command). I'm interested in an estimator of the weighted mean and its standard error, using weights $w_1,\dots,w_k$ with $\sum_i w_i=W$. To compute the weighted mean I simply use: $$\bar{x}_w=\frac{1}{W}\sum_{i=1}^k w_i \bar{x_i}$$ Next, I want to estimate the standard error of the weighted mean by using the (estimated) standard errors of each strata, the weights and $n_1,\dots,n_k$. I imagine I can use the fact that $V(\alpha x + \beta y)=\alpha^2V(x)+\beta^2 V(y)$ and back-up the desired estimator as: $$se(\bar{x}_w) = \frac{1}{W}\sqrt{\sum_{i=1}^k w_i^2 se_i^2}$$ Am I right? I'm a bit confused about the different formulas used by STATA and your help is appreciated.


title: Interpretation of interaction effects in linear mixed models with numeric and factorial IVs

text: I hope my question is not too basic for this community. I can*t figure out how to interpret the output of my linear mixed model, especially the interaction effects. I do my analysis in R using lme from the nlme package. My model consists of 3 factorial IVs (exp.type [2 levels - field vs. lab], org.type [3 levles - autotroph, heterotroph, mixed], system [2 levels - marine vs. freshwater]) and 2 numeric IVs (even, duration). > lme1 = lme(reist ~ even + exp.type + even:exp.type + org.type + > even:org.type + system + duration, random =~1|authors.year, data = data) ![output of the linear model](http://i.stack.imgur.com/ZhcZK.jpg) My question is: How do I interpret the interaction effects? How do I get the slopes of the linear regressions for each level of the IVs that have significant interactions


title: Cluster large boolean dataset

text: I got a dataset with about 5,000 columns and about 135,000 rows - all fields contain boolean (binary) data. I am looking to classify each of these columns into one of 50 groups, based on similarity. Can anyone point me in the right direction where to start? Will a k-mean style algorithm suit me? Any library in `R` that I can easily access? Note: data is quite asymmetric in that about 95% of values are 0, 5% are 1.


Let $(X_1,Y_1), . . . , (X_n,Y_n)$ be a random sample from the discrete distribution with joint probability mass function $$ f_{X,Y} (x,y) = \frac{\theta}4 , \space (x,y) = (0,0)\space or \space(1,1)$$ $$ = \frac{2 - \theta}4, \space \space (x,y) = (0,1)\space or \space(1,0),$$ with $0 \le \theta \le 2$. Find the maximum likelihood estimator of $\theta$.


title: Normalizing rating in a group of people [finding effectiveness]

text: I'm a new guy here. Hopefully, I'm asking this question to right forum. ## Problem: We have data of a group of people (P1, P2, P3). They rank their expertise (1-10, where higher number is better) in a list of components (G1, G2, G3).                   P1  P2  P3     --------------     G1 | 8   4   7     G2 | 7   3   7     G3 | 9   6   5      Also, we have some data regarding work done by each person in each component. Example: For P1,                    W  WD     ----------     G1 | 0   0     G2 | 2   0     G3 | 8   2      where `W` is total work allotted to user P1, and `WD` is actual work done. W >= WD >= 0. We have similar data for P2, and P3 users. **Point to note:** The user might have some level of expertise regardless of work done in a component. Example: P1 has ranked himself 8/10 even though he has not been given any task in G1 component (W = 0). P1 also has ranked himself 7/10 for G2 even though he has not finished any task in that component (WD = 0). Now, we want to calculate effective rating of all users relative to the group of users, not self-ranking, considering their work data and self-ranking. Can anyone suggest some mechanism to achieve this? Thanks much in advance!


title: Multiple imputation questions for multiple regression in SPSS

text: I am currently running a multiple regression model using imputed data and have a few questions. **Background:** Using SPSS 18. My data appears to be MAR. Listwise deletion of cases leaves me with only 92 cases, multiple imputation leaves 153 cases for analysis. All assumptions met - one variable log transformed. 9 IV's 5 - 5 categorical, 3 scale, 1 interval. DV-scale. Using the enter method of standard multiple regression.   * My DV is the difference of scores between a pre- score and a post score measure, both of these variables are missing a number of cases - should I impute missing values for each of these and then work out the differnce between them to calculate my DV (how do I go about doing this), or can I just impute data for my DV? Which is the most appropriate approach?   * Should I run imputations on transformed data or skewed untransformed data?   * Should I enter all variables into the imputation process, even if they are not missing data, or should I just impute data for the variables missing more than 10% of cases? I have run the regression on the listwise deleted cases and my IV's account for very little of the variance in my DV, subsequently I have run the regression on a complete file following multiple imputation - The results are very similar, in that my 9 IV's still predict only approx 12% of the variance in my DV, however, now one of my IV'S indicates that it is making a significant contribution (this happens to be a log transformed variable)...   * Should I report original data if there is little difference between my conclusions - i.e my IV's poorly predict the dv, or report the complete data?


title: Logit with ordinal independent variables

text: In a logit model, is there a smarter way to determine the effect of an independent ordinal variable than to use dummy variables for each level?


title: Assumptions for Spearman's rho test for correlation

text: If I were to use Spearman's rho as a measure of correlation, must the data I'm working with be normally distributed/homoscedastic? Also, if I am plotting the data against time, must the data be time- independent?


title: Correlation coefficient contradicts t-test

text: I am looking at the probability of opting into a program based on a continuous variables $X_1$, $X_2$, $X_3$, etc. When I divide the sample into people who opted in and did not opt in and do a t-test I find that the mean of all the $X$'s for those opting in is significantly different from the mean for those not opting in. However, when I look at a correlation matrix for the $X$'s and opting in I find really low correlation coefficients -- from between $.04$ and $.13$. But my t-tests are significant at the 95%+ confidence level. How do I square these two results?


title: Determining if a function is additive

text: I have data collected on a grid of values $f(x,y)$. I have a hypothesis that the data is "additively separable", i.e. $$ f(x,y) = h_1(x) + h_2(y) + g(x,y) $$ where $g(x,y)$ is small and only represent the noise in my measurements. Is there a way to test for this? If I had infinite sampling resolution I could measure $$ \frac{\partial^2 f(x,y) }{\partial{x}\partial{y}} = \frac{\partial^2 h_1(x) }{\partial{x}\partial{y}} + \frac{\partial^2 h_2(y) }{\partial{x}\partial{y}} + \frac{\partial^2 g(x,y) }{\partial{x}\partial{y}} = g_{xy}(x,y) $$ Then I could test if for some threshold value $\epsilon$, $|g_{xy}(x,y)| > \epsilon$. However, my data is real, with noise and limited samples. What statistical test could I use to check for additive separability?


title: Why are Winsorized random variables independent?

text: While studying trimmed mean I understood that if I have some random variables $X_1, X_2, .., X_n$ by ordering them and trimming, the variables are no longer independent. However it is said that "by re-writing the trimmed mean as the average of independent variables" one can re-use the standard estimation of variance. This re-writing I know that proposes Winsorizing the values. The trimming consists of eliminating a given percentile of values from the sample. The Winsorized sample is the one which has values less then a given percentile censored, and to an upper percentile censored also. It is clear that the trimmed mean is not the same as the Winsorized mean. But how can one conclude formally that the variables after Winsorizing the variables are independent? [Later edit] The question arise during the study of the book "Fundamentals of Modern Statistical Methods, Substantially Improving Power and Accuracy, 2nd edition - Rand R. Wilcox". The following quotes are from pages: 149-151. > [my note: as an introduction we have a $X_1,X_2,..,X_n$ i.i.d. standard > normal variables, they are trimmed and a trimmed mean with $20\%$ is built. > It follows a discussion on estimating variance of trimmed mean.] > > Consequently, the method we used to determine an expression for the variance > of the sample mean is an average of dependent variables. (In fact, there is > yet one more complication: The variances of the ordered variables are not > equal to $\sigma^2$. > > A major theoretical advance during the 1960s was the derivation of a > mathematical technique that leads to a convenient and practical method for > estimating the variance of the sample trimmed mean. The method consists of > rewriting the trimmed mean as the average of independent variables so that > the strategy used to determine the variance of the sample mean can again be > employed. The mathematical details are not given here; [..] > > [my note: next follows an description of how to build an Winsorized sample > variance and concludes with the estimand] $$ VAR(\bar{X_t}) = > \frac{s^2_w}{(1-2\gamma)^2n}$$ [my note: where $X_t$ is trimmed variables, > $s^2_w$ is Windsorized sample mean, $\gamma$ is the trimmed quantity ($0.2$ > in this case), $n$ number of sample values]


title: neural network always have the same output

text: I am using the 'neuralnet' package in R. The example on http://gekkoquant.com/2012/05/26/neural-networks-with-r-simple-example performs good. But when I try the following, I find for the training data, the output are always the same and the training error is huge. Also I find that if setting the dimension of Input as 1, then sum(Output)=sum(Input). This is very strange. I GOOGLE 'neural network always same output', there are several pages concerning the same problem, but no answer. I also coded the BP algorithm in MATLAB, and find the situation also exists. Here is the code in both R and MATLAB R code:               install.packages('neuralnet')     library("neuralnet")          install.packages('mvtnorm')     library("mvtnorm")          #Going to create a neural network to perform sqare rooting     #Type ?neuralnet for more information on the neuralnet library          #Generate 50 random input vectors, dimension of each being 10     #And store them as a dataframe          dimension<-10     data <- rmvnorm(n=50, mean=seq(1:dimension), sigma=diag(dimension), method="chol")     traininginput <-  as.data.frame(data)          #Output equals the norm of first half of the vector     trainingoutput <- apply(traininginput[,1:(dimension/2)]^2,1,sum)          #Column bind the data into one variable     trainingdata <- cbind(trainingoutput,traininginput)          #Column names as 'Output', 'X1', ..., 'X10'         colnames(trainingdata) <- c('Output',paste(rep('X',dimension),seq(1:dimension),sep=''))          #Train the neural network     #Threshold is a numeric value specifying the threshold for the partial     #derivatives of the error function as stopping criteria.          model1.form.text <- paste("Output ~",paste(colnames(trainingdata)[-1],collapse=" + "),collapse=" ")     model1.form <- as.formula( model1.form.text )          net.sqrt <- neuralnet(model1.form, trainingdata, hidden=c(4), threshold=1e-2)          #Output of the network for each input ARE ALWAYS THE SAME     #I tried several kinds of hidden layers, such as c(4,2), c(5,4,3,1) etc, the situation always happen     print(net.sqrt$net.result)          print(net.sqrt)          #Plot the neural network     plot(net.sqrt)      MATLAB code (written according to http://www.amlbook.com/slides/iTunesU_Lecture10_May_03.pdf for the definition of w, delta, s and x)               %%This code tries to realize the function of a neural network          %%The activation function is set as tanh. dactivate is the derivative of tanh     activate=@(s) (exp(s)-exp(-s))./(exp(s)+exp(-s));     dactivate=@(s) 1.-activate(s).^2;               %%n: number of observations; dimension: dimension of inputs;     %%three layers, each with neurons 7,3,1     n=50; dimension=5; layer=[7 3 1]';     X=randn(n,dimension); y=sum(X,2).^2;          %%add 1s to the last column as the bias     X=[X,ones(size(X,1),1)];              %%w: weight of each layer     w1=1e-5*rand(dimension+1,layer(1));     w2=1e-5*rand(layer(1)+1,layer(2));     w3=1e-5*rand(layer(2)+1,layer(3));          %%wlast is the weight of last transform, so that the output ini_slast=[ini_x3, ones(size(ini_x3,1),1)]*ini_wlast;          wlast=1e-5*rand(2,1);%% [w0 w1]';          %%iteration of forward-backward algorithm, set learning rate as 1e-3. Using BATCH learning     lrate=1e-3;     for i=1:1000          %%set initial values     ini_w1=w1;     ini_w2=w2;     ini_w3=w3;     ini_wlast=wlast;          %%s_i=X*w_i, bias already added into X and w_i     ini_s1=X*ini_w1;     ini_x1=activate(ini_s1); ini_s2=[ini_x1, ones(size(ini_x1,1),1)]*ini_w2;     ini_x2=activate(ini_s2); ini_s3=[ini_x2, ones(size(ini_x2,1),1)]*ini_w3;     ini_x3=activate(ini_s3);     ini_slast=[ini_x3, ones(size(ini_x3,1),1)]*ini_wlast;          %%calculate delta     delta3=2*(ini_slast-y) *ini_wlast(1) .*dactivate(ini_s3);         delta2=bsxfun(@times,dactivate(ini_s2), delta3*ini_w3(1:end-1,:)');         delta1=bsxfun(@times,dactivate(ini_s1), delta2*ini_w2(1:end-1,:)');          %%calculate the gradient gw     gwlast=2*(ini_slast-y)'*[ini_x3, ones(size(ini_x3,1),1)]/n ;     gw3=delta3'*[ini_x2, ones(size(ini_x2,1),1)]/n;  %%size(x2)=[1,?]     gw2=delta2'*[ini_x1, ones(size(ini_x1,1),1)]/n;     gw1=delta1'*X/n; %%+2*lambda*ini_w1';          %%update the weights     w1=ini_w1-lrate*gw1';     w2=ini_w2-lrate*gw2';     w3=ini_w3-lrate*gw3';     wlast=ini_wlast-lrate*gwlast';     %%pause;     display(norm(y-ini_slast));          end          display([y ini_slast]);     sum(([y ini_slast]))


title: Experiment Design - Randomized Block or Completely Randomized

text: We have Five litters of rats, each containing four rats. The five litters were randomly selected from a population of litters. Investigators are interested in comparing the mercury content from fish digested by rats, as well as the overall performance of measurement techniques in identify mercury volumes. Within each of the five litters, the four rats were assigned a species of fish (there are 4 types of fish). So, what rat per species, randomly assigned. After a month, a random selection of the rat's kidney (randomly choose left or right) was chosen, and the kidney was split into three sections. Each section measured the mercury content using one of the 3 measurement techniques. The technique used was randomly assigned to ea. 1/3 of the kidney. My Take: I beleive this to be a Complete Block Design with Split-Plot. Block: Litter My question here, if we randomly choose litters from a pouplation... is it really a block or a CRD? I would only think it was a CRD if a TREATMENT was applied at this stage.. but it is not. The treatment is applied at the next experimental unit level on the rats. Rat: Experimental unit for Fish Species (Whole-Plot EU) Kidney: Random Factor (2 Levels) 1/3 Kidney: Experimental Unit for Measurement Technique (Split-Plot EU) So, Litter = Block, Rat = EU for Fish (Rat is also a nested factor within Litter), Kidney is a Random Factor of Rat (Not Nested, simply "left or right"), and 1/3 of Kidney = EU for Measurement (Split-Plot) and is also a factor with 3 fixed levels, [B]not [/B] nested within kidney (simply 1, 2, or 3 piece). Thoughts?


title: Binary Classifier evaluation when Precision is more important than Recall

text: **Problem statement** A set of points is given. We want to classify those points in two distinct classes with labels $\\{0,1\\}$. We count the "hits" of the classifier based on the $1$ class. Assume that we have some classifiers in hand with the following properties   * For some classifiers, the classification decision depends on some thresholds, which we may vary and plot a Precision/Recall or RoC curve   * For others, the decision does not depend on any threshold but rather on a Property test. That is, we check if point $x_i$ has some property, and classify it accordingly. For illustration purposes one can imagine such a test as testing for example if $x_i$ is odd, and classifying all odd data points in the same class.(However this is not the problem I am working on!) Additionally, we care more about misclassifying Label-$0$ points as Label-$1$ points (errors that would affect Precision), than misclassifying Label-$1$ points as Label-$0$ points(errors that affect Recall). To put it simply, we want the classified $1$ set to contain few outliers(Label-$0$ points) and we can sacrifice, to some extent, the size of the retrieved $1$-set. **Question** Which measures could I use to evaluate the performance of each classifier and pick the classifier $c_i$ or the classifier/threshold pair $(c_i,T_i)$ that best solves the problem. Assume we have a dataset with known labels on which we evaluate the classifiers.


title: Deconvolution with fourier transform or characteristic function?

text: Let us consider the following model: $$Y_j = X_j + \epsilon_j \hspace{15pt} j=1, ..., n$$ Where $Y_j$ is a noisy signal, $\epsilon_j$ is the noise which is independend from the signal $X_j$. We have only i.i.d. samples of $Y_j$ and $\epsilon_j$ and are interested in the distribution of $X_j$. The density $f_{\epsilon}$ is assumed to be unknown. Comte and Lacour suggest a method based on fourier transform to solve this problem (see section 2.2). Let's call $\varphi_X(t)$ the characteristic function for $X$ and $f^{*}_X$ the fourier transform of the density $f_X$. Here is my question: The main idea in deconvolution is to use the independence assumption for $X$ and $\epsilon$ and then use fourier transform to solve the equation $f^{*}_X = f^{*}_Y/f^{*}_{\epsilon}$. Applying the inverse fourier transform leads to an estimate for $f_X$. Can I use the characteristic function instead of the fourier transform? Does this give me any advantages or disadvantages? I assume that both fourier transform or characteristic function could be used but would like to know what other people think.


title: Where does $\sqrt{n}$ come from in central limit theorem (CLT)?

text: A very simple version of central limited theorem as below $$ \sqrt{n}\bigg(\bigg(\frac{1}{n}\sum_{i=1}^n X_i\bigg) - \mu\bigg)\ \xrightarrow{d}\ \mathcal{N}(0,\;\sigma^2) $$ which is Lindeberg–Lévy CLT. I do not understand why there is a $\sqrt{n}$ on the left handside. And Lyapunov CLT says $$ \frac{1}{s_n} \sum_{i=1}^{n} (X_i - \mu_i) \ \xrightarrow{d}\ \mathcal{N}(0,\;1) $$ but why not $\sqrt{s_n}$? Would anyone tell me what are these factors, such $\sqrt{n}$ and $\frac{1}{s_n}$? how do we get them in the theorem?


title: Does anyone know any good open source software for visualizing data from database?

text: Recently I came across Tableau and tried to visualize the data from database and csv file. The user iterface enables the user to visualize time and spatial data and create plots in an instant. Such tool is really useful as it enables to observe the data graphically without writing the code. As there are many data sources from which I have to retrieve and visualize the data it would be very useful to have a tool which enabled to generate charts by simply dragging columns on axes and additionally modify the visualization with dragging the column names as well. Does anyone know any free or open source software of that kind?


title: Improving accuracy of total variance through variance components?

text: I'm following a measurement guidance that specifies a technique to improve the estimate of total variance by extracting individual variance components and performing a sum of squares. I'm not a statistician, so I'm looking for somebody to validate the methodology defined by the guidance. This is best explained by an example: We sell a product that needs to meet particular performance requirements and is composed of two separate pieces which may potentially be interchanged outside of our control. As such, we cannot test every possible combination. To show that the product meets the requirements, we perform a study that samples the performance of pairings in a matrix, with two measurements on each pair to capture measurement error. Here's a data set for this: ![Data matrix](http://i.stack.imgur.com/doQwC.png) The guidance states that the sample variance under-estimates the total variance and the solution is to extract the individual variance components and then combine them separately. The calculations, _as prescribed by the guidance_ , are below: ![Calculations on groups](http://i.stack.imgur.com/12zZq.png) ![Calculations on the total](http://i.stack.imgur.com/EM8YK.png) Most of this technique actually makes a lot of sense to me. I see the value in applying the components separately to account for the different sample sizes, but there are also a few aspects of this that I don't fully understand. In particular:   1. The calculation of the measurement variance (Eq 3) is just taking the average of the variance at each pairing. Is this the proper technique to extract measurement variance?   2. Why do the component variances (Eqs 4 & 5) subtract out the variance in the measurement mean? This is particularly troublesome where the component variance is less than the subtracted variance in the measurement mean (as with Eq 5). Should we just set that component to zero in those cases?   3. In the calculation of the total variance of the mean (Eq 7), shouldn't the divisor for $S^2_{meas}$ actually be pq, not pqr? Aside from this guidance, I've been unable to find any other examples that use this exact method. The data setup looks like a two-way crossed ANOVA but the analysis is different. For reference, the guidance refers to the book by Poduri S.R.S Rau, "Variance Components Estimation: Mixed Models, Methodologies, and Applications" [Chapman & Hall, 1997, pp 12-38] as the source of the technique. I'd appreciate it if anybody can validate this technique and help me understand the questions I posed. Thank you for your assistance. PS: Sorry about the use of images instead of formulas - I had already created the images for another board where this question has gone unanswered.


title: Empirical coverage probabilities of confidence intervals

text: Can someone tell me how the type of distribution is related to an empirical coverage probability? (for some distributions the Central-limit-theorem seems to be more easily satisfied)


title: Get Characteristic Function from MGF?

text: I've been calculating characteristic functions and MGF's and was wondering whether we can always get the characteristic function simply by substituting $it$ instead of $t$ in the resulting equation. For example the normal MGF is $ \exp{(\mu t + \frac{\sigma^2 t^2}{2})} $ but if we substitute $it$ into that we get $ \exp{(\mu it - \frac{\sigma^2 t^2}{2})} $ which is the characteristic function. Can we always do this if the MGF exists in a neighborhood of 0?


title: How to analyze repeated A/B choice questions?

text: I feel a bit silly, but I'm stuck on what I thought would be an easy analysis. Haven't had this type of data in a while, so I'm a bit rusty. Tried Googling is, but as I'm not a native English speaker I'm not entirely sure I'm using the right terms.. Anyway, so much for the apology in advance. Participants in my study answered 24 questions where they simply had to choose between option A and option B. The questions are all "the same" in the sense that H0 is that they will answers 50%/50%; 12 times A and 12 times B. How do I aggregate and analyze?   1. Do I aggregate the data into a percentage A (or B, for that matter) per person and then run a one-sample t-test against 50%?   2. Do I aggregate into counts of A (again, or B of course) and then somehow do a Chi-square? <\- this is where it gets fuzzy!


title: Probability of detection and false alarm ratio are wrong in a simple example, don't know why

text: I'm trying to understand how POD (Probability of Detection) and FAR (False Alarm Ratio) are calculated in this paper. On page 234 they are defined as usually and the tree data structure in page 239 shows POD and FAR values for a simple example but I can't figure out how they came up with those values. $$ POD = \frac{TP}{(TP+FN)} \\\ FAR = \frac{FP}{(TP+FP)} $$For this example, let:   * TP = There's tornado (+ label in the series) and there's the word in the series    * FP = There's no tornado and there's the word in the series    * TN = There's no tornado and the word is not in the series    * FN = There's tornado and there's not word in the series Then for the first word (aa), I understand that:   TP = 2; FN=0; FP=2; so POD=1 and FAR= 0.5 which is what is in the paper, But for the second word (ab):   TP= 2; FN=0; FP=1; so POD=1, which is in the paper, but shouldn't FAR= 0.333, instead of 1/4 in the paper? And for the third word (ac):   TP= 1; FN=1; FP=1; so POD=1/2, which is in the paper, but shouldn't FAR= 1/2, instead of 1/4 in the paper? I don't understand what's wrong with my FAR, and if my POD are OK, as they seem to be.


title: Should train and test datasets have similar variance?

text: If variance of test dataset is lower than the one of the train dataset is it worth splitting the data? Since we know our dataset will always be limited is it fair to select models under the above condition? Thanks


title: Find residual sum of squares

text: Let $Y_1, Y_2,Y_3$ be uncorrelated random variables with common variance $\sigma^2>0$ such that   $$ \begin{aligned} E[Y_1]&=\beta_1+\beta_2, \\\ E[Y_2]&=2\beta_1 \\\ E[Y_3]&=\beta_1-\beta_2 \end{aligned} $$ where $\beta_1$ and $\beta_2$ are unknown parameters. Find the residual (error) sum of squares under the above linear model. I proceed in this way: Note that: $$ \begin{aligned} E[Y_1-\beta_1-\beta_2]&=0, \\\ V[Y_1-\beta_1-\beta_2]&=\sigma^2, \\\ E[Y_2-2\beta_1]&=0, \\\ V[Y_2-2\beta_1]&=\sigma^2, \\\ E[Y_3-\beta_1+\beta_2]&=0, \\\ V[Y_3-\beta_1+\beta_2]&=\sigma^2 \end{aligned} $$ So the residual(error) sum of squares under the above linear model is: $$ \sum{e_i^2}=(Y_1-\beta_1-\beta_2)^2+(Y_2-2\beta_1)^2+(Y_3-\beta_1+\beta_2)^2 $$ Am I proceeding in the correct way? Can someone please explain what is meant by linear model here?


title: Сonfidence interval of histogram probability density function estimator

text: There is a circle of radius $R$ and a sequence of points within it. I'm going to estimate a PDF of appearing at elementary area at the distance $D$ from the centre of the circle using a realization of $N$ points $X_0$, $X_1$, $\dots$, $X_{N-1}$. So I divide interval $[0, R]$ into $n$ bins $b_i=[r_i,r_{i+1}]$ and find $N_i$ as the number of points of $D$ appearing in each bin and divide it into the area of the ring of each bin $S=\pi\left(r_{i+1}^2-r_{i}^2\right)$. (For some reason bins might have different lengths $r_{i+1}-r_{i}$.) I'm going to find error bars for values of the histogram $N_i/S$, or better say the function given by a table of values $(r_i+r_{i+1})/2\longrightarrow N_i/S$ as a PDF estimator accurate within a factor of $N$. They usually use something like the "standard Poisson error" as $\sigma=1/\sqrt{n}$, so am I to use $\sigma_i=(1/\sqrt{N_i})/S$ or $\sigma_i = 1/\sqrt{N_i/S}$? Or is there a more correct formula to find a confidence interval? Also, how should I compute a confidence interval if I convert the histogram into logarithmic coordinates $\log((r_i+r_{i+1})/2)\longrightarrow \log (N_i/S)$?


title: Intuition and uses for coefficient of variation

text: I'm currently attending the _An Introduction to Operations Management_ course in Coursera.org. At some point in the course, the professor started to deal with variation in the operations' time. The measurement he uses is the _Coefficient of Variation_ , the ratio between the standard deviation and the mean: $c_v = \frac{\sigma}{\mu}$ Why would this measurement be used? What are the advantages and disadvantages of working with _CV_ besides working with, say, standard deviation? What is the intuition behind this measurement?


title: Logistic regression: big difference in predicted values and highly significant, but poor goodness of fit

text: I have a logistic regression model with a dichotomous response variable and predictors coded from $1$ to $10$ and from $0$ to $18$. When I fit the model, I get these results:               Intercept     2.467 (p=2e-16)     Predictor#1  -0.181 (p=1.76e-07)     Predictor#2  -0.098 (p=3.34e-14)          Null deviance 2252.3 on 1741 DF     Residual dev. 2113.1 on 1739 DF     (1276 obs. deleted due to missingness) — AIC: 2119.1      And a Nagelkerke $R^2$ of $0.10$. The $R^2$ and the Deviance show the model as a very bad one, but when I calculate the predicted values I get a possibility of $24.8\%$ when both predictors are $1$, and of $89.9\%$ when the predictors are at their maximum values ($10$ and $18$, respectively). How can I get such bad measures of the goodness of fit, and, at the same time, get such a big difference (statistical significance, $p<0.01$) in the predicted values?


title: Are there packages for fitting ordinal logistic/probit mixed models with random slopes in R?

text: I'm looking for a way to fit an ordinal logistic and/or probit mixed model that includes random slopes. The only package for R I could find that allows for random effects at all in ordinal mixed models are the **ordinal** and **ordinal2** packages, but only fitting of random intercept models is currently implemented in those packages.


title: Free data set for comparing kernel logistic regression and regular logistic regression

text: I'm looking for a data set that is easily accessible for comparing Kernel Logistic Regression (KLR) and regular logistic regression. All the paper that I find using KLR use synthetic data sets. I'm most interested in data sets with continuous variables, so I can apply the squared exponential kernel. Data sets with a binary target and spatial predictor variables would also be interesting.


title: Evaluation method when using a large training set and a small test set

text: I am facing the evaluation of two text classifiers. I have a **large training dataset** (to be used for training only), **and a separated small test set** (to be used for testing only), **both being balanced**. Which one of the following methods should be most appropiated one and why? 1) Stratified Repeated held-out evaluation (repeated subsampling): Sample k times without repetition from the training dataset, each sample being balanced. For every sample, classifier is trained with the sample, and accuracy is tested with the full test dataset. Results are averaged. 2) Stratified Cross-validation: Divide the full training dataset in k slices with equal size, each slice being balanced. For every slice, classifier is trained with the slice, and accuracy is tested with the full test dataset. Results are averaged. 3) Stratified Bootstrapping: Sample k times with repetition from the large training dataset, each sample being balanced. For every sample, classifier is trained with the sample, and accuracy is tested with the full test dataset. Results are averaged.


title: Hypothesis testing for Pareto distributions

text: I wish to to some simple hypothesis testing of the form provided by T-Tests and ANOVA. However, my data is not normally distributed (it follows a Pareto distribution). My understanding is that T-Tests make the assumption that the data is normally distributed and hence I won't be able to use them - is that correct? Is there something else I can do? **EDIT** Here is some more info about my problem. I'm trying to do some quality analysis on software defects, and am having trouble knowing where to start. One basic question I want to answer is: > Does software produced in department X have more defects than department Y? As some background, we group changes to software as "patches", in which case the question becomes > Does the average patch from department X have more defects than department > Y? Here is a histogram of bugs / patch, N = 3700. ![bugs/patch](http://i.stack.imgur.com/h1Vju.png) There is a philosophical issue of what it means to have "more" defects that I don't have a great answer to. The obvious choice for one of my limited knowledge is to compare the mean defects in each group, but as others have pointed out that's not clearly the best choice. The measure linked to by Procrastinator ($P(X<Y)$) seems like it captures my intuition well.


title: How to perform a repeated measure CFA

text: My purpose is to verify that a specific scale is unidimensional, and that the one dimensional structure is persistent over time. I have a 10 item scale, about 200 cases and score values at two different time points. It would be simple if it was only one measure: I would split the data in half, I would perform PCA in one half and CFA in the other half, an approach that indeed is working well for the first measure. However, I would like to somehow provide a measure that the unique factor structure is persistent over time. **How can this be done?**


title: Bayesian hypothesis testing and Bayes factors

text: Supposing a parameter of interest is assumed to follow an exponential distribution with rate $\lambda$. Two competing prior hypotheses are formed such that the first prior is given by $\text{Ga}(\alpha, \beta)$ and the second prior is given by $\text{Ga}(\alpha+1, \beta)$. Furthermore, data is collected on the parameter such that $x_{1}, x_{2}, \dots, x_{n}$ represent $n$ observations of recorded data. The individual responsible for forming these priors suggests a degree of belief such that $p(\lambda \sim \text{Ga}(\alpha, \beta)) = 0.75$ and $p(\lambda \sim \text{Ga}(\alpha+1, \beta)) = 0.25$. From the above:   * How would a Bayes factor for the two prior distributions be formed? Can this be achieved by calculating the ratio between the two prior distributions? i.e. $$\text{Bayes} factor = \frac{p(x) \sim \text{Ga}(\alpha, \beta)}{p(x) \sim \text{Ga}(\alpha+1, \beta)} = \frac{\frac{\beta^{\alpha}}{\Gamma(\alpha)}x^{\alpha-1}e^{-\beta x}}{\frac{\beta^{\alpha+1}}{\Gamma(\alpha+1)}x^{\alpha}e^{-\beta x}}$$ If so, given that there are multiple observations ($x_{i}$), how should each observation be inserted into the distribution to give a composite distribution (is it the sum of $p(x) \sim \text{Ga}(\alpha, \beta)$ such that $\sum_{i=1}^{n} p(x) \sim \text{Ga}(\alpha, \beta)$)? In which case the Bayes factor is given by: $$\frac{\sum_{i=1}^{n} p(x_{i}) \sim \text{Ga}(\alpha, \beta)}{\sum_{i=1}^{n} p(x_{i}) \sim \text{Ga}(\alpha+1, \beta)}$$   * How can an estimate of the observed data be combined with the model probabilities to achieve a posterior ratio of model probabilities? Can this be achieved by getting the product of the prior and likelihood distributions? $$p(\lambda|{\bf x}) = p(x|\lambda)p(\lambda)$$ If so, assuming there is only a collection of observations and no suggested likelihood, how can the correct likelihood be determined? This a generic homework-like template question, as opposed to a specific exercise hence the absence of specific parameter and data values.


title: How to account for a regressand affecting a regressor?

text: I forget the terminology, but this happens when you regress, say, $Y$ on a list of variables, and you suspect that $Y$ affects, say, $x_3$ in addition to $x_3$ affecting $Y.$ I forget how this is called, and how to deal with it -- something along the lines of adding an extra variable? Thanks!


title: Are random forests Bayesian?

text: Hopefully my naivety doesn't shine through in this question: Are random forests Bayesian?


title: Dirichlet posterior

text: I have a question about the Dirichlet posterior distribution. Given a multinomial likelihood function it's known that the posterior is $Dir({\alpha_i + N_i})$, where $N_i$ is the number of times we've seen $i^{th}$ observation. What happens if we start to decrease $\alpha$s for a given fixed data $D$? It seems from the form of the posterior that after some point $\alpha$s will stop affecting posterior at all. But wouldn't it be right to say that when we make $\alpha$s very small the probability mass moves to the corners of the simplex and the posterior must be affected to a greater extent? What statement is the correct one?


title: Can Inferential Statistics be carried out in a Retrospective Case Study?

text: Just wanted to check with regard to a retrospective case study, if various data is collected for a set of patients on different variables can this be interpreted with inferential stats tests? Even if no power calculation could be made to determine the required sample size before the study began i.e. just the patients that were available with a certain pathology were used? Many thanks for any thoughts.


title: Cox regression with time varying covariates

text: I am looking for resources on Cox proportional hazards model with time varying covariates. I'm new to survival analysis so I'm looking for something not overly mathematical. I would also appreciate any information about software implementations that can deal with this problem.


title: What are complete sufficient statistics?

text: I have some trouble understanding complete sufficient statistics? Let $T=\Sigma x_i$ be a sufficient statistic. If $E[g(T)]=0$ with probability 1, for some function $g$, then it is a complete sufficient statistic. But what does this mean? I've seen examples of uniform and Bernoulli (page 6 http://amath.colorado.edu/courses/4520/2011fall/HandOuts/umvue.pdf), but it's not intuitive, I got more confused seeing the integration. Could someone explain in a simple and intuitive way?


title: Analysing wind data with R

text: Hi i am analaysing wind data for estimating energy from a wind turbine.   I have taken 10 years of wind data and graphed a histogram;   my second stage was to fit a Weibull distribution to the data.   I used R with the package `lmom` to compute the Weibul shape and scale this is the code i used:               >library(lmom)         wind.moments<-samlmu(as.numeric(pp$WS))           moments<-pelwei(wind.moments)          x.wei<-rweibull(n=length(pp$WS), shape=moments["delta"], scale=moments["beta"])         hist(as.numeric(pp$WS), freq=FALSE)         lines(density(x.wei), col="red", lwd=4)          It seems like there is some lag between the data and the density function; can you help me with this? Another question is can you help me in calculating the anual energy from the density function? ![enter image description here](http://i.stack.imgur.com/AqVij.png)   thank you


title: Location test under a bounded non-stationarity?

text: Suppose I draw $n$ observations $X_1,X_2,\ldots,X_n$ independently from a distribution where $X_i \sim \mathcal{N}(\mu_i,\sigma^2)$, where the mean is assumed to be Lipschitz: $\left| \mu_i - \mu_{i+1}\right| \le \gamma,$ with $\gamma$ known. I want to test the null hypothesis: $$H_0: \mu_n = 0$$ against the local alternative $H_1: \mu_n > 0$. A few questions:   1. Is this a known problem? If so, for what terms should I google search?   2. What is a reasonable way to perform this test? I can imagine there might be some way to do it by downweighting less recent observations--maybe throw the innovations into the noise term, then use a weighted t-test? The problem with that approach is that the change in mean may be deterministic, so I should pose the problem as a minimax problem maybe. Perhaps another approach would be to ignore the non-stationarity and use a regular t-test on some subset of the observations $X_j,X_{j+1},\ldots,X_n$, with $n-j$ chosen by some function of $\gamma, \sigma$ (probably their ratio?), and accept that the test will not maintain the nominal Type I rate. What are some better ideas? **Edit** : per the excellent comments so far, I believe this question can be reformulated as follows: Let $V$ be the convex set in $\mathbb{R}^n$ defined by: $$ V = \left\\{\mathbf{x} | x_n = 0, x_{i-1} \le x_i + \gamma, x_{i-1} \ge x_i - \gamma, i = 2,\ldots,n \right\\} $$ Observing vector $X = f + \sigma\epsilon$, where the elements of $\epsilon$ are drawn independently from a standard Gaussian, test the null hypothesis: $$H_0: f \in V$$ In this sense, it is similar to the paper of Baraud _et. al._ mentioned by @Robin Girard, except the Baraud paper provides a method for the case where $V$ is a linear subspace of $\mathbb{R}^n$, _not_ a convex polytope. Under this formulation, it seems that, under the null, the squared distance from $X$ to $f$ should be distributed as a Chi-squared, up to the scaling factor of $\sigma^2$. Then the distance from $X$ to $V$ should be no greater than that. The proposed test would then be to find the projection of $X$ onto $V$, and test the distance from $X$ to that projection based on the Chi- squared bound. I imagine this test will have fairly low power. Perhaps truncating the observation to a final subset of the $X_i$ would have higher power...


title: Choosing an appropriate model - What does this phrase mean?

text: I have a question where I am given a short story: > An investigator wishes to evaluate the effect of medicine administered for > colds on reaction time. Data represent the reaction times (in seconds) on a > stimulus one hour after each of four randomly chosen brands of this type of > medicine was administered to 5 people. They then ask me to: _Give an appropriate model_ After spending some time searching for what this means, I am mostly getting results for "choosing the right test". What I am thinking of choosing is between:   * Analysis of Variance   * Covariance Analysis   * Regression Analysis Can these be considered "models" or is a "model" defined as something else (in statistical terms) ?


title: Can the AICc (corrected Akaike information criterion) measure be used for post-hoc analysis?

text: I am collecting a lot of data from neurons and will be analyzing each neuron separately. Is it valid to pick the best of several models based on the AICc _for each neuron_ and then use the ANOVA p-values from the best model? Or should should I only run the AICc on preliminary data and then use the best model analyze every neuron? My impression was that it is only statistically valid to select the model in the preliminary experimental design phase, but I'm now having difficulty finding this explicitly stated. References would be greatly appreciated.


title: If a factor variable is to be dropped in model selection, should all levels be dropped simultaneously? If so, why?

text: In answer to a previous question factor pooling in model selection was discussed. If a factor or categorical variable is to be dropped in model selection, should all levels be dropped simultaneously? If so, why? The motivation for dropping factors is to aid model interpretation. For example, I might be interested in explaining the factors that influence customer behaviour when visiting a store and have a categorical variable "travel mode" with factors "walking, bus, private car, taxi, etc." In this context, I can remove all the dummy variables except "private car" because they have an insubstantial estimated magnitude and are not significant predictors of behaviour. I then end up with a "travelled in private car" vs "didn't travel in private car" variable and don't have to worry about troubling the reader with interpreting the other largely uninteresting variables.


title: Searching for time series inside another time series

text: I have a time series "A" and another one "B". I would like to find occurrences of "B" inside "A". Typically, "A" is much bigger (magnitude: millions of points) than "B" (magnitude: hundreds of points)   * I don't necessarily want exact matching (obviously, since the values we deal are real numbers not, say, text) The values are preferably close, though.   * Trends/shapes are more important than exact time matching (i.e. somehow Dynamic Time Warping could be useful)   * With regards to the first two points, what happens if I don't care about the absolute values at all but just the shape (maybe taking a log of the values?)   * What happens if I would like to scale the number of points to billions from millions? How to scale the system? How do I index a time series? (maybe use Mahout, Lucene, somehow?) What's inside above parentheses are the results of my short research, not necessarilly correct. All in all, what kind of an approach should I take? What are your general suggestions/tips? I'd like to develop a system, so any practical examples/suggestions are even more welcome. Thanks!


title: One factor repeated measures where first time point is a control

text: Patients were randomized to two groups (control and intervention). Each patient is evaluated at 3 time points. The first time point is taken pre- intervention, so it should have the same mean response in both intervention and control groups. The idea is to test for different slopes in the two groups. I can do this with a repeated measures anova, testing for the group x time interaction. I am wondering if a stronger test would involved constraining the mean effects for time 1 to be equal, but I am unsure how to set this up (either in R or SPSS). $$y_{gkti}=\mu+\alpha_g+\tau_{gt}+\eta_{k}+\epsilon_{gkti}$$ subject to $\tau_{11}=\tau_{21}$ Where the $\alpha$ are the group effects for groups 1 and 2, and the $\tau$ are the time effects and the $\eta$ are the patient effects. Any ideas?


title: Beta coefficients from stratified analysis when there are covariates?

text: Suppose I have a regression model shown below **Model 1:** $$ Y = \beta_0^\ + \beta_1SEX\ + \beta_2ALCOHOL\ + \beta_3SEX*ALCOHOL\ $$ The predictors I am interested in are SEX (binary: 0 female, 1 male) and Alcohol consumption (binary: drinker, non-drinker). Suppose that I found a significant interaction between SEX and ALCOHOL and decided to stratify the data by sex. So I would have two new models: **Model 2a:** $$ \text{Female: } Y_F = \beta_0^\ + \beta_2ALCOHOL\ $$ So for the female subset, the intercept is still $ \beta_0$ and the slope for ALCOHOL is $ \beta_2$ **Model 2b:** $$ \text{Male: } Y_M = (\beta_0^\ + \beta_1) + (\beta_2\ +\beta_3) ALCOHOL $$ For the male subset, the intercept is now $\beta_0^\ + \beta_1$ and the slope for ALCOHOL is $\beta_2^\ + \beta_3$ This is pretty straightforward. If you fit a model like this in any statistical package, you would get this kind of result. However, if say, in the model, I actually included an additional variable AGE, which is a covariate (assuming that it does not interact with either SEX or ALCOHOL), the original model would be the one below: **Model 3:** $$ Y = \beta_0^\ + \beta_1SEX\ + \beta_2ALCOHOL\ + \beta_3SEX*ALCOHOL\ + \beta_4AGE\ $$ Further suppose that we still have a significant interaction between SEX and ALCOHOL and I would like to stratify the data again. I would get two models below if I followed the logic above: **Model 4a:** $$ \text{Female: } Y_F = \beta_0^\ + \beta_2ALCOHOL\ + \beta_4AGE\ $$ **Model 4b:** $$ \text{Male: } Y_M = (\beta_0^\ + \beta_1) + (\beta_2\ +\beta_3)ALCOHOL + \beta_4AGE\ $$ However, the actual beta coefficients obtained using a computer program can be very different from what you obtain using the updated equations above. The difference in $\beta_0$ makes sense, because the stratified models in 4a and 4b still assume the pooled mean for age; namely, it estimates $\bar{Y}$ at the mean of all subjects' age, whereas in the stratified analyses done by a computer program, the intercepts of the models estimate $\bar{Y}$ at the mean age of a sub-group. **However, I wonder why the slopes are different. In other words, why are the slopes in Models 4a and 4b different from those produced by a statistical package.**


title: K-means finds less than K groups

text: According to the $K$-means algorithm,   1. Randomly assign a group or cluster to each point. # initial step   2. Compute the centroids of all the $K$ groups.   3. Reassign each point to the closest groups by measuring the distance to the centroid.   4. Iterate over 2~3 until convergence. I think my understanding about this algorithm is right, if not, please correct me. Here is my problem. I have a data set, and I want to find $3$ clusters using $K$-means, and now suppose after the initial step (step-1), it looks like this:               1112            33233      , and the randomly selected cluster-2 is a tragedy, because the centroid of cluster-2 is in the middle of cluster-1 and cluster-3, and after step-2 and step-3, all the points (just two 2's) will be reassigned to cluster-1 and cluster-3, leaving cluster-2 empty. For real? I want 3 cluster, but end up with 2? Just because inappropriate initial state? What do I do with this?


title: Dichotomised data-how to test for a difference(potentially independent t-test) SPSS

text: I am using a 2X3X4 design, which has the following ivs: Social Anxiety Scores(formed two groups,essentially high and low scores) 3 emotions and 4 gaze directions. I need to test if the groups are essentially different(before I do my mixed anova). I thought I could do this through an independent t-test. As I have tried to do one, it has tested the differences between each condition individually. How can I test them for difference, is it a matter of how I present the data in the data viewer? Or is the another method to do it?


title: Procedure to devise Questionnaire

text: I understand that the reliability and validity are key while constructing a questionnaire. Kindly advise if the following procedures are appropriate to devise a questionnaire:   1. After creating all the items in the questionnaire, I plan to test the question with a sample of 25 participants. Next, I will re-test the same questionnaire with the same sample group, at some other time. To check internal consistencies, I will run a regression analysis to find the value of Cronbach alpha or Rho. Please advise.


title: Why does fixed effect p-value in a mixed model act unintuitively?

text: I have a problem with calculating significance of a fixed effect in simple nested design experiment. Suppose we have a simple dataset:               individual  class   ind_mean    value   replicate     1   A   -0.37651119 -1.23495005 1     1   A   -0.37651119 1.24107015  2     1   A   -0.37651119 -0.41977074 3     1   A   -0.37651119 -1.09239410 4     2   A   0.01519211  0.15858308  1     2   A   0.01519211  0.69513257  2     2   A   0.01519211  0.05390856  3     2   A   0.01519211  -0.84685576 4     3   A   -0.19937009 0.49203770  1     3   A   -0.19937009 0.06218100  2     3   A   -0.19937009 -1.12546256 3     3   A   -0.19937009 -0.22623652 4     4   A   0.62853718  1.63792462  1     4   A   0.62853718  1.26212834  2     4   A   0.62853718  -0.52029892 3     4   A   0.62853718  0.13439466  4     5   B   2.38620939  2.67706019  1     5   B   2.38620939  1.71849546  2     5   B   2.38620939  2.31713740  3     5   B   2.38620939  2.83214451  4     6   B   0.91834030  -0.13221319 1     6   B   0.91834030  2.41841628  2     6   B   0.91834030  0.25367274  3     6   B   0.91834030  1.13348538  4     7   B   2.12694664  2.22003887  1     7   B   2.12694664  1.31492118  2     7   B   2.12694664  2.53495794  3     7   B   2.12694664  2.43786856  4     8   B   1.52867659  0.65959322  1     8   B   1.52867659  1.23188018  2     8   B   1.52867659  1.58976692  3     8   B   1.52867659  2.63346604  4      Now if I fit the lme model to this data, I get               > summary(lme(value ~ class, .~1|individual, data = d))     ...     Random effects:      Formula: . ~ 1 | individual             (Intercept)  Residual     StdDev:   0.3643563 0.8430941          Fixed effects: value ~ class                     Value Std.Error DF  t-value p-value     (Intercept) 0.016962 0.2785935 24 0.060884  0.9520     classB      1.723081 0.3939908  6 4.373405  0.0047     ...      Everything seems to be fine. However, I get exactly the same fixed parameter estimates, statistics and p-values if I would fit the model on the column `ind_mean` that instead of the original values has the averages for each individual.               > summary(lme(ind_mean ~ class, .~1|individual, data = d))     ...     Random effects:      Formula: . ~ 1 | individual             (Intercept)     Residual     StdDev:   0.5571871 1.854137e-16          Fixed effects: ind_mean ~ class                     Value Std.Error DF  t-value p-value     (Intercept) 0.016962 0.2785935 24 0.060884  0.9520     classB      1.723081 0.3939908  6 4.373405  0.0047     ...      Actually this p-value coincides with the ordinary t-test p-value, where I would run the comparison on the averages of individuals. This of course makes sense, since the additional observations in this case **do not** add any information. However, in the case of original model the additional observations **do** add information. In fact, the data was generated without any individual influence at all, so the correct statistics would be.               > summary(lm(value ~ class, data = d))     ...     Coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept)  0.01696    0.22597   0.075    0.941         classB       1.72308    0.31957   5.392  7.7e-06 ***     ...      My statistical intuition would suggest that the p-value from the first model should be somewhere between the ones of third and second model, not exactly the same as in the second model. Is my intuition wrong? Am I using the method in a wrong way? Are there some methods that would give me more appropriate p-values?


title: How can a t-test be used to compare the distributions between groups of data?

text: I understand that the t-test is used to test the difference in means for two populations when the populations have relatively similar variances, the units are independent, and they are normal (especially with smaller sample sizes). However, I was wondering how t-tests are used to look at the difference in the distributions of data between two groups? I am asking this because it's basically the phrasing used by the question I'm trying to answer. It asks to compare whether the distributions of the item of interest are different using a t-test. The reason that I am confused is, even though I understand the mean is a product of the distribution and t-tests may be strongly affected by outliers so the t-test might give some information about two distributions, there could be a case where the two distributions were very similar but the effect size was large simply because they were centered at different means, and there could be a case where the two distributions looked funky with difference variances and whatnot, and that could lead to the same t-stat. So how would to be able to tell anything from a t-test?


title: Help with designing the analysis of a radio telemetry study

text: I was wondering if there are any spatial ecologists out there who can provide some insight on the trouble I'm having with my Masters research. I am looking at the association a bat species to high density aggregations of a flowering plant. The bats were tracked throughout the night with a point taken at approximately 10 minute intervals. The bats were tracked over consecutive nights until an asymptote in their space use was reached. So now I have point patterns for each bat that pretty much covers the entirety of their space use at that period. When the data is displayed visually there is a clear pattern of the bats that were tracked during the flowering period moving into the aggregations of the plant however they also spend a lot of time foraging around their roost sites. I have simulated the point pattern of each bat 999 times by randomly shifting and rotating the points within the study window (the internal structure of the observed pattern was maintained in the simulations to account for autocorrelation). I then did randomisation tests to see if the average distance to the nearest plant aggregation for the observed data is within the lower 5% of the distribution of the averages for the simulated point patterns. However because the majority of the points collected are highly associated with the roosts the average distance is primarily effected by the position of the roost in relation to any of the plant aggregations. I've been thinking of ways to disentangle the association with roosts and the association with the plant. First I was thought about putting a buffer around the roosts and removing the points falling within these buffers from the analysis. The problem is deciding on a buffer. In the other areas I have been using a buffer equal to the circular error probable of the radio tracking error but placing this around the roost will only incorporate the points where the bat was actually inside the roost not foraging around it. All the bats have different foraging patterns around their roosts and I can't see how to assign a buffer distance that is not completely arbitrary. The second method we have been looking into is taking out the points that are more closely associated (closer in distance) with the roost than the nearest plant aggregation and then running the randomisation to look at the association of the remaining points to the plant. I want to get an idea of whether or not removing the association with the roosts in this way will potentially bias the analysis. Any thoughts or further suggestions?! Thanks in advance!


title: proof of convergence of EM algorithm with continuous hidden variables

text: I want to find the proof of convergence for EM algorithm in a case when hidden variables are continuous. Do you know the proof of have a reference? Thank you!


title: Statistical significance of infection site

text: I work in a medical environment. I have a situation where I do not know which statistical test would fit well with what I want to find out and thought you might be able to give me guidance. I am investigating the mortality linked to an infection picked up by some patients, and whether the likelihood of site (blood, urine, and sputum) where the bug was found, is significant enough amongst those who died and those who survived. I am not sure whether a chi-squared test is the right one, not that I remember how to use it either! Essentially my population size is 74, and there have been 22 deaths, and most of these were ill with other co-morbidities. The infection is essentially an opportunistic infection. Literature suggests however that those who have it in their bloodstream or part of a pneumonia are more likely to die than those who might have it as a wound infection etc. I am keen to see if my small sample size correlates with this finding.                                 Death   (22)                                 Discharged (52)     Sputum            22%      (5)                                 15%        (15)     Blood             31%      (7)                                 13%         (7)     Urine             22%      (5)                                 31%        (16)        Is there any statistical significance linked to site of infection?


title: 2x2 chi-square test vs. binomial proportion statistic

text: Suppose I'm doing binary classification, and I want to test whether using feature X is significant or not. (For example, I could be building a decision tree, and I want to see whether I should prune feature X or not.) I believe the standard method is to use a chi-square test on the 2x2 table                            X = 0    X = 1     Outcome = 0      A        B     Outcome = 1      C        D      A "simpler" (IMO) test, though, would be to calculate a statistic on the probability that X gives the correct outcome: take p = [(x = 0 and Outcome = 0) + (x = 1 and Outcome = 1)] / [Total number of observations], and calculate the significance that p is far from 0.5 (say, by using a normal approximation or a Wilson score). What are the disadvantages/advantages of this approach, compared to the chi- square method? Is it totally misguided? Are they equivalent?


title: Bayesian inference adaptable to not mutually exlclusive events?

text: Say I have $n$ possible events that lead to $m$ observable effects. The Bayesian inference hypothesis is that events are mutually exclusive and jointly exhaustive. Could I still be able to use Bayesian inference (modified, eventually) when more than one event occurs at some observation? Real example: $n$ possible faulty components on a single board with $m$ on- board test points that can result "in range" or "out of range," and I cannot _a priori_ exclude that two components are faulty at the same time.


title: How to use DLM with Kalman filtering for forecasting

text: Could someone walk me through an example on how to use DLM Kalman filtering in R on a time series. Say I have a these values (quarterly values with yearly seasonality); how would you use DLM to predict the next values? And BTW, do I have enough historical data (what is the minimum)?               89  2009Q1       82  2009Q2       89  2009Q3       131 2009Q4       97  2010Q1       94  2010Q2       101 2010Q3       151 2010Q4       100 2011Q1       ?   2011Q2      I'm looking for a R code cookbook-style how-to step-by-step type of answer. Accuracy of the prediction is not my main goal, I just want to learn the sequence of code that gives me a number for 2011Q2, even if I don't have enough data. Thanks!


title: How to move from some arbitrary "distance" to a probability distribution?

text: I'm doing some object recognition, and when I compare two images, I get some unbounded "distance" between the two images, representing how similar they are. This is somewhat useful, but it seems like it would be more useful if I could move from this "distance" space into a probability space, where instead of `" The distance between image A and image B is X "`, I could get `" P(image is of class B | image A) = X' "`. What is a reasonable way to move from a distance space to a probability space?


title: Discrepancy measures for transition matrices

text: I'm doing some work on modelling transition matrices, and for this I need a measure of discrepancy or lack of fit: that is, if I have a matrix $T$ and a target matrix $T_0$, I want to be able to calculate how far $T$ is from $T_0$. Would anyone be able to provide pointers on what measure I should be using? I've seen some references to using an elementwise squared-error measure, ie sum up the squared differences of the elements of $T$ and $T_0$, but this seems rather ad-hoc.


title: How to choose number of dummy variables when encoding several categorical variables?

text: I'm building a logistic regression, and two of my variables are categorical with three levels each. (Say one variable is male, female, or unknown, and the other is single, married, or unknown.) How many dummy variables am I supposed to make? Do I make 4 in total (2 for each of the categorical variables, e.g., a male variable, a female variable, a single variable, and a married variable) or 5 in total (2 for one of the categorical variables, 3 for the other)? I know most textbooks say that when you're dummy encoding a categorical variable with k levels, you should only make k-1 dummy variables, since otherwise you'll get a collinearity with the constant. But what do you do when you're dummy encoding several categorical variables? By the collinearity argument, it sounds like I'd only make k-1 dummy variables for one of the categorical variables, and for the rest of the categorical variables I'd build all k dummy variables.


title: Loss function of kernlab

text: I'm looking for the default loss function the ksvm()-function from R package 'kernlab' is using. My guess is the Hinge Loss but I cant find any reference or citation for it (There's nothing explained in the documentation). Can anyone tell me the loss function and better show me a link for a citation?


title: Can I use a polytomous categorical independent variable in logistic regression?

text: I found a few similar questions, but none quite had the answer I was looking for. I aim to use logistic regression and have two IVs ("Independent Variables") that each have three categories. I have only seen IVs in ordinal logistic regression that have been continuous or binary categorical previously. I also cannot convert these IVs to a continuous format, as this will change the interpretability of the variables. Am I able to include these IVs in their current state, or should I be looking to use another algorithm here? (E.g. decision tree?) Thanks.


title: Reference/certified data set to compare custom software results

text: I am currently developing a software application that runs and generates results several statistical tests using .Net libraries. I have been comparing my results to R and another (non free) statistical package. For ANOVA test I can compare against NIST data sets, however I was wondering if there were any similar data sets accompanied with "industrial standard" certified p-values I could use for T-Tests and other statistical tests?


title: Can I represent a non-boolean function with a TREE?

text: I have a question in respect to **tree representation** of a **non-Boolean function**. Is there any problem in representing a function from space Boolean-> **Limited** Integer or   **Limited** Integer-> **Limited** integer by a (decision) **tree**? I don’t see any problem, but don’t find any remark about that on the net. Thanks for your contributions


title: Which type of statistical information will need Error Bar charts for presentation?

text: When presenting statistical information using bar charts, when will you need to use Error Bar charts?


title: Mixed model weight matrix

text: Consider Henderson’s classic mixed model equations: ![enter image description here](http://i.stack.imgur.com/u4KLX.gif) What is the nature of the weight matrix, S-1? I believe that this is composed of the covariance matrix, but at the same time it is diagonal. How can that work? If X is 100 X 5, then S-1 is necessarily 100 X 100. However, the covariance matrix is 5 X 5. Does this mean that S-1 has a series of 5 X 5 matrices on the diagonal? If the off-diagonal elements are zero, then are they all 5 X 5 matrices? If this is correct, how does the multiplication work?


title: How to tell stationarity from a sample path?

text: Given a sample path, we can roughly tell whether the mean changes over the time, and, when it doesn't, whether the deviation from mean changes over the time. (Correct me if I am wrong.) But that is not enough for telling if the process is stationary (here, stationarity is defined in terms of the first and the second moments), as we still need to check the sample covariance for every time lag. So can we do that by directly looking at a sample path without estimating the covariance? Thanks!


title: Categorial features in linear machine learning algorithms

text: How can I deal with categorial features inside linear ML algorithms? The obvious solution is to represent each value as a binary feature. For example, if the categorical feature color has three possible values: red, blue and green, then we replace it with three binary features [color == red], [color == green], [color == blue]. Are there better solutions? Specifically, I am using Vowpal Wabbit's logistic regression. As far as I understood from the tutorial VW can operate with categorial features, but I did not find out how does in work.


title: What does weighted cumulative frequency distribution mean?

text: I have two sets of data (temperature and catch) and I am following a proposed method in an article I am reading on the empirical cumulative function (ECDF) analysis. Firstly, I have derived the ECDF for my temperature data using the ecdf function in R. The second step was to get the catch-weighted cumulative distribution for temperature, which I honestly cannot understand the concept. $$W(t)=\frac1n\sum_{i=1}^n\frac{y_i}{y}I_{\\{x_i\\}}$$ where $y_i$ is the catch for day $i$; $y$ mean is the average of catch for all days and I_{{x_i}} is the indication function   How do I compute for this curve? Thank you very much.


title: What is the difference between 'Laplace approximation' and 'Modified harmonic mean'?

text: this question is about Bayesian and computational statistics. I am learning them right now, I have two very common output from my software, one is **Laplace approximation** and the other is **Modified harmonic mean**. Both of them are used for approximation of log marginal likelihood, usually they have very close value. What are their differences? And could anyone provide me some background knowledge no these methods? Thanks!


title: Variable effect in Weibull and Cox models

text: Is it possible to compare effect of variable on the survival time in the Cox model and in the Weibull regression model?


title: Is a vague prior the same as a non-informative prior?

text: This is a question about terminology. Is a "vague prior" the same as a non- informative prior, or is there some difference between the two? My impression is that they are same (from looking up vague and non-informative together), but I can't be certain.


title: What is the practical difference between association rules and decision trees in data mining?

text: Is there a really simple description of the practical differences between these two techniques?   * Both seem to be used for supervised learning (though association rules can also handle unsupervised).   * Both can be used for prediction The closest I've found to a 'good' description is from the Statsoft Textbook. They say Association Rules are used to: > ...detect relationships or associations between specific values of > categorical variables in large data sets. Whilst Decision Tree classifiers are described as being used to: > ...predict membership of cases or objects in the classes of a categorical > dependent variable from their measurements on one or more predictor > variables. However, over at R Data Mining, they give an example of Association Rules being used with a target field. So both can be used to predict group membership, is the key difference that decision trees can handle non-categorical input data whilst association rules can't? Or is there something more fundamental? One site (sqlserverdatamining.com) says that the key difference is: > The decision trees rules are based on information gain while association > rules are based on popularity and/or confidence. So (possibly answering my own question) does that mean that association rules are evaluated purely on how often they appear in the dataset (and how often they are 'true') whilst decision trees are actually trying to minimise variance? If anyone knows of a good description they'd be willing to point me towards then that would be great.


title: Are there R packages to develop spatiotemporal CAR models?

text: I would like to develop a spatio-temporal model with the North Carolina SIDS data set. I would use a hierachical-bayesian model as in http://www.sciencedirect.com/science/article/pii/S1572312711000098 I would like to use the R software: are there packages to perform my analysis? I found _spdep_ for 'pure' spatial CAR & SAR models and _spBayes_ for spatio- temporal bayesian analysis, but this latter seems to work with points and not with regional (areal) data. Any hint is appreciated, both on alternative software to use and on methodologies to study too!