title: Combining contingency tables derived from background and foreground dataset

text: As some of you may now, various statistical tests are used in the task of finding collocations in given text(either a single document or corpus of texts). Chi-square, t-test, likelihood ratio(with binomial distribution assumed). These tests are used to check if two words occur together significantly more often that apart from each other, that's how the strenth of collocation is measured. A-a-anyway, one good thing to do is put your hands onto a large corpus of texts and save contingency tables(or just frequencies, really) of all the bigrams(two words) and unigrams(one word, obviously) from this corpus. We call this information background. Little introduction for you there, now to the main part: when you receive a text(call it foreground), you extract all bigrams from it and have to decide whether this is a collocation or just some pair of words, like "of the" or "and yes", which do not bear any meaning and not at all important in this particular text. Now, to the question: when you stumble upon a bigram that has in fact been seen in background, it means that both unigrams it consists of can be found in background too, all frequencies are in place, everything else is a piece of cake. But what if this bigram has never occurred before? The only place we can find it in is our foreground, but there's not enough information to infer the distribution and collocation strength. I mean, you can apply any of these statistical tests to foreground contingency table, but it won't be the same, you know? Not that reliable, I suppose. So, at last, to the question itself, all two parts of it: 1) What to do with bigrams that have never occurred in background corpus(and both unigrams are not there either, just two words we haven't seen before) 2) What to do with bigrams only one part of which has occurred in background corpus? 3) What to do with bigrams that have never occured in background corpus, but both words from it are present in background corpus? How to make it smooth, so some scores of collocation strength can be created for all these cases and were at least a little reliable. Thank you in advance, kick this question to whatever site is appropriate if I've misunderstood the purpose of this particular part of stackexchange. _**UPD_** : my 9th attempt to post it ended up with trimmed question body, I'm sorry, now it's intact.


title: Unsupervised anomaly detection with factor analysis (in R)

text: The basic idea i'm trying is to model the data with factor analysis, **assuming a latent variable structure** that underlies the observations. Labels for "real" anomalies are available and used for validation. Another important note is that the data does not have a "very" Gaussian nature. Then, projecting the observations onto the new subspace found, and **looking for outliers** \- i.e. finding the observations that digress from the underlying structure, assuming the anomalies i'm looking for will be more prevalent there. Using R's _factanal_ , with _varimax_ (or _oblimin_ ) rotations, i have encountered the following result, which i'm not sure about :   * When projecting the **normalized** data (z-score) using the FA model found, classification/detection is poor ; but since factanal scales the data in advance, i believe this is the correct way(?)   * When projecting the **original** data (not normalized) using the same FA model, detection is much improved.  Could this be due to projecting the observations onto a subspace that assumes a low variance, much more "normal" distribution, thereby "aggravating" any anomalous points? Would appreciate help/insight as to what i'm missing! Thank you. **Update:** \- LX yields much better classification than LX*, where L is the loadings matrix, X is the data matrix, and X* is the standardized data matrix. LX also yields better classification than in the following. https://stat.ethz.ch/pipermail/r-help/2002-April/020278.html Could the LX transformation be related to LDA?


title: Estimating specific variance for items in factor analysis - how to achieve the theoretical maximum?

text: _(Remark: this is a "scholastic" question - I'm reviewing my implementation of factor analysis procedures; I'm not looking for good approximations for an actual survey/actual data or the like.)_ There are different methods of estimating individual variances for the items in a covariance matrix $C$ with $m$ rows and columns; I know that method of using $D$, the reciprocal of the diagonal of the inverse of $C$. Well, it shall result in a Heywood-case/in negative definiteness of the remaining matrix if I simply remove that variance from the diagonal of the covariance matrix ($C - D$ is surely negative definite); but I can iteratively determine the greatest possible part in $D \cdot 1/r$ to be removed which still keeps the covariance $C - D/r$ positive semidefinite. This gives then a certain sum of that itemspecific variances ($s_1=sum(D/r)$).   Another method is to get the least principal axis $A_m$ , take the eigenvalue $\lambda _m $ , then norm the other axes to the same length $B_k = A_k \cdot \lambda_m / \lambda_k $ and in the diagonal of $ B \cdot B^T$ we get (equal) itemspecific variances. Note that again $ C - B \cdot B^T $ has reduced rank, and thus "all individual variance" is removed - however, the sum of all that itemspecific variances $s_2$ is usually much smaller than than $s_1$.   After that two different solutions, already leading to different amounts of overall itemspecific variance removed, I experimented with further different methods and one gives $s_3$ which is even greater than $s_1$. And having now a handful of further methods with different values $s_j$, the question naturally occurs: Q: is there a specific method, which allows to extract the maximally possible _sum of individual variances_ of a covariance matrix, and if there is a special method, how is it defined? * * * To see, that the differences between the methods are not simply neglectable I add an example with some test-covariance matrix.   Overview, comparision of 4 methods: ![enter image description here](http://i.stack.imgur.com/MmgKu.png) Detail 1:   ![enter image description here](http://i.stack.imgur.com/kldJ9.png) Detail 2: I'm surprised that the shape of the approaching of the maximum has such a spike - I'd expect some smooth "top of a normal-curve" here:   ![enter image description here](http://i.stack.imgur.com/zKLp8.png)


title: Different ANCOVA results in R using different functions

text: before my actual ANCOVA I would like to test whether there is a significant interaction between the IV and the CV as this is one assumption for ANCOVA. I found 3 different ways in R to perform an ANCOVA. However, the result for one solution differs from the other two and I do not understand why. Here is my working code snippet:               library(lattice)          data <- data.frame(group = c(rep("CTRL", 10), rep("P", 10)),                              response = c(10,11,14,16,17,17,19,20,21,22, 10,11,11,11,12,13,14,14,15,16),                             age = c(40,41,45,43,50,51,55,57,60,62, 30,32,34,35,40,41,42,44,43,46))          xyplot(response ~ age, data=data, groups=group, type=c("p","r"))          # 1. ANCOVA         anova(lm(response ~ group + age + group : age, data = data))          # 2. ANCOVA     summary(aov(response ~ group + age + group : age, data = data))          # 3. ANCOVA     summary(lm(response ~ group + age + group : age, data = data))      I don't understand why the two p-values for group and age are identical for ANCOVA 1. and 2. but different for 3. even though the interaction p-value (group : age) is the same for all three. Doing the same thing in SPSS results in exactly the same p-values as ANCOVA 3. Now, as you can imagine, I am pretty unsure what's right and what's wrong or what's actually the difference between them? Can anyone help?


title: When am I allowed to make linear regression on a small sample?

text: I'm trying to evaluate the linear correlation between to continuous variables (% value calculated from EEG datas and the area of an anatomical region of the brain). I have a sample of 18 right now. From what I understand, I need to make sure that the x and y variables vary together in a joint distribution that is normally distributed. Or at least that the distribution of each of these variables is normal. My question is what do I have to do to make sure that the assumptions needed to draw valid conclusions about the correlation coefficient are met? I though about evaluating the normality of the distribution with the Shapiro–Wilk test for each variable and then simply calculate the Pearson R. Would that be valid? Would it be better to draw a normal probability plot? Thanks!


title: Random forests visual introduction-level reference or tutorial

text: I have seen several references, but am looking for something easy to follow that illustrates Random Forests in regression and feature importance applications. I want to make sure that I explain this at the 101 level. Something like this, but more pictorial like this pdf. I am having a hard time getting across a concrete high-level example of how predictors are being chosen and making them tangable to a general audience.


title: Absolute values -- parametric or non-parametric test?

text: I am looking at the absolute difference in values between pairs of numbers. I am then comparing these absolute values between two groups (e.g. Males and Females). Since these are absolute values I am not sure whether it is appropriate to use parametric statistics (e.g. Independent t-test) or whether it is better to be using non-parametric tests (i.e. Mann-Whitney U)? Specifically I am worried that being absolute values (created by ignoring the signs) they cannot fit the assumption of normality needed? The largest sample size I'm looking at is 19 in Group 1 and 17 in Group 2 (however there are also some smaller groups). Any help would be appreciated (and explanations why would be great, as I'm trying to wrap my head around all this). Thanks!


title: PCA or method of principal component for FA?

text: I have just read one of the answer from a member as follow: "One of the biggest reasons for the confusion between the two has to do with the fact that one of the factor extraction methods in Factor Analysis is called "method of principal components". However, it's one thing to use PCA and another thing to use the method of principal components in FA. The names may be similar, but there are significant differences. The former is an independent analytical method while the latter is merely a tool for factor extraction" see this answer here is there any body explain a little more about this, especially about PCA and method of principal component in FA?


title: Factor vs. Covariate in SPSS MIXED

text: I am using the MIXED command in SPSS to run a multilevel model and was wondering what the difference is when using a dummy coded factor either as factor (via the BY command) or as a covariate (via the WITH) command. A similar question was already posted here: BY and WITH command in mixed model SPSS In this thread the issue, however, was the NOINT option. In my case I am using the intercept and running a simple actor-partner model with 2 metric covariates and one dummy variable (0/1) that distinguishes between two groups. What I am doing is to interact the dummy with both metric variables. Now, if I use the dummy as a factor (via BY) I get different results for the two main effects, as compared to when I use the dummy as a covariate (via WITH). Aso the estimate and t-value for the Intercept differ. Can anybody shed light onto this issue and explain why this is the case? I am puzzled.


title: Cross validation for lasso logistic regression

text: I am writing a routine for logistic regression with lasso in `matlab`. So the problem is to minimize the negative log-likelihood function with the penalty term $$\sum \left(\log(1 + e^{X_i' \beta}) - y_i X_i' \beta\right) + \lambda \sum |\beta_i|$$ where $\beta$ is the model parameter, $X_i$ is the $i$th row of matrix $X$, and $y_i$ is the value of observation $i$. My first question is for a 5-fold cross-validation, which criterion should I use to pick the best value of $\lambda$? Should I use the value of the logit function on the validating data set or mis-classification rate on the validating data?


title: How to associate daily rainfall probability with historic data

text: Given a rainfall forecast of `d days` ahead and historic data collected over `y years`, what's a simple (but correct) way to associate the rainfall probability for each forecast day? I know this is a complex problem, but my purpose is to quantify the probability of a given rainfall volume 'happen' (ie, observed rainfall >= forecast for that day). My first take was to use a simple frequency approach over the years:               target_day = 15     target_month = 1      target_day_forecast = 5          data = select over years such that day = target_day AND month = target_month     n = number of days where forecast >= target_day_forecast     total = number of points in data          return n/total      On the other hand, i had trouble fitting a `gamma` distribution:                require(MASS)      fit = fitdistr(serie, "gamma")           #Error in optim(x = c(0, 0, 0, 0, 0, 0, 0.3, 0, 0.1, 0, 0.8, 0, 0, 0, 0,  :       #initial value in 'vmmin' is not finite      If a sum a small quantity to `serie`, it works:               fit = fitdistr(serie + 0.1, "gamma")      My idea is to fit a `gamma` for each month, then i can extract probabilities easily. So, how can i handle the `gamma` fit problem above? I also appreciate any ideas on this, Thanks!


title: Probability of a random outcome given a sample size

text: Say you gave an IQ test to a random sample of $n$ people, but instead of having the actual test results, you have only the normalized data for this sample (mean=0, $\sigma$=1). What is the probability that person $A$ is in the top 5% of the total population (1.64 standard deviations above mean), given that their normalized IQ score $S$ is $X$ standard deviations above the sample mean? I'm approaching that question from a Bayesian perspective (which is new to me), so I define it as: $$ P(A>1.64 | S=X)= \frac{P(S=X|A>1.64) \times 0.05}{P(S=X|A>1.64) \times 0.05 + P(S=X|A<1.64) \times 0.95} $$ The probabilities I'm after are: $P(S=X|A>1.64)$ and $P(S=X|A<1.64)$. I'm trying to use frequentist statistics to determine probability of a sampling error in the sample population. For example: For $X>1.64$, $P(S=X|A<1.64) = $ probability that the sample mean $x$ would be inflated by $X-1.64$, or the sample standard deviation would be dampened by $X-1.64$. Is this the right approach? And if so, how would you determine the same probability for $X<1.64$? Here are my calculations for the probability of the mean and standard deviation being wrong: mean: $$ \mu = X - 1.64, x = 0, \sigma=1 \\\ z = \frac{x - \mu}{\frac{\sigma}{\sqrt{n}}} \\\ z = \frac{1.64 - X}{n^\frac{1}{2}} $$ standard deviation: $$ df = n - 1, \sigma = 1.64 - X, s = 1 \\\ z = \frac{(n-1) \times s^2}{\sigma^2} \\\ z = \frac{n - 1}{(1.64 - X)^2} $$


title: R: Weights in glm-function

text: I'm analyzing a vote result in my country and I'm using a logistic model for this task (glm). The population was asked to vote yes or no, which resulted in an approval-value from 0% to 100% for each municipality. I would like to weight this result based on the people who have voted, since in small municipalities only 100 people voted and in other more than 50'000. Is it sensible to put the amount of votes as a weight in the glm function like this:               vote.glm <- glm(vote_result ~ explaning_variable_1 + explaning_variable_2, weights=amount_of_votes, family=binomial(logit), data=data)      I'm asking because the summary and anova-function both show that all variables are significant, which makes me suspicious. Thanks!


title: Time series prediction: visualising path uncertainty region

text: I am predicting a time series' future evolution and am evaluating the path uncertainty using bootstrapping. Is there a good way to visualise the uncertainty that goes beyond simply plotting a pair of confidence bands, e.g. 90% lines? I thought of something that gives a better feeling for the distribution at each future time point: e.g. something like a contour plot or a plot where the level of transparency indicates the likelihood. Also, is there a convenient way to compare two path predictions from different models visually? The only thing that comes to mind is two overlaid contour plots.


title: How to ensure properties of covariance matrix when fitting multivariate normal model using maximum likelihood?

text: Suppose I have the following model $$y_i=f(x_i,\theta)+\varepsilon_i$$ where $y_i\in \mathbb{R}^K$ , $x_i$ is a vector of explanatory variables, $\theta$ is the parameters of non-linear function $f$ and $\varepsilon_i\sim N(0,\Sigma)$, where $\Sigma$ naturally is $K\times K$ matrix. The goal is the usual to estimate $\theta$ and $\Sigma$. The obvious choice is maximum likelihood method. Log-likelihood for this model (assuming we have a sample $(y_i,x_i),i=1,...,n$) looks like $$l(\theta,\Sigma)=-\frac{n}{2}\log(2\pi)-\frac{n}{2} \log\det\Sigma-\sum_{i=1}^n(y_i-f(x_i,\theta))'\Sigma^{-1}(y-f(x_i,\theta)))$$ Now this seems simple, the log-likelihood is specified, put in data, and use some algorithm for non-linear optimisation. The problem is how to ensure that $\Sigma$ is positive definite. Using for example `optim` in R (or any other non-linear optimisation algorithm) will not guarantee me that $\Sigma$ is positive definite. So the question is how to ensure that $\Sigma$ stays positive definite? I see two possible solutions:   1. Reparametrise $\Sigma$ as $RR'$ where $R$ is upper-triangular or symmetric matrix. Then $\Sigma$ will always be positive-definite and $R$ can be unconstrained.   2. Use profile likelihood. Derive the formulas for $\hat\theta(\Sigma)$ and $\hat{\Sigma}(\theta)$. Start with some $\theta_0$ and iterate $\hat{\Sigma}_j=\hat\Sigma(\hat\theta_{j-1})$, $\hat{\theta}_j=\hat\theta(\hat\Sigma_{j-1})$ until convergence. Is there some other way and what about these 2 approaches, will they work, are they standard? This seems pretty standard problem, but quick search did not give me any pointers. I know that Bayesian estimation would be also possible, but for the moment I would not want to engage in it.


title: Confidence interval of a half-circular parameter

text: I am modelling two dimensional ovals on a two dimensional plane. Each oval is defined by five parameters, namely `x, y, s1, s2, t` where `x` and `y` are its position on the plane, `s1` and `s2` are length of its long and short axis respectively and `t` is the angle that the long axis of the oval makes with the y-axis i.e. `t` is the direction of the oval. As the ovals are symmetric `t = t+pi` so, `t` is set to a range of -0.5*pi to 0.5*pi. To acquire the mean and 95% confidence interval of most of these parameters bootstrapping is used and the 2.5 50 and 97.5 percentile are taken. However, as `t` is a circular parameter it may occur that some of the outcomes of the bootstrap are slightly above -0.5*pi and others are slightly below 0.5*pi. This prevents the use of a percentiles for the confidence interval. While a mean can easily be acquired by converting the angles to a cartesian plane and doing the math from there, I am clueless as to how to find a measure for 95% confidence or a standard error.


title: General method for deriving the standard error

text: I can't seem to find a general method for deriving standard errors anywhere. I've looked on google, this website and even in text books but all I can find is the formula for standard errors for the mean, variance, proportion, risk ratio, etc... and not how these formulas were arrived at. If any body could explain it in simple terms or even link me to a good resource which explains it I'd be grateful.


title: Doing comparisons from two null distributions

text: I have a situation where I have an observed statistic computed from data, and then I have approximated the null distribution by some sort of resampling. I have used this to calculate p-values for a one-sided hypothesis test. Now, if I want to do a comparison between two statistics (are they significantly different?), but I already have the samples from the two null distributions, is there some way to use these to approximate the null distribution for the difference?


title: Probability of object collision

text: I've been searching for a few days on a number of sites but I can't seem to find a good answer for this. I'm developing a collision detection program using Unscented Kalman Filter and predicting possible positions of different objects. These objects have a certain area but my predictions are of course only propagating the objects center of mass and I get the respective variance for their positions. That is, I get the position of their **center of mass** is given by the bivariate normal distributions $N_1$~$(\mu_1,\Sigma_1)$ for object 1 and $N_2$~$(\mu_1,\Sigma_2)$ for object 2 in Cartesian coordinates. Now the objective is finding out the probability of them colliding with each other. If both the objects have area $A$, what is the probability of them colliding ? I know how to assess this threat when only using the objects center of mass, but that distribution doesn't really contain the objects area. For example in cases when $\Sigma_1$ extremly small, I "know" that the objects center of mass is there and from that the object's edges (using the area, or if it's better, the objects $Width$ and $Length$) So, long story short, is there a neat way to transform $N_1$ and $N_2$ to incorporate the area ?


title: Reconstruct a vector with a known vector and residual

text: I observe $\vec y$ and know $\vec x$. I assume that $\vec y$ mostly consists of $\vec x$, with some added residual $\vec r$. This gives me the problem $\vec y = a\vec x + \vec r$, where $a \in [0, 1]$. Furthermore, I assume that we should maximize $a$. Otherwise, $a = 0, \vec r = \vec y$ would be a trivial solution. Another way to express the problem for vectors of length 3 would be: $\begin{pmatrix} y_1 \\\ y_2 \\\ y_3 \\\ 1 \end{pmatrix} = \begin{bmatrix} a & 0 & 0 & r_1 \\\ 0 & a & 0 & r_2 \\\ 0 & 0 & a & r_3 \\\ 0 & 0 & 0 & 1 \\\ \end{bmatrix} \begin{pmatrix} x_1 \\\ x_2 \\\ x_3 \\\ 1 \end{pmatrix}$ In this version, the matrix would need to be calculated, but we need to do it in a way that maximizes $a$. How to best approach this problem?


title: Equivalent definitions of a sufficient statistic

text: I'm trying to understand the definition of a sufficient statistic for continuous random variables given in Introduction to Mathematical Statistics by Hogg and Craig (7th edition). Let $X_1,X_2,...,X_n$ be a random sample with joint pdf $f(x_1,x_2,...,x_n;\theta)$, $\theta \in \Omega$, and $T(X_1,...,X_n)$ be a statistic with pdf $f_T(y;\theta)$. Wikipedia, among other sources, defines $T$ to be sufficient for $\theta$ if and only if the conditional distribution of $X_1,X_2,...,X_n$, given $T=t$, does not depend on $\theta$. On the other hand, Hogg and Craig defines $T$ to be sufficient for $\theta$ if and only if $$\frac{f(x_1,x_2,...,x_n;\theta)}{f_T(T(x_1,...,x_n);\theta)}$$ does not depend on $\theta$. Here's my question: are the two definitions equivalent, and if so, how does one prove this?


title: What happens if a survival curve doesn't reach 0.5?

text: Does this mean you can't compute the median?


title: Site effect in GAM model

text: I am trying to build GAM model to see the effect of several environmental variables on the total abundance of one species. I have collected samples from three sites with three replicates from each site. How is the best for me to see the effect of each variable on the model after removing the variation cause by the three sites and considering the three replicates as random effect. My try is the following :               m01<- gam(TotalInd ~ s(dayinyear, bs="cr") + Site2 + s(Temp) + s(Mud) + s(Chlorophyll) + s(Salinity) ,random=list(Replicate2=~ 1), data=Data)      Do you think the format of this model is correct?


title: Statistical Test for seeing differences in a single sample

text: So I've been trying to analyze 3 repeats for some data (essentially, a single independent variable measured at three different times), which would look like something like this: Test 1 [Value] Test 2 [Value] Test 3 [Value] Now, I want to see if these values are different from one another, but I don't exactly know how to do that when each test only had a single data point taken (no ANOVA). Is there a statistical test to determine if these values are different, without having a known mean?


title: multiple predictors in decision tree model

text: I am using like 10 predictors in my decision tree, but the rpart function uses only like 8 of them. does it mean that rest 2 are not needed or they become redundant? like i have age and child coded based on age as 1/0. but the tree uses only age and not child.


title: Generating random matrices with specific equality constraints

text: Suppose I want to generate a nonnegative $n \times n$ matrix $\mathbf A$ for an odd $n$ (say, $n=5$ for a good enough example), such that   * the individual elements are drawn from a uniform distribution but with the equality constraints   * the sum of all the elements in the rows above the mid row = the sum of all the elements in the rows below the mid row   * the sum of all the elements in the columns left of the middle column = the sum of all the elements in the columns right of the middle column   * the sum of all elements = 1.0 That is, in effect I want to generate matrices "balanced" in a specific way (I intend to use them as convolution matrices and I don't care about translation). For my current purposes I think I have a good enough ad hoc (not entirely correctly distributed) solution that generates a uniform matrix and then enforces the equality constraints by correcting the individual elements using a least squares approach. However I'm curious, mainly for the learning experience, if there would be some easy-ish way to sample the resulting conditional distribution. So far I've had a quick peek on Wikipedia into the Bates distribution (mean of uniform random variables) and into Monte Carlo Markov Chain models, but I'm not sure if they help here or if there would be an easier solution. I'm not necessarily looking for a complete solution; pointers to insightful or helpful material are also appreciated. **Edit:** I'm looking for the conditional probability distribution, i.e. the elements do not need to be uniformly distributed; rather, I'm looking for a distribution where each sample that satisfies the equality constraints has the same probability (hence conditional uniform distribution given the equalities).


title: How do you plot survival functions for specific values of covariates in SAS?

text: Suppose I do the following:                   proc phreg data = new;               model time*censor(0) = x y;         run;      Also suppose $x$ is a binary variable and $y$ is a continuous variable. How would I plot the survival function for $x = 1$ and $y = 100$ (for example) in SAS?


title: Why do these 2 approaches to applying mixed models yield different results?

text: I'm re-analyzing a colleague's data. The data and R code are here. It's a 2x2x2x2x3 completely within-Ss design. One of the predictor variables, `cue`, is a two-level variable that when collapsed to a difference score reflects a theory-pertinent value. She previously collapsed `cue` to a difference score within each subject and condition, then computed an ANOVA, yielding an MSE that she could then use for planned comparisons of each condition's mean difference score against zero. You'll have to trust me that she wasn't fishing and did indeed have good theoretical basis for doing all 24 tests. I thought I'd see if there was any difference when instead using mixed effects models to represent the data. As shown in the code, I took two approaches: Method 1 - Model the data as a 2x2x2x2x3 design, obtain _a posteriori_ samples from this model, compute the `cue` difference score for each condition within each sample, compute the 95% prediction interval for the cue difference score within each condition. Method 2 - Collapse `cue` to a difference score within each subject and condition, model the data as a 2x2x2x3 design, obtain _a posteriori_ samples from this model, compute the 95% prediction interval for the cue difference score within each condition. It appears that method 1 yields broader prediction intervals than method 2 with the consequence that if one uses overlap with zero as a criterion for "significance", only 25% of the cuing scores are "significant" under method 1 while 75% of cuing scores are "significant" under method 2. Noteably, the patterns of significance obtained by method 2 are more akin to the original ANOVA-based results than are the patterns obtained by method 1. Any idea what's going on here?


title: Adjusting existing algorithm - likelihood for presence-only data

text: Logistic regression fits a model that predicts a binary variable whilst performing a logit transformation of the linear combination (LC) of predictors: 1/1 + exp(-LC). I have a working machine learning algorithm which does the same but fits a non linear combination of predictors by maximising the log likelihood. I would like to apply this algorithm to presence only data as suggested here: Presence-only data and the em algorithm Could I just use formula (4) (its log likelihood version)? How do I use the resulting fitted model to perform predictions if I do not know the value of z but just the predictors/co-variates? I also do not understand the difference between observed and full likelihood? Any clarification would be very much appreciated.


title: Include a distance matrix as a random effect in a mixed-effects model

text: I am trying to construct a mixed effects model and to include in it a measure of the (genetic) distance between observations, so that closer observation from different main effect levels will **outweigh** distant observations from different main effect levels.   I have constructed a similarity matrix for this purpose but the question now is how to incorporate it (or a transformation of it) in the model. From my understanding of the `lmer` function in `R`, random effects can only be a `data frame` column not a matrix. Any ideas how to get that incorporated?


title: Finding optimum point of parameters

text: I have an algorithm with 3 parameters and sum of these parameters is equal to one; $a_1+a_2+a_3=1$ and each of them must be between $0$ and $1$. I want to find the optimum point for this parameters. The equation is a linear combination of these three parameters, $y=a_1 g_1+a_2 g_2+a_3 g_3$ and $g_1$, $g_2$ and $g_3$ are some known constants. If I use k fold cross validation to find the optimum point the space of search is very large because the dataset is large. Is there any test method for such cases? If I want to search the optimum point for example by k-fold and by taking two of them constant and change the other, the space become large. Any suggestion? Is there better method for this? Thanks.


title: Why does increasing the number of bootstrapped cases make PLS coefficients significant?

text: I am running a PLS model with a low number of observations ($n=50$). While several pieces of academic work argue that this sample size is appropriate to run this type of model, I am quite confused when it comes to the number of cases I take in my bootstrapping procedure. I get some very strong path estimates for certain paths without a significant t-statistic (50 cases, 500 samples). However, when I put the number of cases above n (lets say 100 cases, 500 samples), the t-values are more in line with the path estimates.


title: Implications and interpretation of conditional independence assumption

text: I am currently studying Treatment Effect Analysis and I am reading about the **Conditional Independence Assumption (CIA)** : \begin{align} (Y_1,Y_0)\perp D|X \end{align} So the outcomes are independent of the treatment, conditional on X. First Question: Can someone explain this better? So give a real interpretation of this and not just the mathematical formulation verbally? The CIA implies **mean independence**. Second Question: What does this implication mean? So what is the connection about the CIA and the measurement of the treatment effect? I am not used to the topic and tried to read typical literature, but they are not really giving an intutively interpretation.


title: Is it ok to adjust significance level according to number of comparisons in simple correlations?

text: a referee asked me to "split" the significance level for the number of comparison, although i am not doing repeated comparisons but just spearman correlations between variables. I have set the significance level at p<0.05 but the referee is asking me to reduce it since the number of "correlations" is high. I wonder if it is ok to do so, or if I should reply to him stating that this is wrong (in which case i should need some reference if you have some). When is it necessary to adjust significance level for the number of comparisons? Thank you in advance


title: Reporting main effects and interaction from ANOVA

text: When reporting a two way 2 x 3 mixed ANOVA I have recently heard that the interaction effect should only be reported if the two main effects are significant (APA style). Is this correct?


title: Bivariate Gamma distribution PDF

text: I'm analyzing a set of data, and I like to fit a gamma distribution. I know how to do it in one dimension, but the data that I'm analyzing now are two dimensional. Is there any way that I can have a bivariate gamma function fitted into the data?


title: Why is the conditional mean of the reduced form error zero?

text: For example, we have a simultaneous equation model of supply and demand: Supply: $$s(p)=\alpha_{s}+\beta_{s}p+\epsilon_{s}$$ Demand: $$d(p)=\alpha_{d}-\beta_{d}p+\epsilon_{d}$$ Market clearing condition: $$q = d(p) = s(p)$$ It is obvious that $E(\epsilon_{s}|p) \neq 0$ and $E(\epsilon_{s}|p) \neq 0$, as price $p$ is endogenous. The reduced form of this model is: $$p = \frac{\alpha_{d}-\alpha_{s}}{\beta_{d}+\beta_{s}}+\frac{\epsilon_{d}-\epsilon_{s}}{\beta_{d}+\beta_{s}}$$ $$q = \frac{\alpha_{d}\beta_{s}+\alpha_{s}\beta_{d}}{\beta_{d}+\beta_{s}}+\frac{\beta_{d}\epsilon_{s}+\beta_{s}\epsilon_{d}}{\beta_{d}+\beta_{s}}$$ It follows that the reduced form errors are linear combinations of the structural errors. My question is, how can one derive that the conditional means of the reduced form errors, given the constants, are zero?


title: Using heckman in combination with mi estimate (Stata)

text: This is a Stata specific question, and may be better directed to Stata's own Statalist, but I'm trying here first. Here is the situation. I have a multiply-imputed data set, where the number of imputations = 10 and where imputations are only performed on independent variables. I would like to run the `heckman` command on these data, using the `mi estimate` prefix. The `heckman` command is not a "supported command", but it is possible to force estimation using the `cmdok` option. My questions are these: Why isn't `heckman` a supported command? Are the results that I get when I force Stata to fit the model (using the `cmdok` option) valid? Or, is there some statistical reason why a Heckman-style selection procedure should not be used in a multiply-imputed data setting? Any guidance on this matter would be much appreciated.


title: Variable cost - benefit analysis on random forest

text: I have a random forest being trained with _n_ vectors each with _m_ variables. Each variable has a cost based on how much time it takes to compute it ( _m1_ might take 1 unit while _m2_ might take 100, making it more "expensive"). As far as I understand, random forests can give you the contribution each variable make to the whole forest. My variables are not independent, so there is some correlation between some of them. **What is the correct way to determine the variable importance of the variables making consideration of their cost?**


title: Why use the F distribution and F test?

text: I don't understand why in the F test we calculate the ratio between MSE between subject and MSE within subject. As far as I know, this is due because we want to use the F distribution, which is a rate between two $\chi^2$ distribution divided their degree of freedom. My question is: Why is not used a simple version of this method? Couldn't we just take the ratio of the sum of squares between subjects with the sum of squares within subjects? Then we could compare this outcome with a distribution that is the result of the ratio between two $\chi^2$ distribution, instead of dividing the numerator and denominator for the degrees of freedom. Does it makes sense? Why they use the other way?


title: Will Z-scores help me compare my variables?

text: I have two variables: 1) Percentage of population of different geographic regions which has bought a certain category of product (say hair sprays) in the last MONTH 2) Percentage of population which has bought a specific brand of hair spray in the last YEAR. I am trying to in which region the brand is doing better compared to the general use of hair spray, and in which it is doing worse. Since the first variable is for a month and the second for a whole year, i can't just compare them like that. So I thought i'll transform them into z-scores, and take the difference between the two. Question: Is that a valid approach? (Side question: the distribution of the second variable has quite a strong positive skew, while the first is fairly normal. Is that a problem?) Thanks in advance for any advice! EDITS in response to questions:   * For variable 2 I only have data for that one particular brand, so I can't aggregate across brands. While variable 1 is only for the whole product category encompassing all brands - can't dissaggregate that either.   * Sample size is over 100 geographical regions   * I am attaching a scatter plot of the two variables. One of which shows variable 2 on a log scale.    * Just to clarify, I am not so much interested in the overall relationship between the two variables, but more in arriving at a metric that will tell me for any specific entity (region) how well it is doing on variable 2 (sales of specific brand) compared to variable 1 (sales of general product category). Now, of course, I can see that maybe the regression/scatterplot approaches suggested might be an intermediate step to achieve that goal. But how do I progress after that? Thanks so much for the suggestions so far though! ![scatterplots](http://i.stack.imgur.com/cdR3k.jpg)


title: How can I calculate prediction error rate in logistic regression?

text: Can anybody lead me with a small simulation? Should I calculate prediction error rate from a classification table while splitting the data into two parts then fitting on the training data and predicting on the test data and see the misclassification (i.e. how many 1 predicted as 0 and 0 as 1, calculating the rate of misclassification). I tried the following R code:               # prediction error rate for logit and probit for different sample sizes and sigma:          mu <- rep(0,4)        library(MASS)       pred.err <- function(sam.size,sigma,cut.point.latent){       d <- mvrnorm(sam.size,mu,sigma)       obs <- ifelse(d[ ,1]>cut.point.latent,1,0) # cut point for y*(latent) to be y(observed)       x1 <- d[ ,2]       x2 <- d[ ,3]       x3 <- d[ ,4]       dat <- data.frame(obs,x1,x2,x3)       sam.select <- sample(1:nrow(dat),size=nrow(dat)*0.6) # 60% training data       train.dat <- dat[sam.select,]           test.dat <- dat[-sam.select,]       fit.logit.train <- with(train.dat,glm(obs ~ x1+x2+x3,family=binomial(link="logit")))       fit.probit.train <- with(train.dat,glm(obs ~ x1+x2+x3,family=binomial(link="probit")))       test.covar <- test.dat[ ,2:4]       pred.prob.logit <-predict(fit.logit.train,newdata=test.covar,type="response")             pred.prob.probit <-predict(fit.probit.train,newdata=test.covar,type="response")       obs.res <- test.dat[ , 1]       cut.point.prob <- pnorm(cut.point.latent) # as cut point should be same for logit and probit       pred.res.logit <- ifelse(pred.prob.logit >= cut.point.prob,1,0)       pred.res.probit <- ifelse(pred.prob.probit >= cut.point.prob,1,0)       logit.obs.pred <- data.frame(obs.res,pred.res.logit) # as in pred.res the columns        #are different from obs.res       probit.obs.pred <- data.frame(obs.res, pred.res.probit)       tab.logit <- table(logit.obs.pred)          tab.probit <- table(probit.obs.pred)        logit.pred.err.rate <- sum(tab.logit[1,2]+tab[2,1])/sum(tab.logit)       probit.pred.err.rate <- sum(tab.probit[1,2]+tab[2,1])/sum(tab.probit)       list(logit.pred.err.rate,probit.pred.err.rate)      }          pred.err.logit <-NULL     pred.err.probit <-NULL     result <- function(rep.no,sam.size,sigma,cut.point.latent){       for(i in 1:rep.no){         sam <-pred.err(sam.size,sigma,cut.point.latent)         pred.err.logit[i] <- sam[[1]]                 pred.err.probit[i] <- sam[[2]]       }             list(pred.err.logit,pred.err.probit)           }          sigma.high <-matrix(c(1,.41,.5,-.7,.41,1,0,0,.5,0,1,0,-.7,0,0,1),nrow=4)          output <- result(rep.no=1000,sam.size=500,sigma=sigma.high,cut.point.latent=.53)          ### findings :          t.test(output[[1]],output[[2]])


title: Selecting problems of the appropriate difficulty based for adaptive learning

text: I'm currently working on an adaptive learning system for high school maths. Students complete questions in quizzes and I need to be able to select questions of the appropriate difficulty level (say approximately 75% chance of being correct) and create a competency metric for measuring the skill of a student at a particular subtopic. We don't have explicit information about the difficulty of the questions or the abilities of the students - both need to be inferred from the data. The questions are broken down into topics and sub- topics (each question is assigned a single topic and sub-topic) and we have a dependency graph indicating which sub-topics depend on other sub-topics, so we can take into account the student's performance on areas outside of the subtopic. How can I estimate the chance of a student getting a question correct based on historical data and create metrics for the difficulties of questions and skill levels of students?


title: Obtaining sample variance from grouped data for goodness of fit test

text: This is a practice question I came across when dong some goodness of fit test examples.A company sells cloths by mail order.The size of clothes is defined by hip size; thus the height of customers of a particular customer may vary considerably.   Data set of heights sent in by customers of size 18 is given.Data set includes class interval,class mid mark and frequency. a) Given that $\sum$f=265, $\sum fx$=40735, $\sum fx^2$=6299425 estimate mean and standard deviation of heights. b)The company decides that it is not possible to produce a range of garments for a particular size suitable for customers of different heights.A single height must be chosen and it is proposed that mean height should be chosen for this.Comment on this suggestion as it applies to customers of size 18 and make an alternative proposal. My question is how can I calculate **sample** variance from $\sum$f=265, $\sum fx$=40735, $\sum fx^2$=6299425   From V(X)=E(X$^2$)- E[X]$^2$ what I get is population variance and not sample variance right? For sample variance I should divide by n-1.So how can I get **sample** variance value? Also why does it ask to **estimate** values for standard deviation?This question hasn't specified that a sample is being taken.So won't these be population variance and mean I am calculating? Please help me to answer part b?


title: Variables with different significance levels in linear model (model interpretation)

text: What is the correct interpretation in (multiple) linear regression model when independent variables have different levels of significance? For example, in the regression below the variable $X$ was significant at 0.001, but the intercept was significant just at 0.1. I am not sure if I made myself clear, but saying in other words:   Is the intercept weakening the model?   Should I say the model is significant at 0.1 level (based on the smaller value)?               Call:     lm(formula = Y ~ X, data = df)          Residuals:         Min      1Q  Median      3Q     Max      -2.3182 -0.5247  0.1119  0.4472  1.5385           Coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept)   6.0456     3.1173   1.939   0.0667 .       X             0.6957     0.1169   5.950 8.09e-06 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 0.9465 on 20 degrees of freedom     Multiple R-squared: 0.639,  Adjusted R-squared: 0.6209      F-statistic:  35.4 on 1 and 20 DF,  p-value: 8.087e-06


title: Predict with type='response' for GLM with errorest() function in package ipred

text: Apologies if this is a simple question... I am attempting to use the `errorest` function of the `ipred` package in R to to K-fold CV with GLM models of the binomial family, as well as earth (MARS) models. I have written routines to do CV and can run my GLM and other models through it and it works pretty well. I came across the `errorest()` function and liked it for being a more compact approach than my scripts and the flexibility to work with different models. My problem is that I cannot find a way to have the `predict` function of `errorest` use the flag for `type='response'`. I use this type of prediction because my response variable is presence/absence (1,0). However, I am predicting probabilities, not a binary classification. In the code sample below, I go through a typical GLM and predict with `type='response'`, and then a straight-forward use of `errorest` and finally, a run of `errorest` that calls a custom predict function, `mypredict.glm` that uses the flag `type=response`, however, the results are still not probabilities.               data(mtcars)     require(ipred)          #fit simple GLM     fit <- glm(am ~ mpg + hp, data=mtcars, family='binomial')     #plot results if desired     par(mfrow=c(2,2))     plot(fit)     #predict as type='response'     pred <- predict(fit, newdata=mtcars, type='response')     summary(pred)          # use errorest to do CV and return RMS error and predicted values     errest1 <- errorest(am ~ mpg + hp, data=mtcars, model=glm, estimator="cv",                         est.para=control.errorest(k=5, predictions = TRUE))     errest1$error     summary(errest1$predictions)          #create function for predict flag in errorest(), includes type='response'     mypredict.glm <- function(object, newdata){         predict(object, newdata, type="response")     }          #run errorest with call to predict.glm function     errest2 <- errorest(am ~ mpg + hp, data=mtcars, model=glm,                          predict=mypredict.glm, estimator="cv",                          est.para=control.errorest(k=5, predictions = TRUE))     errest2$error      Variations give:               #predictions are not responses      summary(errest2$predictions)            # with type='response'     summary(pred)          Min.   1st Qu.    Median      Mean   3rd Qu.      Max.      0.0000989 0.0426000 0.2966000 0.4063000 0.6968000 1.0000000           # with errorest() default     summary(errest1$predictions)        Min. 1st Qu.  Median    Mean 3rd Qu.    Max.      -0.3027  0.2273  0.3424  0.3998  0.4437  1.4220           # with errorest and mypredict.glm and type='response'     summary(errest2$predictions)        Min. 1st Qu.  Median    Mean 3rd Qu.    Max.      -0.1941  0.1952  0.3717  0.4056  0.4652  1.3800       Any help would be greatly appreciated. Thank you.


title: Good reference for number of permutations for a simple test?

text: I have two vectors of 100 real values each. I would like to use a permutation test for a few things such as the correlation of the two vectors. I would use the approach of keeping one vector stationary and then permuting the other in relation to it and then calculate the statistic of interest (such as correlation) on the permutation. I have found some guidance for the number iterations for things like bootstrapping, but have not found a general guide for permutations of a simple and similar type. I have code to run the permutation and it is very fast. I am tempted to just run 1,000,000 and have be self-evidently large, but it would be nice to have a reference to justify this, and I haven't found one. Can anyone give direction to one?


title: Rao-Blackwellising state space for a (marginalised) particle filter

text: I am starting to look at particle filtering for a problem that I have. In particular, I would like to reduce the dimensionality of the particles. The model that I have is able to be partitioned. Let my state vector be $\mathbf{x}=\begin{bmatrix}{\mathbf{x}_1}^T && {\mathbf{x}_2}^T \end{bmatrix}^T$. All are continuous variables. My dynamic model may be expressed as: $\begin{bmatrix}{\mathbf{x}_1}(k+1) \\\ {\mathbf{x}_2}(k+1) \end{bmatrix} = \begin{bmatrix}f_1\left( \begin{bmatrix}{\mathbf{x}_1(k)} \\\ {\mathbf{x}_2}(k) \end{bmatrix}\ , \mathbf{u}(k), \mathbf{w}_1(k) \right)\\\ \mathbf{F}_2 \mathbf{x}_2(k) + \mathbf{w}_2(k)\end{bmatrix} \\\ \mathbf{z}(k) = \mathbf{H}_1 \mathbf{x}_1 + \mathbf{H}_2\left(\mathbf{x}_1\right)\mathbf{x}_2 + \mathbf{v}(k)$ where $\left(\bullet\right)$ is "a function of", upper-case bold are matrices and $\mathbf{u}(k)$ is a control vector. Random vectors $\mathbf{w}_1(k)$, $\mathbf{w}_2(k)$ and $\mathbf{v}(k)$ are mutually independent white Gaussian noise. My state space model is not quite the conditionally linear structure as suggested in Gustafsson: $\begin{bmatrix}{\mathbf{x}_1}(k+1) \\\ {\mathbf{x}_2}(k+1) \end{bmatrix} = \begin{bmatrix} f_1\left( \mathbf{x}_1 \right) + \mathbf{F}_1\left( \mathbf{x}_1\right) \mathbf{x}_2 + g_1 \left( \mathbf{x}_1 \right) \mathbf{w}_1(k) \\\ f_2\left( \mathbf{x}_1 \right) + \mathbf{F}_2\left( \mathbf{x}_1\right) \mathbf{x}_2 + g_2 \left( \mathbf{x}_2 \right) \mathbf{w}_2(k) \end{bmatrix} \\\ \mathbf{y}(k) = h_1\left( \mathbf{x}_1 \right) + \mathbf{H}_2\left( \mathbf{x}_1 \right) \mathbf{x}_2 + \mathbf{v}(k)$ I can approximate the state transition with additive noise without too much of an approximation, but I cannot easily split $f_1\left( \begin{bmatrix}{\mathbf{x}_1(k)} \\\ {\mathbf{x}_2}(k) \end{bmatrix}\right)$. Is this model able to Rao-Blackwellised for a particle filter? If so, what papers should I be reading to learn how to do it? * * * **ADDENDUM:** I have managed to split the model in terms of: $\begin{bmatrix}{\mathbf{x}_1}(k+1) \\\ {\mathbf{x}_2}(k+1) \end{bmatrix} = $ $\begin{bmatrix} f_1\left( \mathbf{x}_1(k) , \mathbf{u}(k), \mathbf{w}(k) \right) \\\ \mathbf{F}_2 \left( \mathbf{x}_1(k) \right) \mathbf{x}_2(k) + \mathbf{B}\left(\mathbf{x}_1(k)\right) \mathbf{u}(k) + \mathbf{G}\left( \mathbf{x}_1(k)\right) \mathbf{w}(k) \end{bmatrix} \\\ \mathbf{z}(k) = \mathbf{H}_2\left(\mathbf{x}_1\right)\mathbf{x}_2 + \mathbf{v}(k)$ The (Gaussian) noise vector $\mathbf{w}(k)$ is the same for both the top and bottom parts of the state transition elements, but I can split them if _really_ required at the expense of accuracy. Measurement noise $\mathbf{v}(k)$ is independent. This model is still not conditionally linear in the form above. Is either the original form, or the form above able to be marginalised? * * * **Context:** I have previously used an Unscented Kalman filter for this model, which has worked reasonably well, but I would like to extract a little more performance, if at all possible. If I can partition as above, I can split 18 states into 9 states each for the linear and non-linear parts which, as I understand the theory so far, should drastically reduce the number of particles that I need.


title: Why is the formula for standard error the way it is?

text: So just "why" is $SE = \frac{s}{\sqrt n}$ ? How should one interpret/articulate the reason of having $\sqrt n$ in the denominator. Why do we divide sample mean by the square root of the sample size, intuitively speaking? And how/why is it called standard "error". (Question equally applicable for true standard deviation of the population: $\frac{\sigma}{\sqrt n}$) Is there an intuitive derivation of $SE$ that can make this clear? Please assume you are explaining it to a 6 year old who understands mean and sample size :)


title: MLE for triangle distribution?

text: Is it possible to apply the usual MLE procedure to the triangle distribution? - I am trying but I seem to be blocked at one step or another in the math by the way the distribution is defined. I am trying to use the fact that I know the number of samples above and below c (without knowing c): these 2 numbers are cn and (1-c)n, if n is the total number of samples. However, that does not seem to help in the derivation. The moment of moments gives an estimator for c without much problem. What is the exact nature of the obstruction to MLE here (if indeed there is one)? More details: Let's consider $c$ in $[0,1]$ and the distribution defined on $[0,1]$ by: $f(x;c) = \frac{2x}{c}$ if x < c   $f(x;c) = \frac{2(1-x)}{(1-c)}$ if c <= x Let's take an $n$ i.i.d samples $\\{x_{i}\\}$ from this distribution form the log-likelihood of c given this sample: $\hat{l}(c | \\{x_{i}\\}) = \sum_{i=1}^{n}ln(f(x_{i}|c))$ I am then trying to use the fact that given the form of $f$, we know that $cn$ samples will fall below the (unknown) $c$, and $(1-c)n$ will fall above $c$. IMHO, this allows to decompose the summation in the expression of the log- likelihood thus: $\hat{l}(c | \\{x_{i}\\}) = \sum_{i=1}^{cn}ln\frac{2 x_{i}}{c} + \sum_{i=1}^{(1-c)n}ln\frac{2(1-x_{i})}{1-c}$ Here, I am unsure how to proceed. MLE will involve taking a derivative w.r.t. $c$ of the log-likelihood, but I have $c$ as the upper bound of the summation, which seems to block that. I could try with another form of the log- likelihood, using indicator functions: $\hat{l}(c | \\{x_{i}\\}) = \sum_{i=1}^{n}\\{x_{i}<c\\}ln\frac{2 x_{i}}{c} + \sum_{i=1}^{n}\\{c<=x_{i}\\}ln\frac{2(1-x_{i})}{1-c}$ But deriving the indicators doesn't seem easy either, although Dirac deltas could allow to continue (while still having indicators, since we need to derive products). So, here I am blocked in MLE. Any idea?


title: Comparing regression coefficients across models with standardized dependent variables

text: **The Situation:** There are four *** _identical_ *** spatial regression models, except each uses a different dependent variable. The independent variables consist of a standard set of variables derived from a principal component analysis. The dependent variables have been standardized by using the standard score for each observation. **The Question:** Can the regression coefficients be directly compared across the models? For example, can I say that because the coefficient for an IV is .25 in a model using the first DV and .50 in a model using the second DV, that the impact of the IV in the first model is twice that of the impact of the second model? **A little more background:** I am using four spatial spatial regression models to to look at what socio- demographic factors are associated with poverty in an area. The units of observation are census tracts. The four models are identical, except for the dependent variable used. Each model has the same set of independent variables from the same location and data set. The independent variables were derived using factors from a principal component analysis. The dependent variables used represent four different ways of measuring the well-being of a census tract, which are: 1) the percentage of poor people per census tract using poverty line A, 2) the percentage of poor people per census tract using poverty line B, 3) The Per Capita Income of the census tract, and 4) A ratio of census tract income to average regional income. Instead of using the observations for the DVs however, I have used the standard scores (z scores) of each observation, with the intention of making the regression results comparable across models. **To further complicate things:** Two of the variables directly measure poverty, while the other two measure income, so the direction of the relationships between the IVs and DVs has been reversed. Would this cause any additional problems in directly comparing regression coefficients?


title: Type of a variable built by binning some numeric variable

text: I have some sample data with age and I'd like to put them into bins (for example, ranges 20-24 25-30, etc.). The resulting variable would constitute values from 1-8. What I'm confused about is what type of variable is the new variable?


title: Acceptable scores on a survey

text: I want to know what per cent is acceptable in a training survey? If we get scores of 77% satisfaction of our materials is that acceptable, good, poor?


title: Proving that MGF determines PDF when the PDF is defined for whole real line

text: > If two PDFs have the same moment generating function that converges in an > open set around 0, then the PDFs are same. This is a well known fact, but I can't find its proof. If the PDFs are defined for only non-negative values, the MGF is essentially Lapleace transform and the uniquesness problem is just that of Laplace transform. However, if the PDFs are defined on the whole real line, the MGF is not Laplace transform. I guess the condition that the MDF converges around zero may take care of something in this case. Could someone provide the proof or reference for it?


title: Standard normal distribution - the z-table for the PDF

text: I have searched high and low, and I simply cannot find this table anywhere. I am quite sure that for practical purposes I don't need it - that there are other methods, but it would be extremely useful if anybody could point me in the right direction of finding one online. Also, it needs to have negative values. I have a funny feeling I'm going to get a ribbing for asking this. Have mercy.


title: What are principal component scores?

text: I know this is probably simplistic but what are Principal component scores? This question originates from my attempt to understand this question here.


title: Learning about a population mean

text: ![enter image description here](http://i.stack.imgur.com/HEZRS.jpg)   1. Assuming $H_0$ is true, what is the distribution of the test statistic? t(29)   2. Assuming $H_0$ is true, what is the expected value of the test statistic? 87.70   3. The sample mean of 87.7 was ___ standard errors below the mean of 90. 1.615 Can someone verify that my work is correct? In particular, I'm not entirely sure about the second question. Will the test statistic actually be 0 in this case (as it wouldn't deviate from the mean)?


title: Fitting a probability distribution to non i.i.d. data?

text: I have temperature time series data that I have determined is not independently and identically distributed (from looking at the autocorrelation plots and Ljung-Box tests). However, I am still able to fit PDFs, typically GEV, gamma, and Weibull distributions to my different data sets, with a $\chi^2$ goodness-of-fit test accepting the null hypothesis. So my question is: since these underlying PDFs assume the data is supposed to be i.i.d., and my data is not i.i.d., how should I proceed? How worried should I be about my estimates being biased?


title: Test overfitting of logistic regression with limited volume

text: I have a set of samples with two labels red and black. I can build a logistic regression model to predict the label colour. Once a model is built, I would like to test whether it is overfitting or not. Normally, I will set aside, say, 30% of my sample as out-of sample. Build a logistic regression model on the 70% development sample and test its performance (for example, Gini) on the out-of sample. If I see Gini drops from development sample to out-of sample, I know that the model may be overfitted. However, when my sample has small number of red (or black or even both), it is reluctant to set aside some out-of sample for validation purpose. I'd rather use as much data as possible given a limitation. So what are some effective validation methods other than the one I describe above that can be used? Note that here I am not trying to determine the best model form (as logistic regression is the one to use), or to determine best value of some tuning parameters (I don't need to determine the number of variables to use). I simply have got a model and want to test whether it is overfitting the training sample or not. So I don't think cross-validation would be applicable here.


title: Empirical distribution alternative

text: ## BOUNTY: The full bounty will be awarded to someone who provides a reference to any published paper which uses or mentions the estimator $\tilde{F}$ below. ## Motivation: This section is probably not important to you and I suspect it won't help you get the bounty, but since someone asked about the motivation, here's what I'm working on. I am working on a statistical graph theory problem. The standard dense graph limiting object $W : [0,1]^2 \to [0,1]$ is a symmetric function in the sense that $W(u,v) = W(v,u)$. Sampling a graph on $n$ vertices can be thought of as sampling $n$ uniform values on the unit interval ($U_i$ for $i = 1, \dots, n$) and then the probability of an edge $(i,j)$ is $W(U_i, U_j)$. Let the resulting adjacency matrix be called $A$. We can treat $W$ as a density $f = W / \iint W$ supposing that $\iint W > 0$. If we estimate $f$ based on $A$ without any constraints to $f$, then we cannot get a consistent estimate. I found an interesting result about consistently estimating $f$ when $f$ comes from a constrained set of possible functions. From this estimator and $\sum A$, we can estimate $W$. Unfortunately, the method that I found shows consistency when we sample from the distribution with density $f$. The way $A$ is constructed requires that I sample a grid of points (as opposed to taking draws from the original $f$). In this stats.SE question, I'm asking for the 1 dimensional (simpler) problem of what happens when we can only sample sample Bernoullis on a grid like this rather than actually sampling from the distribution directly. references for graph limits: L. Lovasz and B. Szegedy. Limits of dense graph sequences (arxiv). C. Borgs, J. Chayes, L. Lovasz, V. Sos, and K. Vesztergombi. Convergent sequences of dense graphs i: Subgraph frequencies, metric properties and testing. (arxiv). ## Notation: Consider a continuous distribution with cdf $F$ and pdf $f$ which has a positive support on the interval $[0,1]$. Suppose $f$ has no pointmass, $F$ is everywhere differentiable, and also that $\sup_{z \in [0,1]} f(z) = c < \infty$ is the supremum of $f$ on the interval $[0,1]$. Let $X \sim F$ mean that the random variable $X$ is sampled from the distribution $F$. $U_i$ are iid uniform random variables on $[0,1]$. ## Problem set up: Often, we can let $X_1, \dots, X_n$ be random variables with distribution $F$ and work with the usual empirical distribution function as $$\hat{F}_n(t) = \frac{1}{n} \sum_{i=1}^n I\\{X_i \leq t\\}$$ where $I$ is the indicator function. Note that this empirical distribution $\hat{F}_n(t)$ is itself random (where $t$ is fixed). Unfortunately, I am not able to draw samples directly from $F$. However, I know that $f$ has positive support only on $[0,1]$, and I can generate random variables $Y_1, \dots, Y_n$ where $Y_i$ is a random variable with a Bernoulli distribution with probability of success $$p_i = f((i-1+U_i)/n)/c$$ where the $c$ and $U_i$ are defined above. So, $Y_i \sim \text{Bern}(p_i)$. One obvious way that I might estimate $F$ from these $Y_i$ values is by taking $$\tilde{F}_n(t) = \frac{1}{\sum_{i=1}^n Y_i} \sum_{i=1}^{\lceil tn \rceil} Y_i$$ where $\lceil \cdot \rceil$ is the ceiling function (that is, just round up to the nearest integer), and redraw if $\sum_{i=1}^n Y_i = 0$ (to avoid dividing by zero and making the universe collapse). Note that $\tilde{F}(t)$ is also a random variable since the $Y_i$ are random variables. ## Questions: From (what I think should be) easiest to hardest.   1. Does anyone know if this $\tilde{F}_n$ (or something similar) has a name? Can you provide a reference where I can see some of its properties?   2. As $n \to \infty$, is $\tilde{F}_n(t)$ a consistent estimator of $F(t)$ (and can you prove it)?   3. What is the limiting distribution of $\tilde{F}_n(t)$ as $n \to \infty$?   4. Ideally, I'd like to bound the following as a function of $n$ -- e.g., $O_P(\log(n) /\sqrt{n})$, but I don't know what the truth is. The $O_P$ stands for Big O in probability $$ \sup_{C \subset [0,1]} \int_C |\tilde{F}_n(t) - F(t)| \, dt $$ ## Some ideas and notes:   1. This looks a lot like acceptance-rejection sampling with a grid-based stratification. Note that it is not though because there we do not draw another sample if we reject the proposal.   2. I'm pretty sure this $\tilde{F}_n$ is biased. I think the alternative $$\tilde{F^*}_n(t) = \frac{c}{n} \sum_{i=1}^{\lceil tn \rceil} Y_i$$ is unbiased, but it has the unpleasant property that $\mathbb{P}\left(\tilde{F^*}(1) = 1\right) < 1$.   3. I'm interested in using $\tilde{F}_n$ as a plug-in estimator. I don't think this is useful information, but maybe you know of some reason why it might be. ## Example in R Here is some R code if you want to compare the empirical distribution with $\tilde{F}_n$. Sorry some of the indentation is wrong... I don't see how to fix that.               # sample from a beta distribution with parameters a and b     a <- 4 # make this > 1 to get the mode right     b <- 1.1 # make this > 1 to get the mode right     qD <- function(x){qbeta(x, a, b)} # inverse     dD <- function(x){dbeta(x, a, b)} # density     pD <- function(x){pbeta(x, a, b)} # cdf     mD <- dbeta((a-1)/(a+b-2), a, b) # maximum value sup_z f(z)               # draw samples for the empirical distribution and \tilde{F}     draw <- function(n){ # n is the number of observations       u <- sort(runif(n))        x <- qD(u) # samples for empirical dist       z <- 0 # keep track of how many y_i == 1       # take bernoulli samples at the points s       s <- seq(0,1-1/n,length=n) + runif(n,0,1/n)        p <- dD(s) # density at s       while(z == 0){ # make sure we get at least one y_i == 1         y <- rbinom(rep(1,n), 1, p/mD) # y_i that we sampled         z <- sum(y)       }       result <- list(x=x, y=y, z=z)       return(result)     }          sim <- function(simdat, n, w){       # F hat -- empirical dist at w       fh <- mean(simdat$x < w)            # F tilde           ft <- sum(simdat$y[1:ceiling(n*w)])/simdat$z           # Uncomment this if we want an unbiased estimate.           # This can take on values > 1 which is undesirable for a cdf.           ### ft <- sum(simdat$y[1:ceiling(n*w)]) * (mD / n)       return(c(fh, ft))     }               set.seed(1) # for reproducibility          n <- 50 # number observations     w <- 0.5555 # some value to test this at (called t above)     reps <- 1000 # look at this many values of Fhat(w) and Ftilde(w)     # simulate this data     samps <- replicate(reps, sim(draw(n), n, w))          # compare the true value to the empirical means     pD(w) # the truth      apply(samps, 1, mean) # sample mean of (Fhat(w), Ftilde(w))     apply(samps, 1, var)  # sample variance of (Fhat(w), Ftilde(w))     apply((samps - pD(w))^2, 1, mean) # variance around truth               # now lets look at what a single realization might look like     dat <- draw(n)     plot(NA, xlim=0:1, ylim=0:1, xlab="t", ylab="empirical cdf",          main="comparing ECDF (red), Ftilde (blue), true CDF (black)")     s <- seq(0,1,length=1000)     lines(s, pD(s), lwd=3) # truth in black     abline(h=0:1)     lines(c(0,rep(dat$x,each=2),Inf),              rep(seq(0,1,length=n+1),each=2),              col="red")         lines(c(0,rep(which(dat$y==1)/n, each=2),1),           rep(seq(0,1,length=dat$z+1),each=2),           col="blue")      ![output from the above data](http://i.stack.imgur.com/BQ4ac.png) ## EDITS: EDIT 1 -- I edited this to address @whuber's comments. EDIT 2 -- I added R code and cleaned it up a bit more. I changed notation slightly for readability, but it is essentially the same. I'm planning on putting a bounty on this as soon as I'm allowed to, so please let me know if you want further clarifications. EDIT 3 -- I think I addressed @cardinal's remarks. I fixed the typos in the total variation. I'm adding a bounty. EDIT 4 -- Added a "motivation" section for @cardinal.


title: Maximizing F1-measure, when you have an algorithm for minimizing another loss function

text: Let's say you have an algorithm for minimizing the following loss function: $$ loss = \sum_i l(y_i, f(x_i)) $$ Let's say you are in the binary classification case, and the ratio of negative to positive instances is $r$, where $r >> 1$ (say $r > 20$). Given this algorithm and this setting, we want to find a model which maximizes F1 measures(or at least approximates it). For simplicity let's say we are using squared hinge loss : $\max(0, 1- y_if(x_i))^2$. The problem is if you just train with the loss function, the resulting classifier will almost always give you the negative answer, which will give a very poor F1. \-- Option 1: Train on a weighted loss function: $$ loss = \sum_{i; y_i=+1} l(y_i, f(x_i)) + \alpha \sum_{i; y_i=-1} l(y_i, f(x_i)) $$ for some $\alpha < 1$ Choose the value of $\alpha$ which maximizes F1, on the validation set. This is equivalent to using only $\alpha$ portion of the samples negative label. \-- Option 2 (or plus option 1): Threshold adjustment: Prior to this, prediction was positive if $f(x_i) > 0$, and negative if otherwise. After training the model on the training set, choose $\theta$ on the validation set, as new threshold: output is positive, if $f(x_i) > \theta$, negative if otherwise. Choose $\theta$ such that it maximizes F1 on the validation set. Question: Do you see any alternative solution? Note: No! I can't change the loss function into negative F1, unless I can find a differentiable approximation.


title: chi squared test for various locus

text: I am trying to compare the observed number of DNA sequencing reads versus the expected number. There are loci where I would expect all the sequencing reads to be the same and there are other loci where I would expect them to be of 2 types, each one with an expected frequency. For example:               LOCUS 1     OBSERVED READS    EXPECTED FREQ     type_1 type_2     type_1 type_2     400    0          100%   0%     400    40         100%   0%     222    644        41%    59%     435    154        100%   0%      I thought using the chi square test would fit the test. In this way I calculated the expected read counts and the chi square value, but I am confused as to what would be my degrees of freedom. I dont think it would be 1 but I am stuck. Also how come in the example locus in the 4th row the chi square value between the observed and the expected is 0 (since expected is zero). Thanks so much for all your help


title: Applying ridge regression for an underdetermined system of equations?

text: When $y = X\beta + e$, the least squares problem which imposes a spherical restriction $\delta$ on the value of $\beta$ can be written as \begin{equation} \begin{array} &\operatorname{min}\ \| y - X\beta \|^2_2 \\\ \operatorname{s.t.}\ \ \|\beta\|^2_2 \le \delta^2 \end{array} \end{equation} for an overdetermined system. $\|\cdot\|_2$ is the Euclidean norm of a vector. The corresponding solution to $\beta$ is given by \begin{equation} \hat{\beta} = \left(X^TX + \lambda I\right)^{-1}X^T y \ , \end{equation} which can be derived from the method of Lagrange multipliers ($\lambda$ is the multiplier): \begin{equation} \mathcal{L}(\beta,\lambda) = \|y-X\beta\|^2_2 + \lambda(\|\beta\|^2_2 - \delta^2) \end{equation} I understand that there is a property that \begin{equation} \left(X^TX + \lambda I\right)^{-1}X^T = X^T\left(XX^T + \lambda I\right)^{-1} \ . \end{equation} The right hand side resembles the pseudo-inverse of the regressor matrix $X$ in the underdetermined case (with the added regularization parameter, $\lambda$). Does this mean mean that the same expression can be used to approximate $\beta$ for the underdetermined case? Is there a separate derivation for the corresponding expression in the underdetermined case, as the spherical restriction constraint is redundant with the objective function (minimum norm of $\beta$): \begin{equation} \begin{array} &\operatorname{min.}\ \| \beta \|_2\\\ \operatorname{s.t.}\ X\beta = y \ . \end{array} \end{equation}


title: Population mean and hypothesis testing

text: > The national average price of a gallon of regular gasoline is 3.71 dollars. > Linda would like to assess if the average gas price in her city is > significantly higher than the national average. She samples gas prices at 20 > local gas stations and records the price of a gallon of regular gasoline. > Assume her sample of prices can be considered a random sample > (representative of the larger population of all such prices). > > What additional condition is required to be able to perform the one-sample > t-test? I know one condition is that the sample is random, but I'm having a hard time figuring out what the second should be. Should the distribution be normal? If so, can I just construct a Q-Q plot to test for this? Any help would be greatly appreciated.


title: Goodness of fit of a model

text: My question might not be making much sense but I have been stuck for few days already so really do hope to achieve a solution on here. So basically I have simulated a data set from a power-law distribution with a known parameter, I would like to test which method (binning methods+MLE) performs the best on estimating the parameter of the model. Power-law distribution: f(x)=Cx^(-\mu) Here is how I simulated my model in R:               u<-runif(500)     fakedata500<-((2*(1-u))^(-1))      and fakedata500 is the data set I want to test. I have let the parameter of the distribution to be \mu =2. First method is linear binning, plot the frequency graph and do linear regression after:               a<-hist(fakedata500, breaks=c(0,2,4,8,16,32,64,128,256),plot=FALSE)     > x<-c(0,2,4,8,16,32,64,128,256)     > y<-c(a$counts,0)     > x<-cbind(x,y)     > x.df<-data.frame(x)     > x.sub<-subset(x.df,y>0 & x>0)     > u<-log10(x.sub)     > plot(u, main=expression(paste("500 random numbers generated, " ,mu, " =2.00" )),     xlab="Log (geometric bin width)", ylab="Log (Frequency x)")          > res=lm(u)     > res          Call:     lm(formula = u)          Coefficients:     (Intercept)            y             2.187       -1.060        Second method is Cumulative distribution and linear regression:                load("fakedata500.Rda")     > X<-fakedata500     > p<-ppoints(100)     > plot(log(quantile(X,p=p)),log(1-p),     +       ylab="Log (1-sum of values < x)",xlab="Log (x)",main=expression(paste("500     random numbers generated, " ,mu, " =2.00" )))     > df  <- data.frame(x=log(1-p),y=log(quantile(X,p=p)))     > fit <- lm(y~x,df)     > fit          Call:     lm(formula = y ~ x, data = df)          Coefficients:     (Intercept)            x           -0.6728      -0.9775        Third method is Maximum likelihood estimation (MLE):               > load("fakedata500.Rda")     > library(stats4)     > library(bbmle)     > x<-fakedata500     > pl <- function(u){-length(x)*log(u-1)-length(x)*(u-1)*log(0.5)+u*sum(log(x))}     > mle1<-mle2(pl, start=list(u=2), data=list(x))     > summary(mle1)     Maximum likelihood estimation          Call:     mle2(minuslogl = pl, start = list(u = 2), data = list(x))          Coefficients:       Estimate Std. Error z value     Pr(z)         u 2.001598   0.044793  44.686 < 2.2e-16 ***     ---     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1          -2 log L: 1303.66       Is likelihood ratio test/g-test/goodness of fit test a appropriate method to compare and look for the best method to estimate the value of the parameter? Or g-test is only for testing the MLE one? If so, how to perform g-test/goodness of fit test in R on testing the validity of the MLE? Sorry for this long post and thank you in advance for any help.


title: In non-negative matrix tri-factorization, initialization not possible because matrix is singular

text: I have implemented the non-negative matrix tri-factorization algorithm (link to paper). If is similar to the more widely known NMF (non-negative matrix factorization), but incorporates _prior knowledge_ to represent semi- supervised learning. The basic tenet is that a term-document matrix X can be factorized as $$X = FSG^T$$ where   * $X$ is an $m\times n$ term-document matrix representing $n$ documents and $m$ terms   * $F$ is an $m\times k$ matrix. i'th row represents probability of $t_i$ belonging to $k$ classes   * $G$ is an $n\times k$ matrix. i'th row represents probability of $d_i$ belonging to $k$ classes   * $S$ is a $k\times k$ matrix, provides a low-dimensional condensed view of $X$. The matrix S is initialized as $S = (F^TF)^{-1}F^TXG(G^TG)^{-1}$. I ran my code, which strictly follows the above paper, and it ran fine on one dataset. On another dataset, however, the matrix F leads to a situation where $F^TF$ is not invertible. As a result, I can't initialize the matrix S. I have literally spent days trying to find some information about how to proceed in this scenario, but none of the related papers mention anything. If anyone knows how to handle this situation, please help me out. Thanks!


title: Is $\bar{h}$ a typo of $h$ in Functional Equivariance and Formal Invariance in Lehmann's TPE?

text: On page 161 of Lehmann and Casella's Theory of Point Estimation, they introduce Functional Equivariance and Formal Invariance. $\bar{h}$ appears in two places:   * in "Functional Equivariance":  > $\phi[\delta(X)]$as the estimator of $\phi[\bar{h}(\theta)]$   * in "Formal Invariance": equation (2.16) I didn't find the definition of $\bar{h}$, so I was wondering if $\bar{h}$ is a typo of $h$? Or how is it defined? Thanks and regards! ![enter image description here](http://i.stack.imgur.com/0K17l.png)


title: interval or ratio?

text: I have a question about how to interpret some discrete data. Some are very easy, like the number of children in a household: that is a discrete and numerical, measured as a ratio variable. Other variables are not that easy, for example   * age, measured in whole years: is that interval or ratio?   * number of cigarettes smoked per day: is that interval, or ratio? The difference between interval and ratio is that ratio has an absolute zero (when the measured property is absent). In case of the age measured in whole years, an infant aged 2 months has zero (whole) years, but that does not mean the person doesn't have an age... so is that ratio?


title: Alternatives to constrained non-linear optimization: a case study with Gaussian curve fitting

text: I have a set of data points distributed like a bell-shape curve with varying amplitudes and widths. I need to parameterize these data points and a Gaussian with amplitude, width and offset parameters function is just perfectly suited. My y-values are defined between at 8 different points between [-2.3561 3.1415] and I am interested fitting a Gaussian that is centered at 0.               out      = amp.*exp( - (0.5*(X./sd).^2 )) + offset;      Optimizing the likelihood function following an initial estimation of these 3 parameters works generally very good. And this has been already discussed elsewhere in this site. The problem starts when my data points are described by a very shallow bell- shaped profile, looking like a second order polynomial with a small coefficient, so that the highest and lowest points do not have so much of a difference. In these cases the optimization procedure returns extremely high sigma and amplitude values. It basically tries to fit the data with the tip of an amazingly wide Gaussian, which by far exceeds the range of my x-values. These huge fitted values destroy my statistics, because ultimately I would like to do some statistics with these numbers comparing different groups. In these situations, I would like to have a small amplitude and relatively wide width parameter returned from the optimization procedure and this would just do the job good enough. I tried constrained optimization, trying to limit valid values for the parameters. But this also results in values that are densely accumulated on the boundary limits. I was wondering what would be a general approach to tackle this kind of situations? If you have encountered similar situation I would be happy to hear about your experience... Here are some examples of data: Gray Curve: Initial Estimates Red Curve: Fitted Gaussian Blue Curve: Data data_01: This shows a situation where sigma and amp parameters are huge, even though a very small amplitude value would do the job good. `-2.3562 1.2210 -1.5708 -1.6224 -0.7854 -0.1330 0 -0.1126 0.7854 -0.0745 1.5708 1.8017 2.3562 -0.1533 3.1416 -0.9268` ![](http://i.stack.imgur.com/AWkwl.png)


title: Best model for longitudinal analysis

text: I have one dependent variable and multiple predictors. Dependent is continuous and predictors may be continuous or dichotomous. My sample is of 200 individuals. Each month, for 6 months, I have a measure of predictors for the 200 individuals and I would like to identify those which have a significant effect over the dependent variable. My problem is that I only have one yearly measure for the dependent variable. This variable is a judgment about the service provided and is given once a year. Is there a good model you would suggest to take into account repeated measures over time (independent variables) and their influence on a single measured dependent variable?


title: Confusion matrix for events in a time series

text: I'm not big on statistics, so please excuse my ignorance. I have a video recording that I want to evaluate, I have an algorithm that can transform this video into a time series where I have 0 everywhere except for a couple of frames where an event (A) occurs. Then I have manual annotations of another event (B) that I think is related to the 1st event (A) (i.e. event B appears shortly after event A). I want to construct a confusion matrix like this:   1. Event A is condition, event B is "test"   2. True positive value is when there is an event A and within 50 frames there is also event B   3. False positive is when A is not present, but B is   4. False negative is when A is present, but B is not Now my problem is about True Negatives. If you do vaccine testing, you have a total number of tests, and you can quantify True Negatives without issue. But what about my case? I either have event A or B but True Negative is by definition everything else in the time series? Does it make sense to even use confusion matrix?


title: Posterior distribution of precision for multivariate normal with normal-wishart prior

text: I'm trying to derive the posterior distribution for the precision matrix for the multivariate normal with normal-wishart prior. According to wikipedia and other sources the answer is as follows: $p(S|\mu, X, W, v) \\\ \quad \sim W((W^{-1} + \sum_{i=1}^{N}(x_i-\bar{x})(x_i-\bar{x})^{T} + \frac{rN}{r+N}(\mu_0-\bar{x})(\mu_0-\bar{x})^{T})^{-1}, v+N)$ When I try to derive it myself I see how I could get the result but it confuses me since I get the result as $p(S|\mu, X, W, v) \\\ \quad\sim W( (W^{-1} + \sum_{i=1}^{N}(x_i-\mu)(x_i-\mu)^{T} + \frac{rN}{r+N}(\mu-\mu_0)(\mu-\mu_0)^{T})^{-1}, v+N+1)$ I come to this result by multiplying the likelihood, prior mean and prior precision. $\mu_0$ being the hyper prior of the mean for $\mu$ I'm thankful for any advice on how to derive the posterior correctly and which step I'm doing wrong. I'm sorry but I don't know how to get the equations presented.


title: Pooling asymmetric confidence intervals for proportions?

text: I have several measurements of proportions (values in [0, 1]) $\theta_1,...,\theta_n$, each with an (asymmetric) 95% confidence interval. The $\theta$'s are repeated measurements of the same variable in the world and so I'd like to pool them together to get an estimate of $\hat{\theta}$ with its own confidence intervals (which can be asymmetric.) The important thing is that I'd like the estimate to downweight measurements in proportion to the width of their confidence interval. If I have $\theta_{1} = 0.6$ with CIs $(0.5,0.7)$ and $\theta_{2} = 0.8$ with CIs $(0.3,0.9)$ then I'd like the estimate $\hat{\theta}$ to weight $\theta_{1}$ more highly and thus give an estimate closer to 0.6 than to 0.8 (but also respect the asymmetry of $\theta_{2}$'s confidence intervals.) What's a general way to do this? One simple minded way to get $\hat{\theta}$ is to use the size of each CI in a weighted sum: $s_1 = 0.7 - 0.5, s_2 = 0.9 - 0.3$ $\hat{\theta} = \frac{\theta_{1}(1/s_{1}) + \theta_{2}(1/s_{2})}{(1/s_{1}) + (1/s_{2})}$ But that doesn't yield a CI for $\hat{\theta}$ and does not use the asymmetry of the CIs of $\theta_{1},\theta_{2}$, it only uses their size. How can this be done correctly? Does one have to transform the proportions to a variable that can be modeled by a more well-behaved distribution like a Gaussian? Final variant of this question: each proportion $\theta_i$ was obtained by monte carlo sampling. Each $\theta_i$ is just the mean value of that distribution (which was a posterior distribution for $\theta_{i}$) and the 95% intervals are derived from the samples used to estimate the distribution. That means I have more information about each $\theta_i$ than just the posterior mean and CIs; I actually have an entire distribution of values for it and so are there general methods to use these here as well? thanks. p.s. I'm happy with the most minimal assumptions about how the $\theta$'s are related to each other. Theoretically, they are all measurements of the same variable but each measurement has distinct amounts of noise. I have no strong priors on the noise process.


title: How to detect outliers in randomForest regression models?

text: The randomForest package in R software includes outlier function for the detection of outliers. This function uses proximity matrix or randomForest object for the outlier detection. The manual says that the type of the randomForest object can not be regression? Why is it so that this function can be used for classification models only? Is there a reasonable way to detect outliers in regression models?


title: Call volume: time series regression model from 52 weeks a year and lagged predictors

text: I have been adamantly searching the web to learn how to successfully implement a dynamic regression time series in the forecast package for R. The time series data that I am using is weekly data (frequency=52) of incoming call volume and prediction variables are mailers sent out every now and then. They are a significant predictor of the data for the week that they hit, the following week, and the week after that. I have created lagged variables and use these three as the predictors. My main concern is that the arima model is not taking into account the time series frequency. When I tell it to recognize the ts with a frequency of 52 it has an error. I have looked at the fortrain function but do not understand it. I also have looked at the tbats suggested but found that those will not work with prediction variables. The Zoo function recognizes 52 frequency but it is not advised to use with the forecast package. Here is the basic code. The problem is that the time series calwater[,5] is not recognized as such. It is imputed as a simple vector as an integer...               #this works without taking into acount the ts     fit2 <- auto.arima(calwater[6:96,5], xreg=calwater[6:96,6:8], d=0)     fccal <- forecast(fit2, xreg=calwater[97:106,6:8], h=10)     fccal     plot(fccal, main="Forecast Cal Water", ylab="Calls")          #to form a ts object     calincall<-ts(calwater[1:106,5],start=c(2011,23),frequency=52)          #once the ts is added to the model this dispalys     #Error in `[.default`(calincall, 2:100, 1) : incorrect number of dimensions      Maybe the error is because there is just a little over two years of data.               #Time Series: Start = c(2011, 23), End = c(2013, 24),Frequency = 52       I would be very grateful for any guidance in for this particular issue. I am using the forecast package and prefer to continue within the package but I am open to suggestions.


title: Order statistic expectations and the variance of the original distribution

text: Is there a proof or theorem out there that says that the difference between the expectations of the order distributions of $n$ draws from a distribution $X$ decreases as the variance of $X$ decreases? Specifically, I am asking: If X and Y are probability distributions defined over the same range and with the same mean, does Variance(X) < Variance (Y) imply that $\mathbb{E}X_{kn} - \mathbb{E}X_{sn} < \mathbb{E}Y_{kn} - \mathbb{E}Y_{sn}$ where $X_{in}$ is the $i_{th}$ greatest of $n$ draws from $X$, and $k>s$. In fact, I can be even more specific, as I am interested in the difference between the expectation of the mean and the extreme values. So, does $Var(X) < Var(Y) \to \mathbb{E}X - \mathbb{E}X_{1n} < \mathbb{E}Y - \mathbb{E}Y_{1n}$ (and the same for the absolute value of the difference between EX and EXnn) Also, for the general case: would it be true that for any distributions X and Y defined over the same range, that Variance(X) < Variance (Y) implies $\mathbb{E}X_{kn} - \mathbb{E}X_{sn} < \mathbb{E}Y_{kn} - \mathbb{E}Y_{sn}$


title: What descriptive statistics should be reported in tables and graphs when using Friedman's nonparametric test?

text: I'm using Friedman's nonparametric tests for repeated measures for my thesis (I tried transformations on all variables first but data was still highly skewed) and had 2 questions:   1. What are the best descriptive stats to put in a table. Means and SDs are commonly reported but Friedman's tests use rank ordering of the means. Should I report median values and range? I read somewhere that this is more appropriate but I've never really seen it done in empirical papers. Or is it better to report mean, SDs and medians?   2. Second question is what statistics to plot in a graph. I had means with standard error bars but I'm also not sure if this is ok, given the use of nonparametric tests. A previous post said confidence intervals can be used instead of standard error bars with median values, but I'd really just like to know if it's ok to use a graph with means and SDs when your analysis contains nonparametric tests, or whether this is totally unacceptable.  Thanks in advance for any advice anyone can provide!


title: Are both of these generalized additive models?

text: I've been reading about generalised additive models. I've been using this data (which is a reformatted version of data from Hastie's website), and running my code in R. This data essentially consists of whether a patient has or does not have coronary heart disease (the variable to be modelled) and several patient characteristics. For the sake of this question, the two I am interested in are `sbp` (sytolic blood pressure) and `chol.ratio` (the ratio of two types of cholesterol). I've been trying to model chd with logistic regression. Looking at scatterplots `sbp` seems non-linearly related to the logit, so I've modelled it with splines as follows: Using `gam`:               data1 <- read.csv("coris.csv", sep = ",", stringsAsFactors = F,                                 header = T)     require(gam)       gam.object <- gam(chd ~ s(sbp, 5) + chol.ratio, data = data1,                               family = binomial(link = "logit"))      And using the `Design` package:               require(design)     rcs.object <- lrm(chd ~ rcs(sbp, 6) + chol.ratio, data = data1)      I can plot the subsequent models in 3d using the rcl() package, and the outputs are very similar - the model using GAM is: ![enter image description here](http://i.stack.imgur.com/kP5wv.png) The plane represents the fitted model over a range of the predictor variables, and the points are the actual fitted model, the z axis to the right is the cholesterol ratio, and the x axis to the left is the sytolic blood pressure, and the vertical axis is the logit. and the model using lrm with a rcs is : ![enter image description here](http://i.stack.imgur.com/he5LZ.png) So - with the lrm command are you actually fitting a generalized additive model if you specify a spline with rcs()? If not, why not, and how do the two approaches differ?


title: $X_{n}$ converges in probability where $X_{n}=\frac{Y_{n}}{n}$

text: Let $Y_{n}$ be a sequence of independent and identically distributed random variables and $X_{n}=\frac{Y_{n}}{n}$ . show that $X_{n}$ converges in probability. decide whether $X_{n}$ converges a.e or not. thanks for help.


title: Why would one use age-squared as a covariate in a genetic association study?

text: Why would one use age and age-squared as covariates in a genetic association study? I can understand the use of age if it has been identified as a significant covariate, but I am at a loss as to the use of age-squared?! Thanks, Kevin


title: What does k mean in AIC expression?

text: I'm using the `AIC()` function in R, where the argument `k` (default to 2) stands for the number of parameters. I'm not sure to understand well what does it mean. How can you know k? For example, in the case of $y=ax^2+bx+c$, does k equal 3? As a sidenote, I saw that there is a function named `sic()` which comes with the `metomet` package. However, I can't install this package as it doesn't seem to be available for latest release of R (version 2.13.0). Do you know where it can be found?


title: Best way to put two histograms on same scale?

text: Let's say I have two distributions I want to compare in detail, i.e. in a way that makes shape, scale and shift easily visible. One good way to do this is to plot a histogram for each distribution, put them on the same X scale, and stack one underneath the other. When doing this, how should binning be done? Should both histograms use the same bin boundaries even if one distribution is much more dispersed than the other, as in Image 1 below? Should binning be done independently for each histogram before zooming, as in Image 2 below? Is there even a good rule of thumb on this? ![Image 1](http://i.stack.imgur.com/Kqzdr.png) ![Image 2](http://i.stack.imgur.com/BpPwP.png)


title: Data snooping and multiple comparison in ANOVA

text: From Kutner's Applied Linear Statistical Models > The Bonferroni multiple comparison procedure does not lend itself to data > snooping, unless one can specify in advance the family of inferences in > which one may be interested and provided this family is not large. > > The Tukey and Scheffe procedures involve families of inferences that lend > themselves naturally to data snooping. the Tukey and Scheffe procedures, > allow data snooping to be undertaken naturally without affecting the > confidence coefficient or significance level. Why is Bonferroni procedure not working well in data snooping problem, while the Tukey and Scheffe procedures can? How shall I understand the reason of "unless one can specify in advance the family of inferences in which one may be interested and provided this family is not large"? Thanks!


title: Which metric should I trust to evaluate my predictive model

text: I am working on predictive model and when I evaluate it, I find good `accuracy_score`, `precision_score`, `recall_score`, and `f1_score`. But I don't get good results using confusion matrix. What is the metric that should I Trust ?               accuracy_score:  0.987727165447     confusion_matrix: [[338961   3857 ]                        [   252     27  ]]     precision_score: 0.998449790932     recall_score:    0.987727165447     f1_score:        0.993031083087      I am using Scikit-Learn metrics.               ‘precision’ sklearn.metrics.precision_score     ‘recall’    sklearn.metrics.recall_score     ‘f1’        sklearn.metrics.f1_score


title: How can I add the values of some variables given in one data set by the nature of the variables given in another data sets?

text: Please consider two data sets dat1 and dat2 in the following format.               Dat1:           x1  x2  x3  x4  x5      1    3   8   2  10      2   12   2   9   6      and Dat2:               var     type     x1    discrete     x2     contns     x3    discrete     x4     discrete      x5     contns      I want to add the values of the variables in dat1 by the nature from dat2. In other words, I want to have variables adding values from from dat1 XX=x1+x3+x4 and YY=x2+x5. Can you please suggest me how can I do it in R. Actually, I have few more categories or types of variable.


title: What is a null conjunction analysis in an fMRI study?

text: I'm reading an article, The commonality of neural networks for verbal and visual short-term memory (Majerus et al., J Cogn Neurosci 2010 22(11): 2570), about brain imaging in which the results are analysed with multiple analyses. One of them is a null _conjunction analysis_. I've tried googling a bit but I didn't come up with something useful. I did find some results on 'conjunction analysis', but I'm not sure if that's the exact same thing. Can anyone explain to me what this analysis is (Goals, method...)?


title: Best analysis method for three interrelated longitudinal series

text: I have data which consist of three longitudinal series of financial data. The hypothesis is that two of the series are "caused" (in a loose sense) by the third. It is fine to investigate this as two separate hypotheses. All the series are about the relationship between two specific countries. There are data for about 20 quarters on one of the series, and about 80 months on the other two series. A little exploratory work reveals that none of the series are stationary, and that there is no apparent seasonality. So, how best to analysis these hypotheses? Can I simply try regressions with different lags, or should I do a more formal time series analysis, or should I do something else? Peter


title: How to draw a complete statistic from an incomplete statistical model?

text: My textbook gives an example, that normal distribution family $\\{N(0,\sigma^2):\sigma\in R^+\\}$ is not complete, but a complete statistic, $T_n=\sum_{i=1}^n X_i^2$, can still be constructed from samples $(X_1,\cdots,X_n)$. So under what (sufficient / necessary) conditions can we draw complete statistics from an incomplete model?